b2mumps_solver.H

Go to the documentation of this file.

//------------------------------------------------------------------------
// b2mumps_solver.H --
//
// written by Mathias Doreille
//            Harald Klimach <harald.klimach@dlr.de>
//
// (c) 2011-2016 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (c) 2021-2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//          Linder Höhe, 51147 Köln
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR or DLR.
//
//------------------------------------------------------------------------
 
#ifndef B2MUMPS_SOLVER_H_
#define B2MUMPS_SOLVER_H_
 
#include <type_traits>
 
#include "b2sparse_solver.H"
#include "utils/b2timing.H"
 
#ifdef HAVE_UPMUMPS_DMUMPS_C_H
#include "MUMPS/dmumps_c.h"
#else
#ifdef HAVE_MUMPS_DMUMPS_C_H
#include "mumps/dmumps_c.h"
#else
#ifdef HAVE_DMUMPS_C_H
#include "dmumps_c.h"
#endif  // HAVE_DMUMPS_C_H
#endif  // HAVE_MUMPS_DMUMPS_C_H
#endif  // HAVE_UPMUMPS_DMUMPS_C_H
 
#ifdef HAVE_UPMUMPS_ZMUMPS_C_H
#include "MUMPS/zmumps_c.h"
#else
#ifdef HAVE_MUMPS_ZMUMPS_C_H
#include "mumps/zmumps_c.h"
#else
#ifdef HAVE_ZMUMPS_C_H
#include "zmumps_c.h"
#endif  // HAVE_ZMUMPS_C_H
#endif  // HAVE_MUMPS_ZMUMPS_C_H
#endif  // HAVE_UPMUMPS_ZMUMPS_C_H
 
#include "b2blaslapack.H"
#include "b2linalg_solver_slave.H"
#include "utils/b2logging.H"
 
//------------------------------------------------------------------------
//
//------------------------------------------------------------------------
 
