b2pastix_solver.H

//------------------------------------------------------------------------
// b2pastix_solver.H --
//
// written by Mathias Doreille
//
// (c) 2011 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// 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.
//------------------------------------------------------------------------
 
// TODO Pastix is only implemented for sequential runs!
// MPI is not supported!
#ifndef B2PASTIX_SOLVER_H_
#define B2PASTIX_SOLVER_H_
#include "b2sparse_solver.H"
 
extern "C" {
#include <pastix.h>
#include <spm.h>
}
 
namespace b2000 { namespace b2linalg {
 
template <typename T>
class PASTIX_LDLt_seq_sparse_direct_solver : public LDLt_sparse_solver<T> {
public:
    PASTIX_LDLt_seq_sparse_direct_solver() {
        UnimplementedError() << "The Pastix library is not enabled." << THROW;
    }
 
    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) {}
 
    void update_value() {}
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx,
          char left_or_right = ' ') {
        UnimplementedError() << THROW;
    }
};
 
template <typename T>
class PASTIX_LDLt_seq_extension_sparse_direct_solver : public LDLt_extension_sparse_solver<T> {
public:
    PASTIX_LDLt_seq_extension_sparse_direct_solver() {
        UnimplementedError() << "The Pastix library is not enabled." << THROW;
    }
 
    void init(
          size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
          size_t s_ext, const int connectivity, const Dictionary& dictionary) {}
 
    void 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 = ' ') {
        UnimplementedError() << THROW;
    }
};
 
// TODO This might be later needed, if we want to provide something to pastix from command line
// Get iparm and dparm
 
/*
  inline void decode_pastix_param(const Dictionary& d, int* iparam, double* dparam) {
  std::string icntl = d.get_string("PASTIX_IPARAM", "");
    if (icntl != "") {
        std::istringstream iss(icntl);
        for (;;) {
            int k, v;
            char c;
            if (iss >> k && iss >> c && iss >> v) {
                if (c != ':' || (k < 0 && k > 65))
                    Exception() << "Cannot parse the PASTIX_IPARAM analysis directive " << icntl <<
THROW; iparam[k] = v; } else Exception() << "Cannot parse the PASTIX_IPARAM analysis directive " <<
icntl << THROW; iss >> c; if (!iss) break; if (c != ',') Exception() << "Cannot parse the
PASTIX_IPARAM analysis directive " << icntl << THROW;
        }
    }
    //Get the DPARAM string from the directory d
    std::string cntl = d.get_string("PASTIX_DPARAM", "");
    if (cntl != "") {
        std::istringstream dss(cntl);
        for (;;) {
            int k;
            double v;
            char c;
            if (dss >> k && dss >> c && dss >> v) {
                if (c != ':' || (k < 0 && k > 24))
                    Exception() << "Cannot parse the PASTIX_DPARAM analysis directive " << cntl <<
THROW; dparam[k] = v; } else Exception() << "Cannot parse the PASTIX_DPARAM analysis directive " <<
cntl << THROW; dss >> c; if (!dss) break; if (c != ',') Exception() << "Cannot parse the
PASTIX_DPARAM analysis directive " << cntl << THROW;
        }
    }
}
*/
 
template <>
class PASTIX_LDLt_seq_sparse_direct_solver<double> : public LDLt_sparse_solver<double> {
public:
    // constructor
    PASTIX_LDLt_seq_sparse_direct_solver() {
        // initialize iparm and dparm
        pastixInitParam(iparm, dparm);
        // We have LDL^t => symmetric matrix type
        iparm[IPARM_FACTORIZATION] = PastixFactLDLT;
        // Number of threads to use for computation
        // TODO Adjust later for shared memory
        iparm[IPARM_THREAD_NBR] = 1;
        // Sequential run
        // TODO Adjust later for parallel execution
        iparm[IPARM_SCHEDULER] = PastixSchedSequential;
        // Factor to prevent pivoting, criteria relative to
        // vector norm, related to the convergence of the
        // equation system!
        dparm[DPARM_EPSILON_MAGN_CTRL] = -1e-14;
        // Initialize pastix
        pastixInit(&pastix_data, MPI_COMM_WORLD, iparm, dparm);
    }
 
    // destructor
    ~PASTIX_LDLt_seq_sparse_direct_solver() { pastixFinalize(&pastix_data); }
 
