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{};
        std::vector<std::reference_wrapper<T>> m_rhs_vector_elements{};
 
        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)
        {
            m_lower_sub_diag_matrix_elements.emplace_back(
                (*(this->m_matrix_manager->lhs_matrix))(i + 1, i));
            m_upper_sub_diag_matrix_elements.emplace_back(
                (*(this->m_matrix_manager->lhs_matrix))(i, i + 1));
        }
        m_diag_matrix_elements.emplace_back((*(this->m_matrix_manager->lhs_matrix))(i, i));
        m_rhs_vector_elements.emplace_back(this->m_matrix_manager->rhs_vector[i]);
    }
}
 
template <typename T>
void TridiagonalMatrixSolver<T>::solveMatrixSystem(std::string_view field_name)
{
    std::string key{field_name};
    std::valarray<T>& solution{this->sim->m_field_map[key]};
    const int n = (this->m_matrix_manager->lhs_matrix->size()).first;
    std::valarray<double> work_arr_1(n);
    std::valarray<double> work_arr_2(n);
    work_arr_1[0] = m_upper_sub_diag_matrix_elements[0] / m_diag_matrix_elements[0];
    work_arr_2[0] = m_rhs_vector_elements[0] / m_diag_matrix_elements[0];
    double denominator{0};
    for (int i = 1; i < n - 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] =
            (m_rhs_vector_elements[i] - m_lower_sub_diag_matrix_elements[i - 1] * work_arr_2[i - 1])
            / denominator;
        /*if( std::isnan(work_arr_1[i]) || std::isnan(work_arr_2[i]))
        {
            std::cout << "NAN: " << i << " " << denominator << " " << work_arr_1[i] << " " <<
        work_arr_2[i] << "\n";
        }*/
    }
    denominator =
        m_diag_matrix_elements[n - 1] - m_lower_sub_diag_matrix_elements[n - 2] * work_arr_1[n - 2];
    solution[n - 1] =
        (m_rhs_vector_elements[n - 1] - m_lower_sub_diag_matrix_elements[n - 2] * work_arr_2[n - 2])
        / denominator;
    for (int i = n - 2; i > -1; i--)
    {
        solution[i] = work_arr_2[i] - work_arr_1[i] * solution[i + 1];
    }
}
 
#endif