b2gmm_solver.H

//------------------------------------------------------------------------
// b2gmm_solver.H --
//
//
// written by Mathias Doreille
//
// Copyright (c) 2004-2012 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.
//------------------------------------------------------------------------
 
#ifndef __B2GMM_SOLVER_H__
#define __B2GMM_SOLVER_H__
 
// #pragma GCC diagnostic push
// #pragma GCC diagnostic warning "-fpermissive"
 
#include "gmm/gmm_dense_Householder.h"
#include "gmm/gmm_interface.h"
#include "gmm/gmm_precond_diagonal.h"
#include "gmm/gmm_precond_ildlt.h"
#include "gmm/gmm_precond_ildltt.h"
#include "gmm/gmm_precond_ilu.h"
#include "gmm/gmm_precond_ilut.h"
#include "gmm/gmm_precond_ilutp.h"
#include "gmm/gmm_precond_mr_approx_inverse.h"
#include "gmm/gmm_solver_cg.h"
#include "gmm/gmm_solver_gmres.h"
 
// #pragma GCC diagnostic pop
 
#include "b2ppconfig.h"
#include "utils/b2exception.H"
 
namespace b2000 { namespace b2linalg {
 
template <typename T, template <typename MATRIX> class PRECOND>
class gmm_cg_iterative_solver : public LDLt_sparse_solver<T> {
public:
    typedef gmm::csc_matrix<T> gmm_matrix_t;
    typedef PRECOND<gmm_matrix_t> precond_type;
 
    gmm_cg_iterative_solver(precond_type* precond_, const Dictionary& dictionary_)
        : colind(0),
          rowind(0),
          value(0),
          precond(precond_),
          modified_value(true),
          dictionary(dictionary_) {}
 
    ~gmm_cg_iterative_solver() { delete precond; }
 
    void init(
          size_t s_, size_t nnz_, const size_t* colind_, const size_t* rowind_, const T* value_,
          int connectivity, const Dictionary& dictionary_) {
        gmm_matrix.nc = s_;
        gmm_matrix.nr = s_;
        colind = colind_;
        rowind = rowind_;
        value = value_;
        size_t nnz = colind[gmm_matrix.nc];
        nnz = nnz * 2 - gmm_matrix.nc;
#if GMM_VERSION == 4
        gmm_matrix.pr.resize(nnz);
        gmm_matrix.ir.resize(nnz);
        gmm_matrix.jc.resize(gmm_matrix.nc + 2);
#else
        gmm_matrix.pr = new T[nnz];
        gmm_matrix.ir = new typename gmm_matrix_t::IND_TYPE[nnz];
        gmm_matrix.jc = new typename gmm_matrix_t::IND_TYPE[gmm_matrix.nc + 2];
#endif
        update_value();
    }
 
