TridiagonalMatrixSolver.h

Go to the documentation of this file.

#ifndef TRIDIAGONAL_MATRIX_SOLVER_H
#define TRIDIAGONAL_MATRIX_SOLVER_H
 
#include <functional>
#include <memory>
#include <valarray>
#include <vector>
 
#include "GenericSolver.h"
#include "SimulationClassBase.h"
 
template <typename T>
class TridiagonalMatrixSolver : public SolverInterface<T>
{
    private:
        std::vector<std::reference_wrapper<T>> m_lower_sub_diag_matrix_elements{};
        std::vector<std::reference_wrapper<T>> m_diag_matrix_elements{};
        std::vector<std::reference_wrapper<T>> m_upper_sub_diag_matrix_elements{};
        // sub-diagonals are 1 elem smaller than diagonal, but all arrays are expected to be of same
        // size. The reference for this extra filler element is set to this unused element.
        T m_unused_element{0};
 
        void getReferenceToElements();
 
    public:
        TridiagonalMatrixSolver(SimulationClassBase<T>* sim_ptr,
                                std::shared_ptr<MatrixManager<T>> matrix_manager_ptr);
        void solveMatrixSystem(std::string_view field_name) override;
};
 
template <typename T>
TridiagonalMatrixSolver<T>::TridiagonalMatrixSolver(
    SimulationClassBase<T>* sim_ptr, std::shared_ptr<MatrixManager<T>> matrix_manager_ptr)
    : SolverInterface<T>(sim_ptr, matrix_manager_ptr)
{
    getReferenceToElements();
}
 
template <typename T>
void TridiagonalMatrixSolver<T>::getReferenceToElements()
{
    std::pair<int, int> mat_size{this->m_matrix_manager->lhs_matrix->size()};
    for (int i{0}; i < mat_size.first; i++)
    {
        if (i < mat_size.first - 1)
        {
            try
            {
                m_lower_sub_diag_matrix_elements.emplace_back(
                    (*(this->m_matrix_manager->lhs_matrix))(i + 1, i));
            }
            catch (const StaticSparseMatrixAccessError)
            {
                m_lower_sub_diag_matrix_elements.emplace_back(m_unused_element);
            }
 
            try
            {
                m_upper_sub_diag_matrix_elements.emplace_back(
                    (*(this->m_matrix_manager->lhs_matrix))(i, i + 1));
            }
            catch (const StaticSparseMatrixAccessError)
            {
                m_upper_sub_diag_matrix_elements.emplace_back(m_unused_element);
            }
        }
        m_diag_matrix_elements.emplace_back((*(this->m_matrix_manager->lhs_matrix))(i, i));
    }
}
 
template <typename T>
void TridiagonalMatrixSolver<T>::solveMatrixSystem(std::string_view field_name)
{
    // This is needed for CSR matrices
    if (m_lower_sub_diag_matrix_elements.empty())
    {
        getReferenceToElements();
    }
    std::string key{field_name};
    std::valarray<T>& solution{this->sim->m_field_map[key]};
    const int dim = (this->m_matrix_manager->lhs_matrix->size()).first;
    std::valarray<double> work_arr_1(dim);
    std::valarray<double> work_arr_2(dim);
    work_arr_1[0] = m_upper_sub_diag_matrix_elements[0] / m_diag_matrix_elements[0];
    work_arr_2[0] = this->m_matrix_manager->rhs_vector[0] / m_diag_matrix_elements[0];
    double denominator{0};
    for (int i = 1; i < dim - 1; i++)
    {
        denominator =
            m_diag_matrix_elements[i] - m_lower_sub_diag_matrix_elements[i - 1] * work_arr_1[i - 1];
        work_arr_1[i] = m_upper_sub_diag_matrix_elements[i] / denominator;
        work_arr_2[i] = (this->m_matrix_manager->rhs_vector[i]
                         - m_lower_sub_diag_matrix_elements[i - 1] * work_arr_2[i - 1])
                        / denominator;
    }
    denominator = m_diag_matrix_elements[dim - 1]
                  - m_lower_sub_diag_matrix_elements[dim - 2] * work_arr_1[dim - 2];
    solution[dim - 1] = (this->m_matrix_manager->rhs_vector[dim - 1]
                         - m_lower_sub_diag_matrix_elements[dim - 2] * work_arr_2[dim - 2])
                        / denominator;
    for (int i = dim - 2; i > -1; i--)
    {
        solution[i] = work_arr_2[i] - work_arr_1[i] * solution[i + 1];
    }
}
 
#endif