    void init(
          size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const double* value,
          const int connectivity, const Dictionary& dictionary) {
        if (s == 0) { return; }
        logging::Logger logger =
              logging::get_logger("linear_algebra.sparse_symmetric_solver.pastix");
 
        // Assign colind and rowind to colptr_loc and row
        colptr_loc.assign(colind, colind + s + 1);
        row.assign(rowind, rowind + nnz);
 
        // Initialize and allocate space for the sparse matrix
        spm = std::make_unique<spmatrix_t>();
        spmInit(spm.get());
 
        spm->values = const_cast<double*>(value);
        // Number of vertices in the local graph
        spm->n = s;
        // Number of non-zero entries
        spm->nnz = nnz;
        // Compressed sparse column format
        spm->fmttype = SpmCSC;
        // Matrix type, symmetric but not necessarily
        // positive-definite
        spm->mtxtype = SpmSymmetric;
        spm->colptr = colptr_loc.data();
        spm->rowptr = row.data();
        // Basval is 1 for Fortran notation and 0 for C
        spm->baseval = 0;
        // Our datatype is of double, see definition of template
        spm->flttype = SpmDouble;
 
        // Set NON computed values (e.g. gnnz or gN)
        spmUpdateComputedFields(spm.get());
 
#ifdef PASTIX_DEBUG_OUTPUT
        /* Use check and correct only if the matrix is not
         * symmetric. It is a costly call and should be avoided!
         * Useful for debugging! */
        int rc = 0;
        spmatrix_t spm2;
        // Correct spm to correspond to PaStiX standards
        rc = spmCheckAndCorrect(spm.get(), &spm2);
        if (rc != 0) {
            spmExit(spm.get());
            *spm = spm2;
            rc = 0;
        }
 
        // Dump the matrix to a file, used for debugging!
        FILE* myMatrix = fopen("spm.out", "a");
        spmSave(spm.get(), myMatrix);
        fclose(myMatrix);
 
#endif  // PASTIX_DEBUG_OUTPUT
 
        /*Perform all the preprocessing steps:
         * ordering, symbolic factorization,
         * reordering, proportionnal mapping */
        pastix_task_analyze(pastix_data, spm.get());
 
        /*Perform all the numerical factorization
         * steps: fill the internal block CSC and
         * the solver matrix structures, then apply
         * the factorization step. */
        pastix_task_numfact(pastix_data, spm.get());
 
        updated_value = false;
    }
 
    void update_value() { updated_value = true; }
 
    void resolve(
          size_t s, size_t nrhs, const double* b, size_t ldb, double* x, size_t ldx,
          char left_or_right = ' ') {
        if (s == 0) { return; }
        logging::Logger logger =
              logging::get_logger("linear_algebra.sparse_symmetric_solver.pastix");
 
        if (updated_value) {
            pastix_task_numfact(pastix_data, spm.get());
            updated_value = false;
        }
 
#ifdef PASTIX_DEBUG_OUTPUT
        // Dump the rhs to a file, meant for debugging!
 
        FILE* myMatrix = fopen("rhs.out", "a");
        spmPrintRHS(spm.get(), nrhs, x, ldx, myMatrix);
        fclose(myMatrix);
 
#endif  // PASTIX_DEBUG_OUTPUT
 
        if (nrhs == 1) {
            if (b != x) { std::copy(b, b + s, x); }
            pastix_task_solve(pastix_data, nrhs, x, ldx);
 
        } else {
            std::vector<double> x_value(ldx * nrhs);
 
            for (size_t i = 0; i != nrhs; ++i) {
                std::copy(b + i * ldb, b + i * ldb + s, x_value.begin() + i * ldx);
            }
            pastix_task_solve(pastix_data, nrhs, x_value.data(), ldx);
 
            // Copy the data back to x
            for (size_t i = 0; i != nrhs; ++i) {
                std::copy(x_value.begin() + i * ldx, x_value.begin() + i * ldx + s, x + i * ldx);
            }
        }
    }
 
protected:
    // Pointer to a storage structure
    pastix_data_t* pastix_data{nullptr};
    bool updated_value{false};
    // Used to copy rowind and colind
    std::vector<int> colptr_loc;
    std::vector<int> row;
    // Sparse matrix storage
    std::unique_ptr<spmatrix_t> spm;
    spm_int_t iparm[IPARM_SIZE];
    double dparm[DPARM_SIZE];
};
 
template <>
class PASTIX_LDLt_seq_extension_sparse_direct_solver<double>
    : public LDLt_extension_sparse_solver<double>,
      public PASTIX_LDLt_seq_sparse_direct_solver<double> {
public:
    PASTIX_LDLt_seq_extension_sparse_direct_solver()
        : PASTIX_LDLt_seq_sparse_direct_solver<double>(), div(0) {}
 
    void init(
          size_t size_, size_t nnz_, const size_t* colind_, const size_t* rowind_,
          const double* value_, size_t size_ext_, const int connectivity,
          const Dictionary& dictionary) {
        if (size_ext_ != 1) { UnimplementedError() << THROW; }
 
        PASTIX_LDLt_seq_sparse_direct_solver<double>::init(
              size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
        m_ab.resize(size_ * 2);
    }
 
    void update_value() { PASTIX_LDLt_seq_sparse_direct_solver<double>::update_value(); }
 
    void resolve(
          size_t s, size_t nrhs, const double* b, size_t ldb, double* x, size_t ldx,
          const double* ma_ = 0, const double* mb_ = 0, const double* mc_ = 0,
          char left_or_right = ' ') {
        if (mc_ != 0) {
            std::copy(ma_, ma_ + s - 1, &m_ab[0]);
            std::copy(mb_, mb_ + s - 1, &m_ab[s - 1]);
            PASTIX_LDLt_seq_sparse_direct_solver<double>::resolve(
                  s - 1, 2, &m_ab[0], s - 1, &m_ab[0], s - 1);
            div = 1 / (*mc_ - blas::dot(s - 1, ma_, 1, &m_ab[s - 1], 1));
        }
 
        for (size_t i = 0; i != nrhs; ++i) {
            const double x2 = x[ldx * i + s - 1] =
                  div * (b[ldb * i + s - 1] - blas::dot(s - 1, b + ldb * i, 1, &m_ab[s - 1], 1));
            PASTIX_LDLt_seq_sparse_direct_solver<double>::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);
        }
    }
 
protected:
    std::vector<double> m_ab;
    double div;
};
 
}}  // namespace b2000::b2linalg
 
#endif /* B2PASTIX_SOLVER_H_ */