umfpack.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/solver/base.hpp>
#include <kernel/lafem/saddle_point_matrix.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/sparse_matrix_csr.hpp>
#include <kernel/lafem/mean_filter.hpp>
#include <kernel/lafem/tuple_vector.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
 
namespace FEAT
{
  namespace Solver
  {
#if defined(FEAT_HAVE_UMFPACK) || defined(DOXYGEN)
    class Umfpack :
      public SolverBase<LAFEM::DenseVector<double, Index>>
    {
    public:
      typedef LAFEM::SparseMatrixCSR<double, Index> MatrixType;
      typedef LAFEM::DenseVector<double, Index> VectorType;
 
      typedef SolverBase<VectorType> BaseClass;
 
    private:
      const MatrixType& _system_matrix;
      double* _umf_control;
      void* _umf_symbolic;
      void* _umf_numeric;
 
      std::size_t _sym_peak_size;
      std::size_t _sym_mem_size;
      std::size_t _num_mem_size;
      std::size_t _umf_peak_size;
 
    public:
      explicit Umfpack(const MatrixType& system_matrix);
 
      virtual ~Umfpack();
 
      virtual String name() const override
      {
        return "Umfpack";
      }
 
      virtual void init_symbolic() override;
      virtual void done_symbolic() override;
      virtual void init_numeric() override;
      virtual void done_numeric() override;
 
      virtual Status apply(VectorType& vec_sol, const VectorType& vec_rhs) override;
    }; // class Umfpack
 
    inline std::shared_ptr<Umfpack> new_umfpack(const LAFEM::SparseMatrixCSR<double, Index>& matrix)
    {
      return std::make_shared<Umfpack>(matrix);
    }
 
    class UmfpackMean :
      public SolverBase<LAFEM::DenseVector<double, Index>>
    {
    public:
      typedef LAFEM::SparseMatrixCSR<double, Index> MatrixType;
      typedef LAFEM::DenseVector<double, Index> VectorType;
 
      typedef SolverBase<VectorType> BaseClass;
 
    private:
      const MatrixType& _system_matrix;
      const VectorType& _weight_vector;
      MatrixType _solver_matrix;
      VectorType _vec_x, _vec_b;
      Umfpack _umfpack;
 
    public:
      explicit UmfpackMean(const MatrixType& system_matrix, const VectorType& weight_vector);
 
      explicit UmfpackMean(
        const MatrixType& system_matrix,
        const LAFEM::MeanFilter<double, Index>& mean_filter) :
        UmfpackMean(system_matrix, mean_filter.get_vec_dual())
      {
      }
 
      virtual String name() const override
      {
        return "UmfpackMean";
      }
 
      virtual void init_symbolic() override;
      virtual void done_symbolic() override;
      virtual void init_numeric() override;
      virtual void done_numeric() override;
 
      virtual Status apply(VectorType& vec_sol, const VectorType& vec_rhs) override;
    }; // class UmfpackMean
 
    inline std::shared_ptr<UmfpackMean> new_umfpack_mean(
      const LAFEM::SparseMatrixCSR<double, Index>& matrix,
      const LAFEM::DenseVector<double, Index>& weight_vector)
    {
      return std::make_shared<UmfpackMean>(matrix, weight_vector);
    }
 
    inline std::shared_ptr<UmfpackMean> new_umfpack_mean(
      const LAFEM::SparseMatrixCSR<double, Index>& matrix,
      const LAFEM::MeanFilter<double, Index>& filter)
    {
      return std::make_shared<UmfpackMean>(matrix, filter);
    }
 
