LinearSystem.cpp

#include "LinearSystem.h"
#include "cassert"

    LinearSystem::LinearSystem(const Matrix& inputA,const Vector& inputB) : A(inputA), b(inputB), n(inputB.size())
    {
        assert(A.rows() == A.cols());
        assert(A.rows() == b.size());
    }
    Vector LinearSystem::solve()
    {
        Matrix U = A;
        Vector rhs = b;
        Vector x(n);
        //forward elimination
        for (int k = 1; k <= n; k++)
        {
            double pivot = U(k, k);
            assert(pivot != 0);

            for (int i = k + 1; i <= n; i ++)
            {
                double factor = U(i, k)/ pivot;
                for (int j = k; j <= n; j++)
                {
                    U(i, j) -= factor * U(k, j);
                }
                rhs(i) -= factor * rhs(k);
            }
        }
        //back substitution
        for (int i = n; i >= 1; i--)
        {
            double sum = rhs(i);
            for (int j = i + 1; j <= n; j++)
            {
                sum -= U(i, j) * x(j);
            }
            x(i) = sum / U(i, i);
        }
        return x;
    }