Metoda Choleskiego. Rozkład macierzy L i U Doolittle

/*

Autor: Dawid Mocek

2014-01-05

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ą Doolittle), 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>
using namespace std;
// Dla czytelniejszego kodu
typedef vector<vector<double>> Macierz;
typedef vector<double> Wektor;
/*
a) Funkcja czytają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 utworzone w main()
Zwraca: rozmiar 'N' macierzy - potrzebne dla pętli w następnych 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 wektora 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 macierz
                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ą Dolittle) oraz sprawdzającą warunek utworzenia tych
elementów (kontrola operacji dzielenia przez 0)
Wyznaczamy macierze L i U (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);
    for (int i = 0; i < n; i++)
        L[i].resize(n);
    U.resize(n);
    for (int i = 0; i < n; i++)
        U[i].resize(n);
    // tmp - suma iloczynów L[][] i U[][]
    double tmp = 0;
    for (int i = 0; i < n; i++)
        L[i][i] = 1;
    for (int j = 0; j < n; j++)
        U[0][j] = A[0][j];
    // Dzielenie przez zero
    if (U[0][0] == 0)
    {
        return false;
    }
    for (int i = 1; i < n; i++)
        L[i][0] = (A[i][0] / U[0][0]);
    for (int i = 1; i < n; i++)
    {
        for (int j = i; j < n; j++)
        {
            tmp = 0;
            for (int k = 0; k < i; k++)
                tmp = tmp + (L[i][k] * U[k][j]);
            U[i][j] = A[i][j] - tmp;
            // Dzielenie przez zero
            if (U[i][i] == 0)
            {
                return false;
            }
            if (i != j)
            {
                tmp = 0;
                for (int k = 0; k < i; k++)
                    tmp = tmp + L[j][k] * U[k][i];
                L[j][i] = (A[j][i] - tmp) / U[i][i];
            }
        }
    }
    return true;
}
/*
c) Funkcja rozwiązuje dwa układy równań z macierzami trójką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 - przechowuje przy wyznaczaniu X i Y sumę iloczynów L[][] z Y[] oraz U[][] z X[]
    double tmp = 0;
    Y[0] = B[0];
    for (int i = 0; i < n; i++)
    {
        tmp = 0;
        for (int k = 0; k < i; k++)
        {
            tmp = tmp + (L[i][k] * Y[k]);
        }
        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 = 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 (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            ofs_raport << A[i][j] << " ";
        }
        ofs_raport << endl;
    }
    ofs_raport << endl;
    ofs_raport << "--> Wektor B <--" << endl;
    for (int i = 0; i < N; i++)
    {
        ofs_raport << B[i] << endl;
    }
    ofs_raport << endl;
    ofs_raport << "--> Macierz L <--" << endl;
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            ofs_raport << L[i][j] << " ";
        }
        ofs_raport << endl;
    }
    ofs_raport << endl;
    ofs_raport << "--> Macierz U <--" << endl;
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            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 (int i = 0; i < N; i++)
        {
            ofs_raport << 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 (int i = 0; i < N; i++)
        {
            ofs_raport << Y[i] << endl;
        }
        ofs_raport << endl;
    }
    ofs_raport.close();
}
int main(void)
{
    // Plik z zestawem
    string plik_zestawu = "D:\\Choleski\\Zestawy\\Zestaw1.txt";
    // Plik raportu - trzeba go ręcznie utworzyć z powodu praw
    string plik_raportu = "D:\\Choleski\\Raporty\\Raport_zestaw1.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 sprawdzić w funkcji GenerujRaport
    // czy wektory dalej są puste i zapisać do pliku raportu błąd dzielenia przez zero.
    Wektor X(0);
    Wektor Y(0);
    // Jeśli 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;
}