namespace b2000 { namespace b2linalg {
 
 
enum MUMPSICNTL : int {
    MUMPS_ErrStream = 1 - 1,
    MUMPS_DiagStream = 2 - 1,
    MUMPS_GlobStream = 3 - 1,
    MUMPS_LogLevel = 4 - 1,
    MUMPS_MatrixForm = 5 - 1,
    MUMPS_ZeroFreeScale = 6 - 1,
    MUMPS_SymPerm = 7 - 1,
    MUMPS_Scaling = 8 - 1,
    MUMPS_Transposed = 9 - 1,
    MUMPS_MaxIterRefine = 10 - 1,
    MUMPS_ErrStats = 11 - 1,
    MUMPS_OrdStrategy = 12 - 1,
    MUMPS_RootParallel = 13 - 1,
    MUMPS_WsIncrease = 14 - 1,
    MUMPS_nThreads = 16 - 1,
    MUMPS_MatrixDist = 18 - 1,
    MUMPS_SchurComp = 19 - 1,
    MUMPS_RhsForm = 20 - 1,
    MUMPS_SolDist = 21 - 1,
    MUMPS_OOC = 22 - 1,
    MUMPS_MaxWorkMem = 23 - 1,
    MUMPS_NullPivotRows = 24 - 1,
    MUMPS_Deficient = 25 - 1,
    MUMPS_SchurPhase = 26 - 1,
    MUMPS_RhsBlockSize = 27 - 1,
    MUMPS_SeqOrPar = 28 - 1,
    MUMPS_ParOrderTool = 29 - 1,
    MUMPS_Subset = 30 - 1,
    MUMPS_DiscardFacts = 31 - 1,
    MUMPS_ForwardInFact = 32 - 1,
    MUMPS_Determinant = 33 - 1,
    MUMPS_DelFiles = 34 - 1,
    MUMPS_BLR = 35 - 1,
    MUMPS_BlrVariant = 36 - 1,
    MUMPS_Compression = 38 - 1,
};
 
 
enum class MUMPSACTION : unsigned int {
    NewObject = 1,
    Continued = 2,
};
 
template <class T>
struct B2MUMPS_STRUC;
 
template <>
struct B2MUMPS_STRUC<double> {
    using type = DMUMPS_STRUC_C;
    using numtype = DMUMPS_REAL;
};
template <>
struct B2MUMPS_STRUC<std::complex<double>> {
    using type = ZMUMPS_STRUC_C;
    using numtype = mumps_double_complex;
};
template <>
struct B2MUMPS_STRUC<b2000::csda<double>> {
    using type = ZMUMPS_STRUC_C;
    using numtype = mumps_double_complex;
};
 
template <class T>
using B2MUMPS_STRUC_t = typename B2MUMPS_STRUC<T>::type;
 
template <class T>
using B2MUMPS_NUM_t = typename B2MUMPS_STRUC<T>::numtype;
 
 
inline void mumps_error_handler(const int* const info, bool autospc = false) {
    const int error = info[0];
    if (error == 0) {
        if (!autospc && info[27] > 0) {
            SingularMatrixError e;
            e.null_space_size = info[27];
            e << "Numerically singular matrix: " << info[27] << " null pivot found\n";
            e << "Note: You may want to reduce the threshold for null pivots\n";
            e << "      by adding a mumps_cntl(3) statement to your case definition.\n";
            e << "      Try `mumps_cntl '3:1.1e-21'`.\n";
            e << "See the section on the mumps_cntl parameter in the user guide for\n";
            e << "further reference.\n";
            e << THROW;
        }
    } else if (error > 0) {
        logging::Logger logger =
              logging ::get_logger("linear_algebra.sparse_symmetric_solver.mumps");
        logger << logging::warning << "warning code " << error << ": ";
        if (error & 1) { logger << "Index (in IRN or JCN) out of range. "; }
        if (error & 2) {
            logger << "During error analysis the max-norm of the computed "
                      "solution was found to be zero. ";
        }
        if (error & 4) {
            logger << "User data JCN has been modified (internally) by the"
                      " solver. ";
        }
        if (error & 8) { logger << "More than ICNTL(10) iterations refinement are required."; }
        logger << logging::LOGFLUSH;
    } else {
        Exception e;
        e << "error code " << error << ": ";
        switch (error) {
            case -6:
                e << "Matrix is singular in structure.";
                break;
            case -10:
                e << "Numerically singular matrix.";
                break;
            case -13:
                e << "out of memory error";
                break;
            case -40:
                e << "The matrix was indicated to be positive definite but a"
                     " negative or null pivot was encountered";
                break;
            default:
                break;
        }
        e << THROW;
    }
}
 
 
template <typename MUMPS_STRUC_C>
inline void decode_mumps_cntl(const Dictionary& d, MUMPS_STRUC_C& mumps_struct) {
    std::string icntl = d.get_string("MUMPS_ICNTL", "");
    if (icntl != "") {
        std::istringstream iss(icntl);
        for (;;) {
            int k, v;
            char c;
            if (iss >> k && iss >> c && iss >> v) {
                if (c != ':' || (k < 1 || k > 33)) {
                    Exception() << "Cannot parse the MUMPS_ICNTL analysis directive " << icntl
                                << THROW;
                }
                if (MUMPS_MatrixForm == k - 1) {
                    Exception() << "Not allowed to modify the Matrix form "
                                   "(5) in MUMPS_ICNTL: "
                                << icntl << THROW;
                }
                if (MUMPS_MatrixDist == k - 1) {
                    Exception() << "Not allowed to modify the Matrix "
                                   "distribution (18) in MUMPS_ICNTL: "
                                << icntl << THROW;
                }
                mumps_struct.icntl[k - 1] = v;
            } else {
                Exception() << "Cannot parse the MUMPS_ICNTL analysis directive " << icntl << THROW;
            }
            iss >> c;
            if (!iss) { break; }
            if (c != ',') {
                Exception() << "Cannot parse the MUMPS_ICNTL analysis directive " << icntl << THROW;
            }
        }
    }
    std::string cntl = d.get_string("MUMPS_CNTL", "");
    if (cntl != "") {
        std::istringstream dss(cntl);
        for (;;) {
            int k;
            double v;
            char c;
            if (dss >> k && dss >> c && dss >> v) {
                if (c != ':' || (k < 1 || k > 5)) {
                    Exception() << "Cannot parse the MUMPS_CNTL analysis directive " << cntl
                                << THROW;
                }
                mumps_struct.cntl[k - 1] = v;
            } else {
                Exception() << "Cannot parse the MUMPS_CNTL analysis directive " << cntl << THROW;
            }
            dss >> c;
            if (!dss) { break; }
            if (c != ',') {
                Exception() << "Cannot parse the MUMPS_CNTL analysis directive " << cntl << THROW;
            }
        }
    }
}
 
 
template <typename MUMPS_STRUC_C>
inline void b2_mumps_call(
      MUMPS_STRUC_C& mumps_struct, MUMPSACTION action = MUMPSACTION::Continued) {
    constexpr bool is_dmumps = std::is_same<MUMPS_STRUC_C, DMUMPS_STRUC_C>::value;
    constexpr unsigned int numtype = is_dmumps ? 0 : 1;
 
    unsigned int act = static_cast<unsigned int>(action);
    unsigned int sym_job = 1 == act ? mumps_struct.sym : mumps_struct.job;
    auto mumpstimer = obtain_timer("mumps");
    (*mumpstimer).second.start();
    SolverSlave::MumpsSend(act, numtype, &mumps_struct, sym_job);
    if constexpr (is_dmumps) {
        dmumps_c(&mumps_struct);
    } else {
        zmumps_c(&mumps_struct);
    }
    (*mumpstimer).second.stop();
}
 
 
template <typename MUMPS_STRUC_C>
inline void b2_mumps_new(MUMPS_STRUC_C& mumps_struct, MUMPS_INT sym) {
    mumps_struct.sym = sym;
    mumps_struct.par = 1;
    mumps_struct.job = -1;
    mumps_struct.comm_fortran = (MUMPS_INT)SolverSlave::GetFcomm();
    b2_mumps_call(mumps_struct, MUMPSACTION::NewObject);
    mumps_error_handler(mumps_struct.infog);
    mumps_struct.irn = 0;
    mumps_struct.jcn = 0;
 
    logging::Logger logger = logging::get_logger("linear_algebra.sparse_symmetric_solver.mumps");
    if (logger.is_enabled_for(logging::debug)) {
        mumps_struct.icntl[MUMPS_LogLevel] = 3;
        mumps_struct.icntl[MUMPS_ErrStream] = 6;
        mumps_struct.icntl[MUMPS_DiagStream] = 6;
        mumps_struct.icntl[MUMPS_GlobStream] = 6;
        mumps_struct.icntl[MUMPS_ErrStats] = 1;
    } else {
        mumps_struct.icntl[MUMPS_LogLevel] = 0;
        mumps_struct.icntl[MUMPS_ErrStream] = 0;
        mumps_struct.icntl[MUMPS_DiagStream] = 0;
        mumps_struct.icntl[MUMPS_GlobStream] = 0;
    }
}
 
 
template <typename MUMPS_STRUC_C>
inline void b2_mumps_delete(MUMPS_STRUC_C& mumps_struct) {
    mumps_struct.job = -2;
    b2_mumps_call(mumps_struct);
    mumps_error_handler(mumps_struct.infog, true);
    delete[] mumps_struct.irn;
    delete[] mumps_struct.jcn;
}
 
 
template <typename MUMPS_STRUC_C>
inline void b2_mumps_call_wsgrow(MUMPS_STRUC_C& mumps_struct) {
    logging::Logger logger = logging::get_logger("linear_algebra.sparse_symmetric_solver.mumps");
    for (;;) {
        b2_mumps_call(mumps_struct);
        if (mumps_struct.infog[0] == -8 || mumps_struct.infog[0] == -9) {
            mumps_struct.icntl[MUMPS_WsIncrease] *= 2;
            if (4 == mumps_struct.job) { mumps_struct.job = 2; }
            logger << "resolve: work array is too small, increasing size "
                      "of work array and retrying resolution."
                   << logging::LOGFLUSH;
        } else {
            break;
        }
    }
}
 
 
template <typename T>
inline void b2_mumps_init_struct(
      B2MUMPS_STRUC_t<T>& mumps_struct, bool& autospc, bool& updated, size_t s, size_t nnz,
      const size_t* const colind, const size_t* const rowind, const T* const value,
      const Dictionary& dictionary) {
    autospc = dictionary.get_bool("AUTOSPC", false);
    mumps_struct.icntl[MUMPS_NullPivotRows] = 1;
    mumps_struct.icntl[MUMPS_OOC] = dictionary.get_bool("MUMPS_OOC", 0);
    mumps_struct.icntl[MUMPS_Scaling] = dictionary.get_int("MUMPS_SCALING", 77);
    mumps_struct.icntl[MUMPS_MaxIterRefine] =
          dictionary.get_int("MUMPS_MAX_ITERATIVE_REFINEMENT", 0);
    mumps_struct.cntl[2] = dictionary.get_double("AUTOSPC_THRESHOLD", 0.0);
    mumps_struct.cntl[4] = dictionary.get_double("AUTOSPC_VALUE", 0.0);
    decode_mumps_cntl(dictionary, mumps_struct);
 
    if (mumps_struct.icntl[MUMPS_OOC] != 0) {
        std::string tmp = dictionary.get_string("MUMPS_OOC_TMPDIR", "");
        if (!tmp.empty()) {
            if (tmp.size() > 255) {
                Exception() << "The MUMPS_OOC_TMPDIR string is limited to 255"
                               " characters."
                            << THROW;
            }
            strncpy(mumps_struct.ooc_tmpdir, tmp.c_str(), 256);
        }
        tmp = dictionary.get_string("MUMPS_OOC_PREFIX", "");
        if (!tmp.empty()) {
            if (tmp.size() > 63) {
                Exception() << "The MUMPS_OOC_PREFIX string is limited to 63"
                               " characters."
                            << THROW;
            }
            strncpy(mumps_struct.ooc_prefix, tmp.c_str(), 64);
        }
    }
 
    mumps_struct.job = 4;
    mumps_struct.n = s;
    mumps_struct.nz = nnz;
    mumps_struct.irn = new int[nnz];
    mumps_struct.jcn = new int[nnz];
    if (s == 0) { return; }
    size_t ptr = colind[0];
    for (size_t j = 1; j <= s; ++j) {
        const size_t ptr_end = colind[j];
        for (; ptr != ptr_end; ++ptr) {
            mumps_struct.irn[ptr] = rowind[ptr] + 1;
            mumps_struct.jcn[ptr] = j;
        }
    }
    mumps_struct.a =
          const_cast<B2MUMPS_NUM_t<T>*>(reinterpret_cast<const B2MUMPS_NUM_t<T>*>(value));
    b2_mumps_call_wsgrow(mumps_struct);
    if ((mumps_struct.sym > 0) && (mumps_struct.infog[0] == -6)) {
        std::vector<bool> dof_nn(s, false);
        size_t r = colind[0];
        for (size_t j = 0; j != s; ++j) {
            for (size_t r_end = colind[j + 1]; r != r_end; ++r) {
                dof_nn[j] = dof_nn[rowind[r]] = (value[r] != T{0.0});
            }
        }
        SingularMatrixError e;
        for (size_t j = 0; j != s; ++j) {
            if (!dof_nn[j]) { e.singular_dofs.push_back(j); }
        }
        e << "Matrix is singular in structure. The involved " << e.singular_dofs.size()
          << " dof's can be visualized by means of the DISP_NS dataset." << THROW;
    }
    mumps_error_handler(mumps_struct.infog, autospc);
    updated = false;
}
 
 
template <typename T>
void mumps_solve(
      B2MUMPS_STRUC_t<T>& mumps_, bool& updated_value, bool autospc, size_t s, size_t nrhs,
      const T* const b, size_t ldb, T* const x, size_t ldx) {
    if (s == 0) { return; }
    mumps_.job = updated_value ? 5 : 3;
    updated_value = false;
 
    if (nrhs == 1 || mumps_.icntl[MUMPS_ErrStats] == 0) {
        mumps_.rhs = reinterpret_cast<B2MUMPS_NUM_t<T>*>(x);
        mumps_.nrhs = nrhs;
        mumps_.lrhs = ldx;
        if (b != x) {
            for (size_t i = 0; i != nrhs; ++i) {
                std::copy(b + i * ldb, b + i * ldb + s, x + i * ldx);
            }
        }
        b2_mumps_call_wsgrow(mumps_);
        mumps_error_handler(mumps_.infog, autospc);
    } else {
        // solve the multiple rhs separately to get the error analyse on
        // each system
        mumps_.nrhs = 1;
        mumps_.lrhs = s;
        for (size_t i = 0; i != nrhs; ++i) {
            mumps_.rhs = reinterpret_cast<B2MUMPS_NUM_t<T>*>(x + i * ldx);
            if (b != x) { std::copy(b + i * ldb, b + i * ldb + s, x + i * ldx); }
            b2_mumps_call_wsgrow(mumps_);
            mumps_error_handler(mumps_.infog, autospc);
        }
    }
}
 
//
// LDLt
//
 
template <typename T>
class MUMPS_LDLt_seq_sparse_direct_solver : public LDLt_sparse_solver<T> {
public:
 
