BandMatrix.hpp

#pragma once
#include <NDEVR/Buffer.h>
#include <NDEVR/Vertex.h>
namespace NDEVR
{
    template<class t_type>
    class BandMatrix
    {
    private:
        Buffer<Buffer<t_type>> m_upper;  // upper band
        Buffer<Buffer<t_type>> m_lower;  // lower band
    public:
        BandMatrix() {}; // constructor
        ~BandMatrix() {}; // matrix dimension
        uint04 upperCount() const
        {
            return m_upper.size() - 1;
        }
        uint04 lowerCount() const
        {
            return m_lower.size() - 1;
        }
        BandMatrix(uint04 dim, uint04 n_u, uint04 n_l)
        {
            resize(dim, n_u, n_l);
        }
        void resize(uint04 dim, uint04 n_u, uint04 n_l)
        {
            lib_assert(dim > 0U, "dim > 0");
            m_upper.setSize(n_u + 1);
            m_lower.setSize(n_l + 1);
            for (uint04 i = 0; i < m_upper.size(); i++)
                m_upper[i].setSize(dim);
            for (uint04 i = 0; i < m_lower.size(); i++)
                m_lower[i].setSize(dim);
        }
        uint04 dim() const
        {
            if (m_upper.size() > 0)
                return m_upper[0].size();
            else
                return 0;
        }
        // defines the new operator (), so that we can access the elements
        // by A(i,j), index going from i=0,...,dim()-1
        t_type& operator() (uint04 i, uint04 j)
        {
            if (j >= i)
                return m_upper[j - i][i];
            else
                return m_lower[i - j][i];
        }
        t_type operator() (uint04 i, uint04 j) const
        {
            if (j >= i) 
                return m_upper[j - i][i];
            else        
                return m_lower[i - j][i];
        }
        // second diag (used in LU decomposition), saved in m_lower
        t_type savedDiag(uint04 i) const
        {
            lib_assert(i < dim(), "Bad band location");
            return m_lower[0][i];
        }
        t_type& savedDiag(uint04 i)
        {
            lib_assert(i < dim(), "Bad band location");
            return m_lower[0][i];
        }
 
        // LR-Decomposition of a band matrix
        void luDecompose()
        {
            for (uint04 i = 0; i < dim(); i++) 
            {
                lib_assert((*this)(i, i) != 0.0, "Bad decompision");
                savedDiag(i) = 1.0 / (*this)(i, i);
                uint04 j_min = getMax(0U, i - lowerCount());
                uint04 j_max = getMin(dim() - 1, i + upperCount());
                for (uint04 j = j_min; j <= j_max; j++)
                    (*this)(i, j) *= savedDiag(i);
                (*this)(i, i) = 1.0;          // prevents rounding errors
            }
 
            // Gauss LR-Decomposition
            for (uint04 k = 0; k < dim(); k++) 
            {
                uint04 i_max = getMin(dim() - 1, k + lowerCount());  // lowerCount not a mistake!
                for (uint04 i = k + 1; i <= i_max; i++) 
                {
                    lib_assert((*this)(k, k) != 0.0, "get minimum");
                    t_type x = -(*this)(i, k) / (*this)(k, k);
                    (*this)(i, k) = -x;                      // assembly part of L
                    uint04 j_max = getMin(dim() - 1, k + upperCount());
                    for (uint04 j = k + 1; j <= j_max; j++)
                    {
                        // assembly part of R
                        (*this)(i, j) = (*this)(i, j) + x * (*this)(k, j);
                    }
                }
            }
        }
        // solves Ly=b
        template<class t_l_type>
        Buffer<t_l_type> lSolve(const Buffer<t_l_type>& b) const
        {
            lib_assert(dim() == b.size(), "Bad size");
            Buffer<t_l_type> x;
            x.setSize(dim());
            for (uint04 i = 0; i < dim(); i++)
            {
                t_l_type sum(0);
                uint04 j_start = getMax(0U, i - lowerCount());
                for (uint04 j = j_start; j < i; j++) 
                    sum += (*this)(i, j) * x[j];
                x[i] = (b[i] * savedDiag(i)) - sum;
            }
            return x;
        }
        // solves Rx=y
        template<class t_r_type>
        Buffer<t_r_type> rSolve(const Buffer<t_r_type>& b) const
        {
            lib_assert(dim() == b.size(), "dimensions not equal");
            Buffer<t_r_type> x;
            x.setSize(dim());
            for (uint04 i = dim() - 1; IsValid(i); i--) 
            {
                t_r_type sum(0);
                uint04 j_stop = getMin(dim() - 1, i + upperCount());
                for (uint04 j = i + 1; j <= j_stop; j++) 
                    sum += (*this)(i, j) * x[j];
                x[i] = (b[i] - sum) / (*this)(i, i);
            }
            return x;
        }
        template<class t_l_type>
        Buffer<t_l_type> luSolve(const Buffer<t_l_type>& b, bool is_lu_decomposed = false)
        {
            lib_assert(dim() == b.size(), "Bad size");
            if (!is_lu_decomposed)
                luDecompose();
            Buffer<t_l_type> y = lSolve(b);
            Buffer<t_l_type> x = rSolve(y);
            return x;
        }
    };
}