#include <iostream>
#include <cstdlib>


using namespace std;

int n = 5;



double get_random(int fMin, int fMax){
    double f = (double)rand() / RAND_MAX;
    double losowa = fMin + f * (fMax - fMin);
    return losowa;
}


// fill matrix with random numbers between fmin and fmax
void fill_matrix(){
    double matrix_a[n][n];
    for (int i=0; i<n; i++){
        for (int j=0; j<n; j++){
            matrix_a[i][j] = 0;
        }
    }

    int fMin = 1;
    int fMax = 5;

    for (int i=0; i<n; i++){
        for (int j=0; j<n; j++){
            double losowa = get_random(fMin,fMax);
            matrix_a[i][i] = 8*losowa;
            try{
                losowa = get_random(fMin,fMax);
                matrix_a[i][i-1] = losowa;

                losowa = get_random(fMin,fMax);
                matrix_a[i][i+1] = losowa;
            }
            catch(...){
            }
        }
    }
}

//print matrix
void print_matrix(){
    double matrix_a[n][n];
    for (int i=0; i<n; i++){
        for (int j=0; j<n; j++){
            cout << matrix_a[i][j] << "\t\t";
            if (j == (n-1) ){
                cout << endl;
            }
        }
    }
}


//fill vector
void fill_vector(){
    double matrix_b[n];
    int fMin = 10;
    int fMax = 40;
    for (int w=0; w<n; w++){
            double losowa = get_random(fMin,fMax);
            matrix_b[w] = losowa;
    }
}

//print vector
void print_vector(){
    double matrix_b[n];
    for (int i=0; i<n; i++){
        cout << matrix_b[i] << endl;
    }
}


//solve system equation - SOR method
void solve(){
    double matrix_a [n][n];
    double matrix_b [n];

    double SOR_param = 1.0;
    double max_residuum = 0;
    double r0 = 1E-2;
    bool active_comp = true;
    int iter = 0;
    int maxiter = 100;

    // initialize array

    double x_old[n] = {0};
    double x[n] = {0};
    double residuum[n] = {0};

    //fill x[n] with zeros
    for (int k=0; k<=n; k++){
        x[k] = 0;
        x_old[k] = 0;
        residuum[k] = 0;
    }

    fill_matrix();
    fill_vector();
    cout << endl;
    print_matrix();
    cout << endl;
    print_vector();
    cout << endl;

        cout<< "Startowe" << endl;
        cout << "x[" << 0 << "]: " << x[0] << endl;
        for (int l =1;l<n;l++){
            cout << "x[" << l << "]: " << x[l] << endl;}
        cout << "x[" << n << "]: " << x[n] << endl;
        cout <<"Koniec startowych" << endl;

    while (active_comp){


        x[0] = ( matrix_b[0] - x[1] * matrix_a[0][1] )/matrix_a[0][0];
        x[n] = ( matrix_b[n] - x[n-1] * matrix_a[n][n-1] )/matrix_a[n][n];

        //solve
        for (int i=1; i<n; i++){
            //for (int j=1; j<n-1; j++){
                iter += 1;
                //try{
                    x[i] =  ( ( matrix_b[i] - x[i-1]*matrix_a[i][i-1] - x[i+1]*matrix_a[i][i+1] )/matrix_a[i][i] );// + ( (1.0 - SOR_param) * x_old[i] ) );

                    /*
                    residuum[i] = abs((x[i] - x_old[i])/x[i]);

                    if (residuum[i] > max_residuum){
                        max_residuum = residuum[i];
                    }
                    }
                    */

                    if (iter == maxiter){
                        cout << "End of computation,max iter, iter: " << iter << endl;
                        active_comp = false;
                    }

                    /*
                    if (max_residuum < r0){
                        cout << "End of computation, iter: " << iter << endl;
                        active_comp = false;
                    }
                    */

                    else{
                        x_old[i] = x[i];
                        //cout << "Iter: " << iter << endl;
                    }
                //}
                //catch(...){}


            //}
        }
    }
}


int print_solution(){
    double x[n];
    cout << "Sollution: ";
    for (int i=1; i<n; i++){
        cout << x[i] << " ";
    }
}



int main(){

    fill_matrix();
    print_matrix();

    cout << endl << endl;

    fill_vector();
    print_vector();

    cout <<endl<<endl;

    srand (time( 0 ));
    solve();
    print_solution();
    cout << endl;


    return 0;
}