    template<typename DT_, typename IT_, int dim_>
    class SaddleUmfpackMean :
      public SolverBase<LAFEM::TupleVector<LAFEM::DenseVectorBlocked<DT_, IT_, dim_>, LAFEM::DenseVector<DT_, IT_>>>
    {
    public:
      typedef LAFEM::SparseMatrixBCSR<DT_, IT_, dim_, dim_> MatrixTypeA;
      typedef LAFEM::SparseMatrixBCSR<DT_, IT_, dim_, 1> MatrixTypeB;
      typedef LAFEM::SparseMatrixBCSR<DT_, IT_, 1, dim_> MatrixTypeD;
      typedef LAFEM::SaddlePointMatrix<MatrixTypeA, MatrixTypeB, MatrixTypeD> MatrixType;
 
      typedef LAFEM::DenseVectorBlocked<DT_, IT_, dim_> VectorTypeV;
      typedef LAFEM::DenseVector<DT_, IT_> VectorTypeP;
      typedef LAFEM::TupleVector<VectorTypeV, VectorTypeP> VectorType;
 
      typedef LAFEM::SparseMatrixCSR<double, Index> UmfMatrixType;
      typedef LAFEM::DenseVector<double, Index> UmfVectorType;
 
      typedef SolverBase<VectorType> BaseClass;
 
    private:
      const MatrixType& _system_matrix;
      const VectorTypeP& _weight_vector;
      UmfMatrixType _solver_matrix;
      UmfVectorType _vec_x, _vec_b;
      Umfpack _umfpack;
 
    public:
      explicit SaddleUmfpackMean(const MatrixType& system_matrix, const VectorTypeP& weight_vector) :
        BaseClass(),
        _system_matrix(system_matrix),
        _weight_vector(weight_vector),
        _umfpack(_solver_matrix)
      {
      }
 
      explicit SaddleUmfpackMean(
        const MatrixType& system_matrix,
        const LAFEM::MeanFilter<DT_, IT_>& mean_filter) :
        SaddleUmfpackMean(system_matrix, mean_filter.get_vec_dual())
      {
      }
 
      virtual String name() const override
      {
        return "SaddleUmfpackMean";
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
 
        // get sub-matrices
        const MatrixTypeA& matrix_a = _system_matrix.block_a();
        const MatrixTypeB& matrix_b = _system_matrix.block_b();
        const MatrixTypeD& matrix_d = _system_matrix.block_d();
 
        // get velocity/pressure space dimensions
        const Index n_v = matrix_b.rows();
        const Index n_p = matrix_b.columns();
 
        // verify matrix dimensions
        XASSERT(n_v == matrix_a.rows());
        XASSERT(n_v == matrix_a.columns());
        XASSERT(n_p == matrix_d.rows());
        XASSERT(n_v == matrix_d.columns());
        XASSERTM(n_p == _weight_vector.size(), "invalid weight vector/mean filter size");
 
        // get the number of non-zeroes in sub-matrices
        const Index nnze_a = matrix_a.used_elements();
        const Index nnze_b = matrix_b.used_elements();
        const Index nnze_d = matrix_d.used_elements();
 
        // allocate our solver matrix
        _solver_matrix = UmfMatrixType(dim_*n_v+n_p+1, dim_*n_v+n_p+1, dim_*dim_*nnze_a + dim_*(nnze_b+nnze_d) + 2*n_p);
 
        // get our input matrix arrays
        const IT_* irow_ptr_a = matrix_a.row_ptr();
        const IT_* icol_idx_a = matrix_a.col_ind();
        const IT_* irow_ptr_b = matrix_b.row_ptr();
        const IT_* icol_idx_b = matrix_b.col_ind();
        const IT_* irow_ptr_d = matrix_d.row_ptr();
        const IT_* icol_idx_d = matrix_d.col_ind();
 
        // get our output matrix arrays
        Index* orow_ptr = _solver_matrix.row_ptr();
        Index* ocol_idx = _solver_matrix.col_ind();
 
        // assemble the solver matrix structure
        orow_ptr[0] = Index(0);
        // loop over all rows of A and B
        for(Index i(0); i < n_v; ++i)
        {
          Index id = i*Index(dim_);
          // loop over all block dimensions
          for(Index j(0); j < Index(dim_); ++j, ++id)
          {
            Index op = orow_ptr[id];
            // copy input matrix A row pattern
            for(Index ip(irow_ptr_a[i]); ip < irow_ptr_a[i+1]; ++ip)
            {
              // convert block to scalar indices
              for(Index k(0); k < Index(dim_); ++k, ++op)
                ocol_idx[op] = Index(icol_idx_a[ip])*Index(dim_) + k;
            }
            // copy input matrix B row pattern
            for(Index ip(irow_ptr_b[i]); ip < irow_ptr_b[i+1]; ++ip, ++op)
              ocol_idx[op] = n_v*Index(dim_) + Index(icol_idx_b[ip]);
            // set next row pointer
            orow_ptr[id+1] = op;
          }
        }
        // loop over all rows of D
        for(Index i(0); i < n_p; ++i)
        {
          Index op = orow_ptr[Index(dim_)*n_v + i];
          // copy input matrix D row pattern
          for(Index ip(irow_ptr_d[i]); ip < irow_ptr_d[i+1]; ++ip)
          {
            // convert block to scalar value
            for(Index k(0); k < Index(dim_); ++k, ++op)
              ocol_idx[op] = Index(icol_idx_d[ip])*Index(dim_) + k;
          }
          // append a single entry in the last column
          ocol_idx[op] = Index(dim_)*n_v + n_p;
          // set next row pointer
          orow_ptr[Index(dim_)*n_v+i+1] = ++op;
        }
 
        // append an (almost) dense row for pressure
        Index op(orow_ptr[Index(dim_)*n_v + n_p]);
        for(Index j(0); j < n_p; ++j, ++op)
        {
          ocol_idx[op] = Index(dim_)*n_v + j;
        }
        // set last row pointer
        orow_ptr[Index(dim_)*n_v+n_p+1] = op;
 
        // Okay, solver matrix structure assembly complete
 
        // create temporary vectors
        _vec_x = _solver_matrix.create_vector_r();
        _vec_b = _solver_matrix.create_vector_r();
 
        // initialize umfpack
        _umfpack.init_symbolic();
      }
 
      virtual void done_symbolic() override
      {
        _umfpack.done_symbolic();
        _vec_b.clear();
        _vec_x.clear();
        _solver_matrix.clear();
 
        BaseClass::done_symbolic();
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        // get sub-matrices
        const MatrixTypeA& matrix_a = _system_matrix.block_a();
        const MatrixTypeB& matrix_b = _system_matrix.block_b();
        const MatrixTypeD& matrix_d = _system_matrix.block_d();
 
        // get velocity/pressure space dimensions
        const Index n_v = matrix_b.rows();
        const Index n_p = matrix_b.columns();
 
        // get our input matrix arrays
        const IT_* irow_ptr_a = matrix_a.row_ptr();
        const auto* idata_a   = matrix_a.val();
        const IT_* irow_ptr_b = matrix_b.row_ptr();
        const auto* idata_b   = matrix_b.val();
        const IT_* irow_ptr_d = matrix_d.row_ptr();
        const auto* idata_d   = matrix_d.val();
 
        // get input vector array
        const DT_* weight = _weight_vector.elements();
 
        // get output matrix arrays
        const Index* orow_ptr = _solver_matrix.row_ptr();
        double* odata = _solver_matrix.val();
 
        // loop over all rows of A and B
        for(Index i(0); i < n_v; ++i)
        {
          Index id = i*Index(dim_);
          // loop over all block dimensions
          for(int j(0); j < dim_; ++j, ++id)
          {
            Index op = orow_ptr[id];
            // copy input matrix A row data
            for(Index ip(irow_ptr_a[i]); ip < irow_ptr_a[i+1]; ++ip)
            {
              for(int k(0); k < dim_; ++k, ++op)
                odata[op] = double(idata_a[ip][j][k]);
            }
            // copy input matrix B row data
            for(Index ip(irow_ptr_b[i]); ip < irow_ptr_b[i+1]; ++ip, ++op)
              odata[op] = double(idata_b[ip][j][0]);
          }
        }
        // loop over all rows of D
        for(Index i(0); i < n_p; ++i)
        {
          Index op = orow_ptr[Index(dim_)*n_v + i];
          // copy input matrix D row data
          for(Index ip(irow_ptr_d[i]); ip < irow_ptr_d[i+1]; ++ip)
          {
            // convert block to scalar value
            for(int k(0); k < dim_; ++k, ++op)
              odata[op] = double(idata_d[ip][0][k]);
          }
          // copy weight vector entry as last column entry
          odata[op] = double(weight[i]);
        }
 
        // copy weight vector into last row
        Index op(orow_ptr[Index(dim_)*n_v + n_p]);
        for(Index j(0); j < n_p; ++j, ++op)
        {
          odata[op] = double(weight[j]);
        }
 
        // okay, solver matrix data assembly completed
 
        // initialize umfpack
        _umfpack.init_numeric();
      }
 
      virtual void done_numeric() override
      {
        _umfpack.done_numeric();
 
        BaseClass::done_numeric();
      }
 
      virtual Status apply(VectorType& vec_sol, const VectorType& vec_rhs) override
      {
        const Index n_v = vec_sol.template at<0>().size() * Index(dim_);
        const Index n_p = vec_sol.template at<1>().size();
 
        // copy RHS vector
        const DT_* vr_v = vec_rhs.template at<0>().template elements<LAFEM::Perspective::pod>();
        const DT_* vr_p = vec_rhs.template at<1>().elements();
        double* vb = _vec_b.elements();
 
        // copy velocity RHS
        for(Index i(0); i < n_v; ++i, ++vb)
          *vb = double(vr_v[i]);
        // copy pressure RHS
        for(Index i(0); i < n_p; ++i, ++vb)
          *vb = double(vr_p[i]);
        // append Lagrange multiplier RHS
        *vb = 0.0;
 
        // solve system
        Status status = _umfpack.apply(_vec_x, _vec_b);
 
        // copy sol vector
        const double* vx = _vec_x.elements();
        DT_* vs_v = vec_sol.template at<0>().template elements<LAFEM::Perspective::pod>();
        DT_* vs_p = vec_sol.template at<1>().elements();
        // copy velocity solution
        for(Index i(0); i < n_v; ++i, ++vx)
          vs_v[i] = DT_(*vx);
        // copy pressure solution
        for(Index i(0); i < n_p; ++i, ++vx)
          vs_p[i] = DT_(*vx);
 
        // okay
        return status;
      }
    }; // class SaddleUmfpackMean
 