    void update_value() { modified_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 = ' ') {
        if (modified_value) {
            size_t nnz = colind[gmm_matrix.nc];
            nnz = nnz * 2 - gmm_matrix.nc;
            size_t size = gmm_matrix.nc;
#if GMM_VERSION == 4
            typename gmm_matrix_t::IND_TYPE* colptr1 = &gmm_matrix.jc[0];
#else
            typename gmm_matrix_t::IND_TYPE* colptr1 = gmm_matrix.jc;
#endif
            colptr1[0] = colptr1[1] = 0;
            colptr1 += 2;
            for (size_t j = 0; j != size; ++j) { colptr1[j] = colind[j + 1] - colind[j]; }
            for (size_t j = 0; j != size; ++j) {
                const size_t i_end = colind[j + 1];
                for (size_t i = colind[j] + 1; i < i_end; ++i) { ++colptr1[rowind[i]]; }
            }
            --colptr1;
            for (size_t j = 0; j != size; ++j) { colptr1[j + 1] += colptr1[j]; }
#if GMM_VERSION == 4
            typename gmm_matrix_t::IND_TYPE* rowind1 = &gmm_matrix.ir[0];
            T* nzval1 = &gmm_matrix.pr[0];
#else
            typename gmm_matrix_t::IND_TYPE* rowind1 = gmm_matrix.ir;
            T* nzval1 = gmm_matrix.pr;
#endif
            for (size_t j = 0; j != size; ++j) {
                size_t i = colind[j];
                const size_t i_end = colind[j + 1];
                typename gmm_matrix_t::IND_TYPE ii = colptr1[j];
                colptr1[j] += i_end - i;
                std::copy(rowind + i, rowind + i_end, rowind1 + ii);
                std::copy(value + i, value + i_end, nzval1 + ii);
                for (++i; i < i_end; ++i) {
                    int ii = colptr1[rowind[i]]++;
                    rowind1[ii] = j;
                    nzval1[ii] = value[i];
                }
            }
            precond->build_with(gmm_matrix);
            modified_value = false;
        }
 
        if (left_or_right != ' ') { UnimplementedError() << THROW; }
 
        gmm::iteration iter(
              dictionary.get_double("SLS_TOL_RESIDUUM", 1.0E-8),
              dictionary.get_int("SLS_VERBOSITY", 0), dictionary.get_int("SLS_MAX_ITER", -1));
        T* tmp = 0;
        if (b == x) { tmp = new T[s]; }
        for (size_t i = 0; i != nrhs; ++i) {
            if (b == x) { std::copy(b + ldb * i, b + ldb * (i + 1), tmp); }
            gmm::array1D_reference<const T*> B(b == x ? tmp : b + ldb * i, s);
            gmm::array1D_reference<T*> X(x + ldx * i, s);
            gmm::cg(gmm_matrix, X, B, *precond, iter);
            if (!iter.converged()) {
                Exception() << "The iterative sparse linear solver did not converge" << THROW;
            }
        }
        if (b == x) { delete[] tmp; }
    }
 
protected:
    const size_t* colind;
    const size_t* rowind;
    const T* value;
    gmm_matrix_t gmm_matrix;
    precond_type* precond;
    bool modified_value;
    const Dictionary& dictionary;
};
 
template <typename T>
LDLt_sparse_solver<T>* new_gmm_cg_iterative_solver(const Dictionary& dictionary) {
    std::string precond_name = dictionary.get_string("SLS_PRECOND", "ILDLT");
    if (precond_name == "ILDLT") {
        typedef gmm_cg_iterative_solver<T, gmm::ildlt_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "ILDLTT") {
        typedef gmm_cg_iterative_solver<T, gmm::ildltt_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "DIAGONAL") {
        typedef gmm_cg_iterative_solver<T, gmm::diagonal_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    // if (precond_name == "_mr_approx_inverse") {
    //    typedef gmm_cg_iterative_solver<T, gmm::mr_approx_inverse_precond> s_t;
    //    return new s_t(new typename s_t::precond_type(), dictionary);
    //}
    Exception() << "Unknown GMM++ cg preconditioner " << precond_name << "." << THROW;
    return 0;
}
 
template <typename T, template <typename MATRIX> class PRECOND>
class extension_gmm_cg_iterative_solver : public LDLt_extension_sparse_solver<T>,
                                          public gmm_cg_iterative_solver<T, PRECOND> {
public:
    extension_gmm_cg_iterative_solver(
          typename gmm_cg_iterative_solver<T, PRECOND>::precond_type* precond_,
          const Dictionary& dictionary_)
        : gmm_cg_iterative_solver<T, PRECOND>(precond_, dictionary_) {}
 
    ~extension_gmm_cg_iterative_solver() {}
 
    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) {
        if (size_ext_ != 1) { UnimplementedError() << THROW; }
 
        gmm_cg_iterative_solver<T, PRECOND>::init(
              size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
        m_a.resize(size_);
        m_b.resize(size_);
    }
 
    void update_value() { gmm_cg_iterative_solver<T, PRECOND>::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 = ' ') {
        if (ma_ != 0 && mb_ != 0 && mc_ != 0) {
            gmm_cg_iterative_solver<T, PRECOND>::resolve(s - 1, 1, ma_, s - 1, &m_a[0], s - 1);
            gmm_cg_iterative_solver<T, PRECOND>::resolve(s - 1, 1, mb_, s - 1, &m_b[0], s - 1);
            div = 1 / (*mc_ - blas::dot(s - 1, ma_, 1, &m_b[0], 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_b[0], 1));
            gmm_cg_iterative_solver<T, PRECOND>::resolve(
                  s - 1, 1, b + ldb * i, s - 1, x + ldx * i, s - 1);
            blas::axpy(s - 1, -x2, &m_a[0], 1, x + ldx * i, 1);
        }
    }
 
protected:
    std::vector<T> m_a;
    std::vector<T> m_b;
    T div;
};
 
template <typename T>
LDLt_extension_sparse_solver<T>* new_extension_gmm_cg_iterative_solver(
      const Dictionary& dictionary) {
    std::string precond_name = dictionary.get_string("SLS_PRECOND", "ILDLT");
    if (precond_name == "ILDLT") {
        typedef extension_gmm_cg_iterative_solver<T, gmm::ildlt_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "ILDLTT") {
        typedef extension_gmm_cg_iterative_solver<T, gmm::ildltt_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "DIAGONAL") {
        typedef extension_gmm_cg_iterative_solver<T, gmm::diagonal_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    // if (precond_name == "_mr_approx_inverse") {
    //    typedef gmm_cg_iterative_solver<T, gmm::mr_approx_inverse_precond> s_t;
    //    return new s_t(new typename s_t::precond_type(), dictionary);
    //}
 
    Exception() << "Unknow GMM++ cg preconditioner " << precond_name << "." << THROW;
    return 0;
}
 
template <typename T, template <typename MATRIX> class PRECOND>
class gmm_gmres_iterative_solver : public LU_sparse_solver<T> {
public:
    typedef gmm::csc_matrix_ref<const T*, const size_t*, const size_t*> gmm_matrix_ref_t;
    typedef PRECOND<gmm_matrix_ref_t> precond_type;
 
    gmm_gmres_iterative_solver(precond_type* precond_, const Dictionary& dictionary_)
        : precond(precond_), modified_value(true), dictionary(dictionary_) {}
 
    ~gmm_gmres_iterative_solver() { delete precond; }
 
    void init(
          size_t s, size_t nnz, const size_t* colind, const size_t* rowind, const T* value,
          int connectivity, const Dictionary& dictionary_) {
        gmm_matrix_ref = gmm_matrix_ref_t(value, rowind, colind, s, s);
        update_value();
    }
 
    void update_value() { modified_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 = ' ') {
        if (modified_value) {
            precond->build_with(gmm_matrix_ref);
            modified_value = false;
        }
        gmm::iteration iter(
              dictionary.get_double("SLS_TOL_RESIDUUM", 1.0E-8),
              dictionary.get_int("SLS_VERBOSITY", 0), dictionary.get_int("SLS_MAX_ITER", -1));
        T* tmp = 0;
        if (b == x) { tmp = new T[s]; }
        size_t restart = dictionary.get_int("SLS_RESTART", 200);
        for (size_t i = 0; i != nrhs; ++i) {
            if (b == x) { std::copy(b + ldb * i, b + ldb * (i + 1), tmp); }
            gmm::array1D_reference<const T*> B(b == x ? tmp : b + ldb * i, s);
            gmm::array1D_reference<T*> X(x + ldx * i, s);
            gmm::gmres(gmm_matrix_ref, X, B, *precond, restart, iter);
            if (!iter.converged()) {
                Exception() << "The iterative sparse linear solver did not converge" << THROW;
            }
        }
        if (b == x) { delete[] tmp; }
    }
 
protected:
    gmm_matrix_ref_t gmm_matrix_ref;
    precond_type* precond;
    bool modified_value;
    const Dictionary& dictionary;
};
 
template <typename T>
LU_sparse_solver<T>* new_gmm_gmres_iterative_solver(const Dictionary& dictionary) {
    std::string precond_name = dictionary.get_string("SLS_PRECOND", "ILU");
    if (precond_name == "ILU") {
        typedef gmm_gmres_iterative_solver<T, gmm::ilu_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "ILUT") {
        typedef gmm_gmres_iterative_solver<T, gmm::ilut_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "ILUTP") {
        typedef gmm_gmres_iterative_solver<T, gmm::ilutp_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "DIAGONAL") {
        typedef gmm_gmres_iterative_solver<T, gmm::diagonal_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
    if (precond_name == "MR_APPROX_INVERSE") {
        typedef gmm_gmres_iterative_solver<T, gmm::mr_approx_inverse_precond> s_t;
        return new s_t(new typename s_t::precond_type(), dictionary);
    }
 
    Exception() << "Unknown GMM++ gmres preconditioner " << precond_name << "." << THROW;
 
    return 0;
}
 
}}  // namespace b2000::b2linalg
 
#endif