Metoda Choleskiego. Rozkład macierzy L i U Gauss

/*

Autor: Dawid Mocek

2014-01-19

Testowane pod MS VS Express 2013

Napisać program rozwiązujący układ równań AX=B metodą Choleskiego
Program ma zawierać:
a) funk. czytającą macierz A i wektor B z pliku oraz rozmiar n macierzy A i wektorów B i X z kalawiatury
b) funk. tworzącą macierze Li U (metodą Gaussa), oraz sprawdzającą warunek utworzenia tych elementów (kontrola dzielenia przez 0)
c) funkcję rozwiązującą dwa układy równań z macierzami trójkątnymi
d) funkcję generująca raport do pliku....


Obsługa:

Plik tekstowy z zestawem ma mieć  tak poukładane wartość:

--- Zesta1.txt /Start/ ----
10 1 1 1 2 15
3 20 4 2 1 30
5 1 40 9 5 60
3 0 1 10 1 15
6 2 1 1 20 30
---- /End/ ------

Pierwsze 5 kolumn to macierz A a ostatnia Kolumna to wektor B

czyli dla takiego zestawu przy uruchomieniu programu(tak został zaprojketowany) podajemy liczbę 5 i program wygeneruje raport - a jak nie wygeneruje plik raportu utworzyc recznie

plik wynikowy-raportu oraz zestawu nalezy ustawić w funkcji main()

*/

#include <iostream>
#include <string>
#include <vector>
#include <fstream>
#include <iomanip>

// Dla funkcji GenerujRaport() - gdziekolwiek indziej stosować z rozwagą !
#define for_loop(a, b, c) for (int a = b; a < c; a++)
#define for_i_loop for_loop(i, 0, n)
#define for_j_loop for_loop(j, 0, n)

using namespace std;

// Czytelniejszy kod
typedef vector<vector<double>> Macierz;
typedef vector<double> Wektor;


/*
a) Funkcja czytjąca macierz A i wektor B z pliku oraz rozmiar n macierzy A i wektorów B i X z klawiatury

Tu modyfikujemy domyślne Macierze A i B utowrzone w main() także potrzebne są &
Zwraca: rozmiar 'N' macierzy - potrzebne dla pętli w następych funkcjach
*/
int Czytaj(string plik_zestawu, Macierz &A, Wektor& B)
{

    int n = 0;
    cout << "Podaj rozmiar 'N' macierzy A i wektorow B i X: ";
    cin >> n;

    // Poszerzamy wektor B i macierz A do rozmiaru N;
    A.resize(n);
    for (int i = 0; i < n; i++)
        A[i].resize(n);

    B.resize(n);

    // Czytamy wartości macierzy i wektoru z pliku
    ifstream ifs_zestaw = ifstream(plik_zestawu, ios::in);

    for (int i = 0; i < n; i++)
    {
        // Przesuwamy o +1 dla pobrania wartości Wektora - ostatnie liczby w liniach w pliku zestawu
        for(int j = 0; j < n + 1; ++j)
        {

            if (j < n)
                // Ładujemy wartości do macierzy
                ifs_zestaw >> A[i][j];
            else
                // Ładujemy wartości wektora
                ifs_zestaw >> B[i];
        }
    }

    ifs_zestaw.close();

    // Zwaracamy n-kę
    return n;
}

/*
b) Funkcja tworząca macierze L i U  (metodą Gaussa) oraz sprawdzającą warunek utworzenia tych elementów (kontrola operacji dzielenia przez 0)

Wyznaczamy macierze L i Y (utworzone w main()) które w tej funkcji musimy nadpisać - &
Dbamy o nie nadpisanie czegokolwiek w Macierzy A - const

Funkcja zwraca TRUE jeśli wszystko poszło dobrze, FALSE jeśli doszło do dzielenia przez zero.
*/
bool LU(int n, const Macierz A, Macierz &L, Macierz &U)
{
    // Nasze kochane macierze L i U
    L.resize(n);
    U.resize(n);

    for (int i = 0; i < n; i++)
    {
        L[i].resize(n);
        U[i].resize(n);

        for (int j = 0; j < n; j++)
        {
            U[i][j] = A[i][j];
            L[i][j] = 0;
        }
    }

    // Jedynki po przekątnej
    for (int i = 0; i < n; i++)
        L[i][i] = 1;

    for (int k = 0; k < n; k++)
    {
        for (int i = k + 1; i < n; i++)
        {
            // Sprawdź dzielenie
            if (U[k][k] == 0)
                return false;

            // Obliczamy marzcierz L
            L[i][k] = U[i][k] / U[k][k];

            for (int j = k + 1; j < n; j++)
            {
                // Obliczamy macierzy U
                U[i][j] -= ((U[i][k] * U[k][j]) / U[k][k]);
            }

            // Zerowanie U
            U[i][k] = 0;
        }
    }

    return true;
}

