gs i jacobi

#include <iostream>
#include <vector>
#include <stdlib.h>
#include <time.h>
using namespace std;

vector<double> mnozenie_wierszy(vector<double> matylda, double liczba){
    for(int i=0; i<matylda.size(); i++)
        matylda[i] = matylda[i]*liczba;
    return matylda;
}

vector<double> odejmowanie_wierszy(vector<double> matylda, vector<double> wladyslaw){
    for(int i=0; i<matylda.size(); i++)
        matylda[i] = matylda[i]-wladyslaw[i];
    return matylda; //matylda - wladyslaw
}

vector<vector<double>> dodawanie(vector<vector<double>> a, vector<vector<double>> b)
{
    if(a.size()==b.size() && a.front().size() == b.front().size())
    {
        for(int i=0; i<a.size(); i++)
        {
            for(int j=0; j<a.front().size(); j++)
            {
                a[i][j]+=b[i][j];
            }
        }
        return a;
    }
}

int main()
{
    int t = 1;
    srand(time(NULL));
    int liczba_rownan = 3;
    double epsilon = 1e-10;
    int iter = 6;
    for(int iterator=0; iterator<t; iterator++)
    {
        //cin >> liczba_rownan;
        vector<vector<double> > macierzA;
        vector<vector<double> > macierzM;
        vector<vector<double> > macierzL;
        vector<vector<double> > macierzU;
        vector<vector<double> > macierzD;
        vector<vector<double> > macierzN;
        ////macierz L - od poczatku bez przekatnej
        ////macierz D - przekatna
        ////macierz U - od przekatnej do konca

        vector<double> x;
        vector<double> b;
        vector<double> przyblizenie, x1, x2, b_gs;

        for(int i=0; i<liczba_rownan; i++)
        {
            x.push_back((double)(rand()%999));
            //to przyblizenie poczatkowe moze byc inne
            //podobno zwykle daje sie wektor zer.
            przyblizenie.push_back(0);
            x1.push_back(0);
            x2.push_back(0);
            b_gs.push_back(0);
        }

        for(int i=0; i<liczba_rownan; i++)
        {
            vector<double> usagi;
            vector<double> marchewka;
            double btemp=0;
            for(int j=0; j<liczba_rownan; j++)
            {
                //chyba nie ma szans, zeby wyskoczyly tutaj jakies zera na glownej przekatnej.
                //dodam coś, derp. nigdy nie wiadomo z tymi randomami.
                double sailor_moon = (rand()%999) + 1e-10;
                if(i==j)
                {
                    sailor_moon= (rand()%99999)*liczba_rownan;
                    //diagonalnie dominujaca
                }
                usagi.push_back(sailor_moon);
                marchewka.push_back(0);
                btemp = btemp+sailor_moon * x[j];
            }

            b.push_back(btemp);
            macierzA.push_back(usagi);
            macierzM.push_back(marchewka);
            macierzD.push_back(marchewka);
            macierzU.push_back(marchewka);
            macierzL.push_back(marchewka);
            macierzN.push_back(marchewka);
            //teraz jest sytuacja troche inna niz w gausie bo mamy wektor b i wektor x
            //macierz a nie jest polaczona z wektorem b;
            //a do M tez odrazu wpisze bo potrzebuje zeby miala wymiary ustalone
            //w ogole wszedzie oprocz A sa same zera
        }

        //a moze dane z klawiatury? a może frytki do tego?
        //for(int i=0; i<liczba_rownan; i++)
        //{
        //  vector<double> usagi;
        //  for(int j=0; j<liczba_rownan+1; j++)
        //  {
        //      double sailor_moon;
        //      cin >> sailor_moon;
        //      usagi.push_back(sailor_moon);
        //  }
        //  macierz.push_back(usagi);
        //}

        //jacobi
        ////macierz L - od poczatku bez przekatnej
        ////macierz D - przekatna
        ////macierz U - od przekatnej do konca
        for(int i=0; i<liczba_rownan; i++)
        {
            for(int j=0; j<i; j++)
                macierzL[i][j]=macierzA[i][j];
        }

        for(int i=0; i<liczba_rownan; i++)
            macierzD[i][i]=macierzA[i][i];

        for(int i=0; i<liczba_rownan; i++)
        {
            for(int j=i+1; j<liczba_rownan; j++)
                macierzU[i][j]=macierzA[i][j];
        }

        // N = D^(-1)
        for(int i=0; i<liczba_rownan; i++)
            macierzN[i][i] = (1/macierzD[i][i]);

        for (int i=0; i<liczba_rownan; i++)
        {
            for (int j=0; j<liczba_rownan; j++)
            {
                if (i == j)
                    macierzM[i][j] = 0;
                else
                    macierzM[i][j] = - (macierzA[i][j] * macierzN[i][i]);
            }
        }

        for (int k=0; k<iter; k++) {
            for (int i=0; i<liczba_rownan; i++)
            {
                x2[i] = macierzN[i][i]*b[i];
                for (int j=0; j<liczba_rownan; j++)
                    x2[i] += macierzM[i][j]*x1[j];
            }
            for (int i=0; i<liczba_rownan; i++)
                x1[i] = x2[i];
        }

        printf("Wynik jacobi\n");
        for(int i=0; i<liczba_rownan; i++)
            x2[i]=abs(abs(x[i]-x1[i])/x[i]);

        for(int i = 0; i < liczba_rownan; i++)
        {
            //printf("x[%d] = %f\n", (i + 1), x1[i]) ;
            printf("blad wzgledny x[%d] = %f\n", (i + 1), x2[i]) ;
        }

        //to samo ale gauss-seidel czy jak tam.
        // Calculate D^-1 * b
        for (int i=0; i<liczba_rownan; i++)
            b_gs[i] = b[i] * macierzN[i][i];

        //Calculate D^-1 * L
        for (int i=0; i<liczba_rownan; i++)
            for (int j=0; j<i; j++)
                macierzL[i][j] *= macierzN[i][i];

        //Calculate D^-1 * U
        for (int i=0; i<liczba_rownan; i++)
            for (int j=i+1; j<liczba_rownan; j++)
                macierzU[i][j] *= macierzN[i][i];
        //tu tracimy rozwiazania z jacobiego
        for(int i=0; i<liczba_rownan; i++)
            x1[i]=0;

        for (int k=0; k<iter; k++)
            for (int i=0; i<liczba_rownan; i++)
            {
                x1[i] = b_gs[i];                       // x = D^-1*b -
                for (int j=0; j<i; j++)
                    x1[i] -= macierzL[i][j]*x1[j];    // D^-1*L * x -
                for (int j=i+1; j<liczba_rownan; j++)
                    x1[i] -= macierzU[i][j]*x1[j];    // D^-1*U * x
            }

            printf("Wynik gs\n");

            for(int i=0; i<liczba_rownan; i++)
                x2[i]=abs(abs(x[i]-x1[i])/x[i]);

            for (int i=0; i<liczba_rownan; i++)
            {
                printf("blad wzgledny x[%d] = %f\n", (i + 1), x2[i]) ;
            }
                //printf("x[%d] = %f\n", (i+1), x[i]);

        //double  eg_sr= 0;
        //double x_sr=0;
        //for(int i=0; i<x.size(); i++)
        //{
        //  eg_sr += abs(x_eg[i]);
        //  x_sr += abs(x[i]);
        //}
        //eg_sr = eg_sr/x.size();
        //x_sr=x_sr/x.size();

        //double bb= abs(eg_sr - x_sr);
        //double bw = bb/x_sr;


        //cout << "blad bezwzgledny " << bb << "\n";
        //cout << "blad wzgledny " << bw << "\n\n";
        //cout << bw << "\n";
        //if(bb == 0 || bw == 0)
        //{
        //  iterator = 0;
        //}
        //cout << liczba_rownan << "\t"
        //cout << liczba_obiegow; //<< "\n";
        //liczba_rownan+=5;
    }

    return 0;
}