    explicit MUMPS_LDLt_seq_sparse_direct_solver(bool definit_positive = false) {
        b2_mumps_new(mumps_, definit_positive ? 1 : 2);
    }
 
    ~MUMPS_LDLt_seq_sparse_direct_solver() override { b2_mumps_delete(mumps_); }
 
    void init(
          size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
          const int connectivity, const Dictionary& dictionary) override {
        b2_mumps_init_struct(
              mumps_, autospc, updated_value, s, nnz, colind, rowind, value, dictionary);
    }
 
    void update_value() override { updated_value = true; }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx,
          char left_or_right = ' ') override {
        mumps_solve(mumps_, updated_value, autospc, s, nrhs, b, ldb, x, ldx);
    }
 
    size_t get_null_space_size() override { return mumps_.infog[27]; }
 
    void get_null_space(size_t s, size_t nv, T* m, size_t lm) override {
        if (s == 0) { return; }
 
        mumps_.job = updated_value ? 5 : 3;
        updated_value = false;
        mumps_.icntl[MUMPS_NullPivotRows] = 1;
        mumps_.icntl[MUMPS_Deficient] = nv == 1 ? 1 : -1;
        mumps_.rhs = reinterpret_cast<B2MUMPS_NUM_t<T>*>(m);
        mumps_.nrhs = nv;
        mumps_.lrhs = lm;
        b2_mumps_call_wsgrow(mumps_);
        mumps_.icntl[MUMPS_Deficient] = 0;
        mumps_error_handler(mumps_.infog, true);
    }
 
private:
    bool updated_value = false;
    bool autospc = false;
    B2MUMPS_STRUC_t<T> mumps_;
};  // class MUMPS_LDLt_seq_sparse_direct_solver
 
template <typename T>
class MUMPS_LDLt_seq_extension_sparse_direct_solver
    : public LDLt_extension_sparse_solver<T>,
      public MUMPS_LDLt_seq_sparse_direct_solver<T> {
public:
 
    explicit MUMPS_LDLt_seq_extension_sparse_direct_solver(bool definit_positive = false)
        : MUMPS_LDLt_seq_sparse_direct_solver<T>(definit_positive) {}
 
    void init(
          size_t size_, size_t nnz_, const size_t* colind_, const size_t* rowind_, const T* value_,
          size_t size_ext_, const int connectivity, const Dictionary& dictionary) override {
        if (size_ext_ != 1) { UnimplementedError() << THROW; }
 
        MUMPS_LDLt_seq_sparse_direct_solver<T>::init(
              size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
        m_ab.resize(size_ * 2);
    }
 
    void update_value() override { MUMPS_LDLt_seq_sparse_direct_solver<T>::update_value(); }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx, const T* ma_ = 0,
          const T* mb_ = 0, const T* mc_ = 0, char left_or_right = ' ') override {
        if (mc_ != 0) {
            std::copy(ma_, ma_ + s - 1, &m_ab[0]);
            std::copy(mb_, mb_ + s - 1, &m_ab[s - 1]);
            MUMPS_LDLt_seq_sparse_direct_solver<T>::resolve(
                  s - 1, 2, &m_ab[0], s - 1, &m_ab[0], s - 1);
            m_div = T{1.0} / (*mc_ - blas::dot(s - 1, ma_, 1, &m_ab[s - 1], 1));
        }
 
        for (size_t i = 0; i != nrhs; ++i) {
            const T x2 = x[ldx * i + s - 1] =
                  m_div * (b[ldb * i + s - 1] - blas::dot(s - 1, b + ldb * i, 1, &m_ab[s - 1], 1));
            MUMPS_LDLt_seq_sparse_direct_solver<T>::resolve(
                  s - 1, 1, b + ldb * i, ldb, x + ldx * i, ldx);
            blas::axpy(s - 1, -x2, &m_ab[0], 1, x + ldx * i, 1);
        }
    }
 
    size_t get_null_space_size() override {
        return MUMPS_LDLt_seq_sparse_direct_solver<T>::get_null_space_size();
    }
 
    void get_null_space(size_t s, size_t nv, T* m, size_t lm) override {
        MUMPS_LDLt_seq_sparse_direct_solver<T>::get_null_space(s, nv, m, lm);
    }
 
protected:
    std::vector<T> m_ab;
    T m_div{0.0};
};  // class MUMPS_LDLt_seq_extension_sparse_direct_solver
 
//
// LU
//
 
template <typename T>
class MUMPS_LU_seq_sparse_direct_solver : public LU_sparse_solver<T> {
public:
    MUMPS_LU_seq_sparse_direct_solver() { b2_mumps_new(mumps_, 0); }
 