    template<typename DT_, typename IT_, int dim_>
    inline std::shared_ptr<SaddleUmfpackMean<DT_, IT_, dim_>> new_saddle_umfpack_mean(
      const LAFEM::SaddlePointMatrix<
        LAFEM::SparseMatrixBCSR<DT_, IT_, dim_, dim_>,
        LAFEM::SparseMatrixBCSR<DT_, IT_, dim_, 1>,
        LAFEM::SparseMatrixBCSR<DT_, IT_, 1, dim_>>& matrix,
      const LAFEM::DenseVector<DT_, IT_>& weight_vector)
    {
      return std::make_shared<SaddleUmfpackMean<DT_, IT_, dim_>>(matrix, weight_vector);
    }
 
    template<typename DT_, typename IT_, int dim_>
    inline std::shared_ptr<SaddleUmfpackMean<DT_, IT_, dim_>> new_saddle_umfpack_mean(
      const LAFEM::SaddlePointMatrix<
        LAFEM::SparseMatrixBCSR<DT_, IT_, dim_, dim_>,
        LAFEM::SparseMatrixBCSR<DT_, IT_, dim_, 1>,
        LAFEM::SparseMatrixBCSR<DT_, IT_, 1, dim_>>& matrix,
      const LAFEM::MeanFilter<DT_, IT_>& filter)
    {
      return std::make_shared<SaddleUmfpackMean<DT_, IT_, dim_>>(matrix, filter);
    }
 
