For_Ilyusha

#include "SlaeSolver.h"

SlaeSolver::SlaeSolver(const ExtendedMatrix & other)
{
    matrix = new ExtendedMatrix(other);
    roots = other.getColumn();
    N = matrix->getRows();
}

// -------------------- Метод Гаусса с выбором главного элемента --------------------
vector<double> SlaeSolver::gaussMethod()
{
    ExtendedMatrix * matrixCopy = new ExtendedMatrix(*matrix);
    double tmp = 0.0;

    // Прямой ход метода Гаусса
    for (int n = 0; n != N; ++n)
    {
        tmp = matrix->at(n,n);
        if (fabs(tmp) < 0.001)
        {
            double max = fabs(matrix->at(n + 1, n));
            int iMax = n + 1;
            for (int i = n + 2; i != N; ++i)
            {
                if (fabs(matrix->at(i, n)) > max)
                {
                    max = fabs(matrix->at(i, n));
                    iMax = i;
                }
            }
            if (fabs(max) < 0.001)
                throw runtime_error("Error: Cant solve this system!");
            else
                matrix->swap(n, iMax);
        }
        for (int i = N; i >= n; --i)
            matrix->at(n,i) /= tmp;

        for (int i = n + 1; i != N; ++i)
        {
            tmp = matrix->at(i,n);
            for (int j = N; j >= n; --j)
                matrix->at(i, j) -= tmp * matrix->at(n, j);
        }
    }

    // Обратный ход метода Гаусса
    roots.at(N - 1) = matrix->at(N - 1, N);
    for (int i = N - 2; i >= 0; --i)
    {
        roots.at(i) = matrix->at(i, N);
        for (int j = i + 1; j != N; ++j)
            roots.at(i) -= matrix->at(i, j) * roots.at(j);
    }

    matrix = matrixCopy;
    return roots;
}

// -------------------- L-U разложение --------------------
vector<double> SlaeSolver::LuDecompositionMethod()
{
    if (L == NULL && U == NULL) // Если L и U матрицы ещё не рассчитаны
    {
        ExtendedMatrix * matrixCopy = new ExtendedMatrix(*matrix);
        double tmp = 0.0;

        L = new SquareMatrix(N);

        // Прямой ход метода Гаусса
        for (int n = 0; n != N; ++n)
        {
            tmp = matrix->at(n, n);
            if (fabs(tmp) < 0.001)
            {
                double max = fabs(matrix->at(n + 1, n));
                int iMax = n + 1;
                for (int i = n + 2; i != N; ++i)
                {
                    if (fabs(matrix->at(i, n)) > max)
                    {
                        max = fabs(matrix->at(i, n));
                        iMax = i;
                    }
                }
                if (fabs(max) < 0.001)
                    throw runtime_error("Error: Cant solve this system!");
                else
                    matrix->swap(n, iMax);
            }

            L->at(n, n) = tmp; // Заполняем матрицу L (элементы главной диагонали)
            for (int i = N; i >= n; --i)
                matrix->at(n, i) /= tmp;

            for (int i = n + 1; i != N; ++i)
            {
                tmp = matrix->at(i, n);
                L->at(i, n) = tmp; // Заполняем матрицу L (элменты под главной диагональю)
                for (int j = N; j >= n; --j)
                    matrix->at(i, j) -= tmp * matrix->at(n, j);
            }
        }

        // Берём матрицу U как матрицу с единицами на главной диагонали, получившуюся после
        // выполнения прямого хода метода Гаусса
        U = new SquareMatrix(matrix->getMatrix());
        for (int i = 0; i != N; ++i)
            U->at(i, i) = 1.0;

        matrix = matrixCopy;
    }

    // Находем матрицу-столбец Y по формуле
    vector<double> Y(N);
    Y.at(0) = matrix->getColumn().at(0) / L->at(0, 0);
    for (int i = 1; i != N; ++i)
    {
        Y.at(i) = matrix->getColumn().at(i);
        for (int j = 0; j != i; ++j)
            Y.at(i) -= L->at(i, j) * Y.at(j);
        Y.at(i) /= L->at(i, i);
    }

    // Находим корни по формуле
    roots.at(N - 1) = Y.at(N - 1);
    for (int i = N - 2; i != -1; i--)
    {
        roots.at(i) = Y.at(i);
        for (int j = i + 1; j != N; ++j)
            roots.at(i) -= U->at(i, j) * roots.at(j);
    }

    return roots;
}

