b2umfpack_solver.H

//------------------------------------------------------------------------
// b2umfpack_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 __B2UMFPACK_SOLVER_H__
#define __B2UMFPACK_SOLVER_H__
 
#include "b2blaslapack.H"
#include "b2ppconfig.h"
#include "b2sparse_solver.H"
 
extern "C" {
#ifdef HAVE_UMFPACK_UMFPACK_H
#include "umfpack/umfpack.h"
#elif defined HAVE_UFSPARSE_UMFPACK_H
#include "ufsparse/umfpack.h"
#elif defined HAVE_SUITESPARSE_UMFPACK_H
#include "suitesparse/umfpack.h"
#elif defined HAVE_UMFPACK_H
#include "umfpack.h"
#endif
}
 
#ifdef HAVE_UMFPACK
 
namespace b2000 { namespace b2linalg {
 
template <typename T>
class UMFPACK_LU_sparse_direct_solver : public LU_sparse_solver<T> {
public:
    UMFPACK_LU_sparse_direct_solver()
        : Symbolic(0),
          Numeric(0),
          s(0),
          colind(0),
          rowind(0),
          value(0),
          logger(logging::get_logger("linear_algebra.sparse_lu_solver.umfpack")) {
        // set defaults is used for the Umfpack call, which is different
        // depending on the data type (double or complex)
        set_defaults();
        logger.LOG(logging::info, "Using UMFPACK as linear solver! ");
        Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM;
        if (logger.is_enabled_for(logging::debug)) { Control[UMFPACK_PRL] = 2; }
    }
 
    ~UMFPACK_LU_sparse_direct_solver() {
        /* The set_free_linker function allows to use the
         * umfpack call to free Symbolic and Numeric, depending on the
         * data type.*/
        set_free_linker();
    }
 
    /* Initialization function is the same for all data types.
       The only difference is the umfpack free call for Symbolic and Numeric,
       which is different for double and complex. The set_free_linker function
       takes care of this difference. The correct call will be
       included later as required. */
    void init(
          size_t s_, size_t nnz_, const size_t* colind_, const size_t* rowind_, const T* value_,
          int connectivity, const Dictionary& dictionary) {
        s = s_;
        if (s == 0) { return; }
        colind = colind_;
        rowind = rowind_;
        value = value_;
        set_free_linker();
        Symbolic = 0;
        Numeric = 0;
    }
 
    // Function to update the values
    void update_value();
 
    /*Function to resolve the linear system, valid for complex and
    csda data type. */
    void resolve(
          size_t s_, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx,
          char left_or_right = ' ') {
        if (s == 0) { return; }
        if (!Symbolic) {
            int status = umfpack_zl_symbolic(
                  s_, s_, reinterpret_cast<long*>(const_cast<size_t*>(colind)),
                  reinterpret_cast<long*>(const_cast<size_t*>(rowind)),
                  reinterpret_cast<double*>(const_cast<T*>(value)), 0, &Symbolic, Control, info);
            if (status != UMFPACK_OK) {
                Exception() << "umfpack_zl_symbolic return the error code " << status << ": "
                            << get_error(status) << "." << THROW;
            }
        }
        if (!Numeric) {
            int status = umfpack_zl_numeric(
                  reinterpret_cast<long*>(const_cast<size_t*>(colind)),
                  reinterpret_cast<long*>(const_cast<size_t*>(rowind)),
                  reinterpret_cast<double*>(const_cast<T*>(value)), 0, Symbolic, &Numeric, Control,
                  info);
            if (status < UMFPACK_OK) {
                Exception() << "umfpack_zl_numeric return the error code " << status << ": "
                            << get_error(status) << "." << THROW;
            }
        }
        if (left_or_right != ' ') { UnimplementedError() << THROW; }
 
        if (b != x) {
            for (size_t i = 0; i != nrhs; ++i) {
                int status = umfpack_zl_solve(
                      UMFPACK_A, reinterpret_cast<long*>(const_cast<size_t*>(colind)),
                      reinterpret_cast<long*>(const_cast<size_t*>(rowind)),
                      reinterpret_cast<double*>(const_cast<T*>(value)), 0,
                      reinterpret_cast<double*>(x + i * ldx), 0,
                      reinterpret_cast<double*>(const_cast<T*>(b + i * ldb)), 0, Numeric, Control,
                      info);
                if (logger.is_enabled_for(logging::info)) { umfpack_zl_report_info(Control, info); }
                if (status < UMFPACK_OK) {
                    Exception() << "umfpack_zl_solve return the error code " << status << ": "
                                << get_error(status) << "." << THROW;
                }
            }
        } else {
            auto_ptr_array<double> bxcp(new double[s_ * 2]);
            for (size_t i = 0; i != nrhs; ++i) {
                std::copy(
                      reinterpret_cast<const double*>(b + i * ldb),
                      reinterpret_cast<const double*>(b + i * ldb + s_), bxcp.get());
                int status = umfpack_zl_solve(
                      UMFPACK_A, reinterpret_cast<long*>(const_cast<size_t*>(colind)),
                      reinterpret_cast<long*>(const_cast<size_t*>(rowind)),
                      reinterpret_cast<double*>(const_cast<T*>(value)), 0,
                      reinterpret_cast<double*>(x + i * ldx), 0, bxcp.get(), 0, Numeric, Control,
                      info);
                if (logger.is_enabled_for(logging::info)) { umfpack_zl_report_info(Control, info); }
                if (status < UMFPACK_OK) {
                    Exception() << "umfpack_zl_solve return the error code " << status << ": "
                                << get_error(status) << "." << THROW;
                }
            }
        }
    }
 
protected:
    static std::string get_error(int status) {
        switch (status) {
            case UMFPACK_ERROR_out_of_memory:
                return "out of memory";
            case UMFPACK_ERROR_invalid_Numeric_object:
                return "invalid Numeric object";
            case UMFPACK_ERROR_invalid_Symbolic_object:
                return "invalid Symbolic object";
            case UMFPACK_ERROR_argument_missing:
                return "argument missing";
            case UMFPACK_ERROR_n_nonpositive:
                return "n non-positive definite";
            case UMFPACK_ERROR_invalid_matrix:
                return "invalid matrix";
            case UMFPACK_ERROR_different_pattern:
                return "different pattern";
            case UMFPACK_ERROR_invalid_system:
                return "invalid system";
            case UMFPACK_ERROR_invalid_permutation:
                return "invalid permutation";
            case UMFPACK_ERROR_internal_error:
                return "internal error";
            case UMFPACK_ERROR_file_IO:
                return "file IO error";
            default:
                return "unknown error code";
        }
    }
 