    template<typename Matrix_>
    class GenericUmfpack :
      public SolverBase<typename Matrix_::VectorTypeL>
    {
    public:
      typedef SolverBase<typename Matrix_::VectorTypeL> BaseClass;
      typedef Matrix_ MatrixType;
      typedef typename MatrixType::VectorTypeL VectorType;
 
    protected:
      const MatrixType& _matrix;
      typename Umfpack::MatrixType _umf_matrix;
      typename Umfpack::VectorType _umf_vsol, _umf_vrhs;
      LAFEM::DenseVector<typename VectorType::DataType, typename VectorType::IndexType> _tmp_vsol, _tmp_vrhs;
      Umfpack _umfpack;
 
    public:
      explicit GenericUmfpack(const MatrixType& matrix) :
        _matrix(matrix),
        _umfpack(_umf_matrix)
      {
      }
 
      virtual ~GenericUmfpack()
      {
      }
 
      virtual String name() const override
      {
        return "GenericUmfpack";
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
 
        // convert our system matrix to <double, Index> (if necessary)
        typename MatrixType::template ContainerType<double, Index> mat_double;
        mat_double.convert(_matrix);
 
        // convert matrix to obtain the structure
        _umf_matrix.convert(mat_double);
 
        // create umfpack vectors
        _umf_vsol = _umf_matrix.create_vector_l();
        _umf_vrhs = _umf_matrix.create_vector_l();
 
        // convert to temporary vectors; note that these share their data arrays
        // with the umfpack vectors if their datatype is also 'double'
        _tmp_vsol.convert(_umf_vsol);
        _tmp_vrhs.convert(_umf_vrhs);
 
        // factorize symbolic
        _umfpack.init_symbolic();
      }
 
      virtual void done_symbolic() override
      {
        _umfpack.done_symbolic();
 
        _tmp_vrhs.clear();
        _tmp_vsol.clear();
        _umf_vrhs.clear();
        _umf_vsol.clear();
        _umf_matrix.clear();
 
        BaseClass::done_symbolic();
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        // convert our system matrix to <double, Index> (if necessary)
        typename MatrixType::template ContainerType<double, Index> mat_double;
        mat_double.convert(_matrix);
 
        // get the array of our CSR matrix
        const Index* row_ptr = _umf_matrix.row_ptr();
        /*const*/ Index* col_idx = _umf_matrix.col_ind();
        double* val = _umf_matrix.val();
 
        // copy entries into our CSR matrix
        for(Index i(0); i < _umf_matrix.rows(); ++i)
          mat_double.set_line(i, val + row_ptr[i], col_idx + row_ptr[i], 0);
 
        // factorize
        _umfpack.init_numeric();
      }
 
      virtual void done_numeric() override
      {
        _umfpack.done_numeric();
        BaseClass::done_numeric();
      }
 
      virtual Status apply(VectorType& vec_sol, const VectorType& vec_rhs) override
      {
        // convert RHS vector
        _tmp_vrhs.copy(vec_rhs);
        _umf_vrhs.convert(_tmp_vrhs);
 
        // solve
        Status status = _umfpack.apply(_umf_vsol, _umf_vrhs);
 
        // convert sol vector
        _tmp_vsol.convert(_umf_vsol);
        _tmp_vsol.copy_inv(vec_sol);
 
        return status;
      }
    }; // class GenericUmfpack<...>
 
    template<typename Matrix_>
    inline std::shared_ptr<GenericUmfpack<Matrix_>> new_generic_umfpack(const Matrix_& matrix)
    {
      return std::make_shared<GenericUmfpack<Matrix_>> (matrix);
    }
 
#endif // defined(FEAT_HAVE_UMFPACK) || defined(DOXYGEN)
  } // namespace Solver
} // namespace FEAT