// Задать столбец свободных членов
void SlaeSolver::setColumn(const vector<double> column)
{
    this->matrix->setColumn(column);
    roots = column;
}

// Задать СЛАУ целиком
void SlaeSolver::setSlae(const ExtendedMatrix & extendedMatrix)
{
    *(this->matrix) = extendedMatrix;
    L->clear();
    U->clear();
    roots.clear();
    roots = extendedMatrix.getColumn();
}

// -------------------- Метод прогонки --------------------
vector<double> SlaeSolver::sweepMethod()
{
    TridiagonalMatrix * m = new TridiagonalMatrix(this->matrix->getMatrix());

    bool notEnoughPresicionWarning = false;

    L->clear();
    U->clear();

    int n = this->matrix->getRows();
    vector<double> alphaArr;
    vector<double> betaArr;
    double alpha = 0.0, beta = 0.0, gamma = 0.0, a = 0.0, b = 0.0, c = 0.0, d = 0.0;


    // Первая итерация прямого хода
    b = m->getMainDiagonal().at(0);
    c = m->getAboveDiagonal().at(0);
    d = matrix->getColumn().at(0);

    if (!(fabs(b) >= fabs(a) + fabs(c) && fabs(b) > fabs(a)))
        notEnoughPresicionWarning = true;

    gamma = b;
    beta = d / gamma;
    alpha = -c / gamma;
    if (gamma == 0)
        throw runtime_error("Can't solve this! Gamma is 0!");

    alphaArr.push_back(alpha);
    betaArr.push_back(beta);

    // Прямой ход
    for (int i = 1; i < n - 1; ++i)
    {
        a = m->getBelowDiagonal().at(i - 1);
        b = m->getMainDiagonal().at(i);
        c = m->getAboveDiagonal().at(i);
        d = matrix->getColumn().at(i);

        if (!(fabs(b) >= fabs(a) + fabs(c) && fabs(b) > fabs(a)))
            notEnoughPresicionWarning = true;

        gamma = b + a * alphaArr.at(i - 1);
        beta = (d - a * betaArr.at(i - 1)) / gamma;
        alpha = -c / gamma;
        if (gamma == 0)
            throw runtime_error("Can't solve this! Gamma is 0!");

        alphaArr.push_back(alpha);
        betaArr.push_back(beta);
    }

    // Последняя итерация прямого хода
    a = m->getBelowDiagonal().at(n - 2);
    b = m->getMainDiagonal().at(n - 1);
    d = matrix->getColumn().at(n - 1);

    if (!(fabs(b) >= fabs(a) + fabs(c) && fabs(b) > fabs(a)))
        notEnoughPresicionWarning = true;

    gamma = b + a * alphaArr.at(n - 2);
    beta = (d - a * betaArr.at(n - 2)) / gamma;
    if (gamma == 0)
        throw runtime_error("Can't solve this! Gamma is 0!");

    roots.at(n - 1) = beta;


    // Обратный ход
    for (int i = n - 2; i > -1; --i)
    {
        alpha = alphaArr.back();
        alphaArr.pop_back();

        beta = betaArr.back();
        betaArr.pop_back();

        roots.at(i) = alpha * roots.at(i + 1) + beta;
    }

    if (notEnoughPresicionWarning)
        cout << "WARNING! Loss of accuracy is possible using sweep method!" << endl << endl;

    return roots;
}

// -------------------- Тестирование --------------------
void SlaeSolver::showRoots()
{
    cout << "Roots:" << endl;
    for (int i = 0; i != roots.size(); ++i)
        cout << "\t" << roots.at(i) << endl;
}

void SlaeSolver::test()
{
    cout << "Slae matrix:" << endl;
    matrix->show();
    cout << endl;
    showRoots();
    cout << endl << "Checking roots:" << endl;
    Tester * tester = new Tester(this->matrix, roots);
    tester->test();
}