    void set_defaults();
    double Control[UMFPACK_CONTROL];
    double info[UMFPACK_INFO];
    void* Symbolic;
    void* Numeric;
    size_t s;
    const size_t* colind;
    const size_t* rowind;
    const T* value;
    logging::Logger& logger;
    // Function to free Symbolic and Numeric
    void set_free_linker();
};
 
/* Resolve for double data_type, which is defined in the C file.
   By default we define the resolve to be of type complex, as both
   data types (complex and csda) do the same operations, thus the
   function body is identical!Only in the case of data type double
   the operation differs. */
template <>
void UMFPACK_LU_sparse_direct_solver<double>::resolve(
      size_t s, size_t nrhs, const double* b, size_t lbd, double* x, size_t ldx,
      char left_or_right);
 
template <typename T>
class UMFPACK_LU_extension_sparse_direct_solver : public LU_extension_sparse_solver<T>,
                                                  public UMFPACK_LU_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, int connectivity, const Dictionary& dictionary) {
        if (size_ext != 1) { UnimplementedError() << THROW; }
        UMFPACK_LU_sparse_direct_solver<T>::init(
              size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
        modified_value = true;
    }
 
    void update_value() {
        modified_value = true;
        UMFPACK_LU_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 = ' ') {
        if (((s_ > 1 && (ma_ == 0 || mb_ == 0)) || mc_ == 0) && modified_value == true) {
            Exception() << THROW;
        }
 
        if (left_or_right != ' ') { UnimplementedError() << THROW; }
 
        if (modified_value && (ma_ != 0 && mb_ != 0 && mc_ != 0)) {
            mb.resize(s_ - 1);
            std::copy(mb_, mb_ + s_ - 1, mb.begin());
 
            d.resize((s_ - 1) * (nrhs + 1));
            {
                std::vector<T> d_tmp((s_ - 1) * (nrhs + 1));
                typename std::vector<T>::iterator iter =
                      std::copy(ma_, ma_ + s_ - 1, d_tmp.begin());
                for (size_t i = 0; i != nrhs; ++i) {
                    iter = std::copy(b + ldb * i, b + ldb * i + s_ - 1, iter);
                }
                UMFPACK_LU_sparse_direct_solver<T>::resolve(
                      s_ - 1, nrhs + 1, &d_tmp[0], s_ - 1, &d[0], s_ - 1, left_or_right);
            }
 
            div = *mc_ - blas::dot(s_ - 1, &mb[0], 1, &d[0], 1);
            if (div == T(0)) { UnimplementedError() << THROW; }
            div = 1.0 / div;
            if (s_ == 1) {
                blas::scal(nrhs, div, x, s_);
            } else {
                blas::gemv(
                      'T', s_ - 1, nrhs, -div, &d[s_ - 1], s_ - 1, &mb[0], 1, div, x + s_ - 1, s_);
                blas::ger(s_ - 1, nrhs, -1, &d[0], 1, x + s_ - 1, s_, &d[s_ - 1], s_ - 1);
            }
            for (size_t i = 0; i != nrhs; ++i) {
                std::copy(
                      d.begin() + (s_ - 1) * (i + 1), d.begin() + (s_ - 1) * (i + 2), x + ldx * i);
            }
            d.resize(s_ - 1);
            modified_value = false;
        } else {
            UMFPACK_LU_sparse_direct_solver<T>::resolve(
                  s_ - 1, nrhs, b, ldb, x, ldx, left_or_right);
            if (s_ == 1) {
                blas::scal(nrhs, div, x, s_);
            } else {
                blas::gemv('T', s_ - 1, nrhs, -div, x, ldx, &mb[0], 1, div, x + s_ - 1, s_);
                blas::ger(s_ - 1, nrhs, -1, &d[0], 1, x + s_ - 1, s_, x, ldx);
            }
        }
    }
 
private:
    bool modified_value;
    std::vector<T> d;
    std::vector<T> mb;
    T div;
};
 
template <typename T>
class UMFPACK_LDLt_sparse_direct_solver : public LDLt_sparse_solver<T> {
public:
    UMFPACK_LDLt_sparse_direct_solver()
        : s(0),
          colind(0),
          rowind(0),
          value(0),
          colind_lu(0),
          rowind_lu(0),
          value_lu(0),
          update_value_flag(true) {}
 
