EigenLisLinearSolver.cpp

Go to the documentation of this file.

 
#include "EigenLisLinearSolver.h"
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include "BaseLib/ConfigTree.h"
#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
#include "MathLib/LinAlg/Eigen/EigenVector.h"
#include "MathLib/LinAlg/Lis/LisLinearSolver.h"
#include "MathLib/LinAlg/Lis/LisMatrix.h"
#include "MathLib/LinAlg/Lis/LisVector.h"
 
namespace MathLib
{
EigenLisLinearSolver::EigenLisLinearSolver(
    const std::string /*solver_name*/, BaseLib::ConfigTree const* const option)
    : lis_option_(option)
{
}
 
bool EigenLisLinearSolver::solve(EigenMatrix& A_, EigenVector& b_,
                                 EigenVector& x_)
{
    static_assert(EigenMatrix::RawMatrixType::IsRowMajor,
                  "Sparse matrix is required to be in row major storage.");
    auto& A = A_.getRawMatrix();
    auto& b = b_.getRawVector();
    auto& x = x_.getRawVector();
 
    if (!A.isCompressed())
        A.makeCompressed();
    int nnz = A.nonZeros();
    int* ptr = A.outerIndexPtr();
    int* col = A.innerIndexPtr();
    double* data = A.valuePtr();
    LisMatrix lisA(A_.getNumberOfRows(), nnz, ptr, col, data);
    LisVector lisb(b.rows(), b.data());
    LisVector lisx(x.rows(), x.data());
 
    LisLinearSolver lissol;  // TODO not always create Lis solver here
    lissol.setOption(lis_option_);
    bool const status = lissol.solve(lisA, lisb, lisx);
 
    for (std::size_t i = 0; i < lisx.size(); i++)
        x[i] = lisx[i];
 
    return status;
}
 
}  // namespace MathLib