jacobi smierdzi.

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

int main()
{
    int t = 1;
    srand(time(NULL));
    int liczba_rownan = 100;
    double epsilon = 1e-10;
    int iter = 10000;
    cout.precision(26);
    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]);
            }
        }

        //tu wlasciwe obliczanie
        double blad=0;
        int jacobi=0; //k to iteracje ktorych szukamy, albo droidy
        for (jacobi=0; jacobi<iter; jacobi++) {
            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];

            blad=0;
            for(int i=0; i<liczba_rownan; i++)
            {
                double usagi=0;
                for(int j=0; j<liczba_rownan; j++)
                {
                    usagi+=macierzA[i][j]*x1[j];
                }
                usagi=abs(usagi-b[i]);
                blad+=usagi;
            }

            if(blad<epsilon)
                break;
        }

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

        double blad_sredni_j=blad/liczba_rownan;

        //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]) ;
        //  cout << "blad wzgledny x[" << i+1 << "] = " << x2[i] << "\n";
        //}
        cout << "blad sredni wzgledny = " << blad_sredni_j << "\n";
        cout << "iteracje = " << jacobi << "\n";

        //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;

        int gauss_s =0;
        for (gauss_s=0;gauss_s<iter; gauss_s++)
        {
            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
            }

            blad=0;
            for(int i=0; i<liczba_rownan; i++)
            {
                double usagi=0;
                for(int j=0; j<liczba_rownan; j++)
                {
                    usagi+=macierzA[i][j]*x1[j];
                }
                usagi=abs(usagi-b[i]);
                blad+=usagi;
            }

            if(blad<epsilon)
                break;
        }

            printf("Wynik gs\n");

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

            double blad_sredni_gs=blad/liczba_rownan;
            //for(int i=0; i<liczba_rownan; i++)
            //  blad_sredni_gs+=x2[i];

            //blad_sredni_gs=blad_sredni_gs/liczba_rownan;


            //for (int i=0; i<liczba_rownan; i++)
            //{
            //  //printf("blad wzgledny x[%d] = %f\n", (i + 1), x2[i]) ;
            //  cout << "blad wzgledny x[" << i+1 << "] = " << x2[i] << "\n";
            //}
            cout << "iteracje = " << gauss_s << "\n";
            cout << "blad sredni wzgledny = " << blad_sredni_gs << "\n";
                //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;
}