    ~UMFPACK_LDLt_sparse_direct_solver() {
        delete[] colind_lu;
        delete[] rowind_lu;
        delete[] value_lu;
    }
 
    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) {
        s = s_;
        colind = colind_;
        rowind = rowind_;
        value = value_;
        delete[] colind_lu;
        delete[] rowind_lu;
        delete[] value_lu;
        colind_lu = new size_t[s + 2];
        std::fill_n(colind_lu, s + 2, 0);
        colind_lu += 2;
        size_t i = colind[0];
        for (size_t j = 0; j != s; ++j) {
            for (size_t i_end = colind[j + 1]; i != i_end; ++i) {
                ++colind_lu[j];
                if (rowind[i] != j) { ++colind_lu[rowind[i]]; }
            }
        }
        --colind_lu;
        for (size_t j = 0; j != s; ++j) { colind_lu[j + 1] += colind_lu[j]; }
        rowind_lu = new size_t[colind_lu[s]];
        value_lu = new T[colind_lu[s]];
        i = colind[0];
        for (size_t j = 0; j != s; ++j) {
            for (size_t i_end = colind[j + 1]; i != i_end; ++i) {
                rowind_lu[colind_lu[j]] = rowind[i];
                value_lu[colind_lu[j]++] = value[i];
                if (rowind[i] != j) {
                    rowind_lu[colind_lu[rowind[i]]] = j;
                    value_lu[colind_lu[rowind[i]]++] = value[i];
                }
            }
        }
        --colind_lu;
 
        for (size_t j = 0; j != s; ++j) {
            if (colind_lu[j] >= colind_lu[j + 1]) {
                Exception() << "Matrix is singular in structure: dof " << j << THROW;
            }
        }
 
        solver.init(s, colind_lu[s], colind_lu, rowind_lu, value_lu, connectivity, dictionary);
        update_value_flag = false;
    }
 
    void update_value() {
        update_value_flag = true;
        solver.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 = ' ') {
        if (update_value_flag) {
            const size_t nnz = colind_lu[s];
            for (size_t j = s; j != 0; --j) { colind_lu[j] = colind_lu[j - 1]; }
            ++colind_lu;
            size_t i = colind[0];
            for (size_t j = 0; j != s; ++j) {
                for (size_t i_end = colind[j + 1]; i != i_end; ++i) {
                    value_lu[colind_lu[j]++] = value[i];
                    if (rowind[i] != j) { value_lu[colind_lu[rowind[i]]++] = value[i]; }
                }
            }
            --colind_lu;
            if (colind_lu[s] != nnz) { Exception() << THROW; }
            update_value_flag = false;
        }
        solver.resolve(s, nrhs, b, ldb, x, ldx, left_or_right);
    }
 
private:
    size_t s;
    const size_t* colind;
    const size_t* rowind;
    const T* value;
    size_t* colind_lu;
    size_t* rowind_lu;
    T* value_lu;
    bool update_value_flag;
    UMFPACK_LU_sparse_direct_solver<T> solver;
};
 
template <typename T>
class UMFPACK_LDLt_extension_sparse_direct_solver : public LDLt_extension_sparse_solver<T>,
                                                    public UMFPACK_LDLt_sparse_direct_solver<T> {
public:
    UMFPACK_LDLt_extension_sparse_direct_solver() : UMFPACK_LDLt_sparse_direct_solver<T>() {}
 
    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; }
 
        UMFPACK_LDLt_sparse_direct_solver<T>::init(
              size_, nnz_, colind_, rowind_, value_, connectivity, dictionary);
        m_ab.resize(size_ * 2);
    }
 
    void update_value() { UMFPACK_LDLt_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 = ' ') {
        if (left_or_right != ' ') { UnimplementedError() << THROW; }
 
        if (ma_ != 0 && mb_ != 0 && mc_ != 0) {
            std::copy(ma_, ma_ + s - 1, &m_ab[0]);
            std::copy(mb_, mb_ + s - 1, &m_ab[s - 1]);
            UMFPACK_LDLt_sparse_direct_solver<T>::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 T 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));
            UMFPACK_LDLt_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);
        }
    }
 
protected:
    std::vector<T> m_ab;
    T div;
};
 
}}  // namespace b2000::b2linalg
 
#endif
 
#endif