    ~MUMPS_LU_seq_sparse_direct_solver() override { b2_mumps_delete(mumps_); }
 
    void init(
          size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
          const int connectivity, const Dictionary& dictionary) override {
        b2_mumps_init_struct(
              mumps_, autospc, updated_value, s, nnz, colind, rowind, value, dictionary);
    }
 
    void update_value() override { updated_value = true; }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx,
          char left_or_right = ' ') override {
        mumps_solve(mumps_, updated_value, autospc, s, nrhs, b, ldb, x, ldx);
    }
 
protected:
    bool updated_value = false;
    bool autospc = false;
    B2MUMPS_STRUC_t<T> mumps_;
};  // class MUMPS_LU_seq_sparse_direct_solver
 
template <typename T>
class MUMPS_LU_seq_extension_sparse_direct_solver : public LU_extension_sparse_solver<T>,
                                                    public MUMPS_LU_seq_sparse_direct_solver<T> {
public:
    void init(
          size_t size_, size_t nnz_, const size_t* colind_, const size_t* rowind_, const T* value_,
          size_t size_ext_, const int connectivity, const Dictionary& dictionary) override {
        if (size_ext_ != 1) { UnimplementedError() << THROW; }
 
        MUMPS_LU_seq_sparse_direct_solver<T>::init(
              size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
        m_ab.resize(size_ * 2);
    }
 
    void update_value() override { MUMPS_LU_seq_sparse_direct_solver<T>::update_value(); }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx, const T* ma_ = 0,
          const T* mb_ = 0, const T* mc_ = 0, char left_or_right = ' ') override {
        if (mc_ != 0) {
            std::copy(ma_, ma_ + s - 1, &m_ab[0]);
            std::copy(mb_, mb_ + s - 1, &m_ab[s - 1]);
            MUMPS_LU_seq_sparse_direct_solver<T>::resolve(
                  s - 1, 2, &m_ab[0], s - 1, &m_ab[0], s - 1);
            m_div = T{1.0} / (*mc_ - blas::dot(s - 1, ma_, 1, &m_ab[s - 1], 1));
        }
 
        for (size_t i = 0; i != nrhs; ++i) {
            const T x2 = x[ldx * i + s - 1] =
                  m_div * (b[ldb * i + s - 1] - blas::dot(s - 1, b + ldb * i, 1, &m_ab[s - 1], 1));
            MUMPS_LU_seq_sparse_direct_solver<T>::resolve(
                  s - 1, 1, b + ldb * i, s - 1, x + ldx * i, s - 1);
            blas::axpy(s - 1, -x2, &m_ab[0], 1, x + ldx * i, 1);
        }
    }
 
protected:
    std::vector<T> m_ab;
    T m_div{};
};  // class MUMPS_LU_seq_extension_sparse_direct_solver
 
}}  // namespace b2000::b2linalg
 
#endif  // B2MUMPS_SOLVER_H_