/*
c) Funkcja rozwiązuje dwa układy równań z macierzami trjókątnymi

Dbamy o nie nadpisywanie czegokolwiek w L, U i B - const

*/
void XY(int n, const Macierz L, const Macierz U, const Wektor B, Wektor& X, Wektor& Y)
{
    // X i Y poszerzamy
    X.resize(n);
    Y.resize(n);

    // tmp - przeochwuje przy wyznaczaniu X i Y sumę iloczynów L[][] z Y[] oraz U[][] z X[]
    double tmp;

    for (int i = 0; i < n; i++)
    {
        tmp = 0;
        for (int j = 0; j < n; j++)
        {
            tmp += (L[i][j] * Y[j]);
        }
        Y[i] = B[i] - tmp;
    }


    X[n - 1] = Y[n - 1] / U[n - 1][n - 1];

    for (int i = (n - 2); i >= 0; i--)
    {
            tmp = 0;

        for (int k = i; k < n; k++)
        {
            tmp += (U[i][k] * X[k]);
        }

        X[i] = (Y[i] - tmp) / U[i][i];
    }
}
/*
d) Funkcja generująca raport do pliku. Raport ma zawierać dane wejściowe (macierz A, wektor B), dane pośrednie (macierze L i U)
oraz wyniki obliczeń wektorów Y i X. Obydwa wektory mają mieć postać wykładniczą(mantysa i cecha), a mantysa ma mieć 10 cyfr znaczących.
*/
void GenerujRaport(string plik_raportu, int n, const Macierz A, const Wektor B, const Macierz L, const Macierz U, const Wektor X, const Wektor Y)
{
    ofstream ofs_raport;
    ofs_raport.open(plik_raportu);


    ofs_raport << "--> Macierz A <--" << endl;
    for_i_loop
    {
        for_j_loop { ofs_raport << A[i][j] << " "; }

        ofs_raport << endl;
    }
    ofs_raport << endl;


    ofs_raport << "--> Wektor B <--" << endl;
    for_i_loop { ofs_raport << B[i] << endl; }
    ofs_raport << endl;

    ofs_raport << "--> Macierz L <--" << endl;
    for_i_loop
    {
        for_j_loop { ofs_raport << L[i][j] << " "; }

        ofs_raport << endl;
    }
    ofs_raport << endl;

    ofs_raport << "--> Macierz U <--" << endl;
    for_i_loop
    {
        for_j_loop { ofs_raport << U[i][j] << " "; }
        ofs_raport << endl;
    }
    ofs_raport << endl;

    // 10 cyfr znaczących
    ofs_raport.precision(10);

    // Ustawmy notacje wykładniczą dla wektorów
    //ofs_raport.setf(ios::scientific);

    ofs_raport << "--> Wektor X <--" << endl;
    if (X.empty())
    {

        ofs_raport << "Nie istnieje z powodu bledu dzielenia" << endl;
    }
    else
    {
        for_i_loop{ ofs_raport << fixed << X[i] << endl; }
        ofs_raport << endl;
    }


    ofs_raport << "--> Wektor Y <--" << endl;
    if (Y.empty())
        ofs_raport << "Nie istnieje z powodu bledu dzielenia" << endl;
    else
    {
        for_i_loop{ ofs_raport << fixed  << Y[i] << endl; }
        ofs_raport << endl;
    }

    ofs_raport.close();

}

int main(void)
{

    // Plik z zestawem
    string plik_zestawu = "D:\\Choleski\\Zestawy\\Zestaw4.txt";

    // Plik raportu - trzeba go ręcznie utowrzyć z powodu praw
    string plik_raportu = "D:\\Choleski\\Raporty\\Raport_zestaw4.txt";

    // Nasza Macierz i Wektor
    Macierz A;
    Wektor B;

    // Czytamy z pliku
    int N = Czytaj(plik_zestawu, A, B);

    Macierz L;
    Macierz U;

    bool blad = LU(N, A, L, U);

    //  Wektory X i Y ustawiamy na puste(0) - na wypadek zwrócenia FALSE przez funkcje LU.
    //  Dzięki temu możemy łatwo srawdzić w funkcji GenerujRaport
    //  czy wektory dalej są puste i zapisać  do pliku raportu bląd dzielenia przez zero.
    Wektor X(0);
    Wektor Y(0);

    // Jesli wystąpił błąd dzielenia przez zero - XY się nie wykona
    if (blad)
    {
        XY(N, L, U, B, X, Y);
    }

    GenerujRaport(plik_raportu, N, A, B, L, U, X, Y);

    return 0;
}