DecomposeLU.hpp

Go to the documentation of this file.

#pragma once
#include "../vector_matrix/VectorMatrix.hpp"
#include <vector>
#include <utility>
#include <limits>
#include <type_traits>
#include <cmath>
 
namespace LinearAlgebra{
    template<typename T>
    class DecomposeLU {
        static_assert(std::is_floating_point_v<T>, 
          "T must be a floating-point type (float or double)");
    private:
        
        VectorMatrix<T> matrix;
        double detP;
        void elimination(int col);
        void initL() {
            int rows = matrix.getRows();
            int cols = matrix.getColumns();
            L = std::vector<std::vector<T>>(rows, std::vector<T>(cols, 0));
            for (int i = 0; i < rows; i++) {
                L[i][i] = 1;  
            }
            
        };
        void initP() {
            int n = matrix.getRows();
            P = std::vector<int>(n);
            for(int i=0; i < n; i++){
                P[i] = i; 
            }
        };
        std::vector<int> initInvP() const {
            int n = matrix.getRows();
            std::vector<int>invP(n,-1);
            for (int i = 0; i < n; i++)
            {
                int pi = P[i];
                if (pi < 0 || pi >= n) throw std::invalid_argument("invalid permutation");
                invP[pi] = i;
            }
            return invP;
        };
        int pivoting(int col);
        void forwardSubstitution(std::vector<T>& y, const std::vector<T>& b, int n) const;
        void backwardSubstitution(std::vector<T>& x, const std::vector<T>& y, int n) const;
        std::vector<std::vector<T>> L;
        std::vector<std::vector<T>> U;
        std::vector<int> P; //vector of swap
        static constexpr T eps = std::numeric_limits<T>::epsilon() * static_cast<T>(100);
        void decomposition();
    public:
 
        DecomposeLU(const VectorMatrix<T>& m) : matrix(m) { 
            detP = 1;
            decomposition(); 
        }
 
        DecomposeLU(VectorMatrix<T>&& m) : matrix(std::move(m)) { 
            detP = 1;
            decomposition(); 
        }
        T det() const;
        VectorMatrix<T> inv() const;
        const std::vector<std::vector<T>>& getU() const{ return U; }
        const std::vector<std::vector<T>>& getL() const{ return L; }
        const std::vector<int>& getP() const{ return P; }
    };
 
 
    template<typename T>
    int DecomposeLU<T>::pivoting(int col) {
        T pivotVal = std::abs(U[col][col]);
        int pivot = col;
        int rows = matrix.getRows();
        if (col < 0 || col >= rows) {
            throw std::runtime_error("Invalid column index");
        }
        for (int i = col+1; i < rows; i++) {
            auto &Ui = U[i];
            if(std::abs(Ui[col]) > pivotVal) {
                pivotVal = std::abs(Ui[col]);
                pivot = i;
            }
        }
        if (pivotVal <= eps) throw std::runtime_error("Matrix is singular");
        return pivot;
    }
 
 
    template<typename T>
    void DecomposeLU<T>::elimination(int col) {
        int rows = matrix.getRows();
        int cols = matrix.getColumns();
        
        if(std::abs(U[col][col]) <= eps) throw std::runtime_error("Pivot is zero");
        T pivot = U[col][col];
        for (int i = col+1; i < rows; i++)
        {
            auto &Ui = U[i];
            T& Li_col = L[i][col];
            Li_col=Ui[col]/pivot;
            const auto &Ucol = U[col]; 
 
            for (int j = col+1; j < cols; j++)
            {
                Ui[j] -= Li_col * Ucol[j];
            }
            Ui[col] = static_cast<T>(0);
        }
        
    }
 
 
    template<typename T>
    void DecomposeLU<T>::decomposition() {
        if (matrix.getRows() != matrix.getColumns())
        throw std::runtime_error("Matrix must be square for LU decomposition");
        int n = matrix.getRows();
        U = matrix.getMatrix();
        initL();
        initP();
        int swapCount = 0;
        for (int k = 0; k < n; ++k) {
            int pivot = pivoting(k);
         
            if(pivot != k) {
                using std::swap;
                swap(U[pivot],U[k]);
                swap(P[pivot],P[k]);
                if (k > 0) 
                    swap_ranges(L[pivot].begin(), L[pivot].begin() + k, L[k].begin());
                swapCount++;
            }
            elimination(k);
        }
        detP = swapCount % 2 == 0 ? 1 : -1;
    }
 
    template<typename T>
    T DecomposeLU<T>::det() const{
        const auto& U = getU();   
        int rows = U.size();
        T det = T(1);
        for (int i = 0; i < rows; i++) det*=U[i][i];
        
        return det*detP;
    }
 
    template<typename T>
    void DecomposeLU<T>::forwardSubstitution(std::vector<T>& y, const std::vector<T>& b, int n) const{
        for (int j = 0; j < n; ++j) {
            T sum = 0;
            auto& Lj = L[j];
            auto& Ljj = L[j][j];
            for (int k = 0; k < j; ++k) sum += Lj[k] * y[k];
            y[j] = (b[j] - sum)/Ljj;
        }
    }
 
    template<typename T>
    void DecomposeLU<T>::backwardSubstitution(std::vector<T>& x, const std::vector<T>& y, int n) const{
        for (int k = n-1; k >= 0; --k){
            T sum = 0;
            auto& Uk = U[k];
            auto& Ukk = U[k][k];
            for (int s = k+1; s < n; ++s) sum += Uk[s] * x[s];
            x[k] = (y[k] - sum)/Ukk;
        }
    }
 
    template<typename T>
    VectorMatrix<T> DecomposeLU<T>::inv() const{
        int n = matrix.getRows();
        std::vector<std::vector<T>> X(n, std::vector<T>(n));
        std::vector<int> invP = initInvP();
        std::vector<T> b(n), y(n), x(n);
        int prev = invP[0];
        b[prev] = T(1);
        for(int i = 0; i < n; ++i) {
            int curr = invP[i];
 
            if (i > 0) {
                b[prev] = T(0);
                b[curr] = T(1);
            }
            
            forwardSubstitution(y, b, n);
            backwardSubstitution(x, y, n);
 
            for (int k = 0; k < n; ++k) {
                auto& Xki = X[k][i];
                Xki = x[k];
            }
 
            prev = curr;
        }
 
        return VectorMatrix<T>(std::move(X));
    }
}