sparse_direct.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2001 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_sparse_direct_h
#define dealii_sparse_direct_h
 
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/exceptions.h>
 
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_matrix_ez.h>
#include <deal.II/lac/vector.h>
 
#ifdef DEAL_II_WITH_UMFPACK
#  include <umfpack.h>
#endif
 
#ifdef DEAL_II_WITH_MUMPS
#  include <dmumps_c.h>
#endif
 
#ifdef DEAL_II_WITH_TRILINOS
#  include <deal.II/lac/trilinos_sparse_matrix.h>
#  include <deal.II/lac/trilinos_vector.h>
#endif
 
#ifdef DEAL_II_WITH_PETSC
#  include <deal.II/lac/petsc_sparse_matrix.h>
#  include <deal.II/lac/petsc_vector.h>
#endif
 
#include <complex>
 
DEAL_II_NAMESPACE_OPEN
 
namespace types
{
#ifdef SuiteSparse_long
  using suitesparse_index = SuiteSparse_long;
#else
  using suitesparse_index = long int;
#endif
 
#ifdef DEAL_II_WITH_MUMPS
  using mumps_index = MUMPS_INT;
  using mumps_nnz   = MUMPS_INT8;
#else
  using mumps_index       = int;
  using mumps_nnz         = std::size_t;
#endif
 
 
} // namespace types
 
class SparseDirectUMFPACK : public EnableObserverPointer
{
public:
  using size_type = types::global_dof_index;
 
  class AdditionalData
  {};
 
 
  SparseDirectUMFPACK();
 
  ~SparseDirectUMFPACK() override;
 
  void
  initialize(const SparsityPattern &sparsity_pattern);
 
  template <class Matrix>
  void
  factorize(const Matrix &matrix);
 
  template <class Matrix>
  void
  initialize(const Matrix        &matrix,
             const AdditionalData additional_data = AdditionalData());
 
  void
  vmult(Vector<double> &dst, const Vector<double> &src) const;
 
  void
  vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
 
  void
  Tvmult(Vector<double> &dst, const Vector<double> &src) const;
 
  void
  Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
 
  size_type
  m() const;
 
  size_type
  n() const;
 
  void
  solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
 
  void
  solve(Vector<std::complex<double>> &rhs_and_solution,
        const bool                    transpose = false) const;
 
  void
  solve(BlockVector<double> &rhs_and_solution,
        const bool           transpose = false) const;
 
  void
  solve(BlockVector<std::complex<double>> &rhs_and_solution,
        const bool                         transpose = false) const;
 
  template <class Matrix>
  void
  solve(const Matrix   &matrix,
        Vector<double> &rhs_and_solution,
        const bool      transpose = false);
 
  template <class Matrix>
  void
  solve(const Matrix                 &matrix,
        Vector<std::complex<double>> &rhs_and_solution,
        const bool                    transpose = false);
 
  template <class Matrix>
  void
  solve(const Matrix        &matrix,
        BlockVector<double> &rhs_and_solution,
        const bool           transpose = false);
 
  template <class Matrix>
  void
  solve(const Matrix                      &matrix,
        BlockVector<std::complex<double>> &rhs_and_solution,
        const bool                         transpose = false);
 
  DeclException2(
    ExcUMFPACKError,
    std::string,
    int,
    << "UMFPACK routine " << arg1 << " returned error status " << arg2 << '.'
    << "\n\n"
    << ("A complete list of error codes can be found in the file "
        "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
        "\n\n"
        "That said, the two most common errors that can happen are "
        "that your matrix cannot be factorized because it is "
        "rank deficient, and that UMFPACK runs out of memory "
        "because your problem is too large."
        "\n\n"
        "The first of these cases most often happens if you "
        "forget terms in your bilinear form necessary to ensure "
        "that the matrix has full rank, or if your equation has a "
        "spatially variable coefficient (or nonlinearity) that is "
        "supposed to be strictly positive but, for whatever "
        "reasons, is negative or zero. In either case, you probably "
        "want to check your assembly procedure. Similarly, a "
        "matrix can be rank deficient if you forgot to apply the "
        "appropriate boundary conditions. For example, the "
        "Laplace equation for a problem where only Neumann boundary "
        "conditions are posed (or where you forget to apply Dirichlet "
        "boundary conditions) has exactly one eigenvalue equal to zero "
        "and its rank is therefore deficient by one. Finally, the matrix "
        "may be rank deficient because you are using a quadrature "
        "formula with too few quadrature points."
        "\n\n"
        "The other common situation is that you run out of memory. "
        "On a typical laptop or desktop, it should easily be possible "
        "to solve problems with 100,000 unknowns in 2d. If you are "
        "solving problems with many more unknowns than that, in "
        "particular if you are in 3d, then you may be running out "
        "of memory and you will need to consider iterative "
        "solvers instead of the direct solver employed by "
        "UMFPACK."));
 
private:
  size_type n_rows;
 
  size_type n_cols;
 
  void *symbolic_decomposition;
  void *numeric_decomposition;
 
  void
  clear();
 
  template <typename number>
  void
  sort_arrays(const SparseMatrixEZ<number> &);
 
  template <typename number>
  void
  sort_arrays(const SparseMatrix<number> &);
 
  template <typename number>
  void
  sort_arrays(const BlockSparseMatrix<number> &);
 
  std::vector<types::suitesparse_index> Ap;
  std::vector<types::suitesparse_index> Ai;
  std::vector<double>                   Ax;
  std::vector<double>                   Az;
 
  std::vector<double> control;
};
 
 
 
class SparseDirectMUMPS : public EnableObserverPointer
{
public:
  using size_type = types::global_dof_index;
 
  struct AdditionalData
  {
    struct BlockLowRank
    {
      BlockLowRank(const bool   blr_ucfs          = false,
                   const double lowrank_threshold = 1e-8)
        : blr_ucfs(blr_ucfs)
        , lowrank_threshold(lowrank_threshold)
      {}
 
 
      bool blr_ucfs;
 
      double lowrank_threshold;
    };
 
    AdditionalData(const bool          output_details    = false,
                   const bool          error_statistics  = false,
                   const bool          symmetric         = false,
                   const bool          posdef            = false,
                   const bool          blr_factorization = false,
                   const BlockLowRank &blr               = BlockLowRank())
      : output_details(output_details)
      , error_statistics(error_statistics)
      , symmetric(symmetric)
      , posdef(posdef)
      , blr_factorization(blr_factorization)
      , blr(blr)
    {}
 
    bool output_details;
 
    bool error_statistics;
 
    bool symmetric;
 
    bool posdef;
 
    bool blr_factorization;
 
    BlockLowRank blr;
  };
 
  SparseDirectMUMPS(const AdditionalData &additional_data = AdditionalData(),
                    const MPI_Comm       &communicator    = MPI_COMM_WORLD);
 
  ~SparseDirectMUMPS();
 
  DeclException0(ExcInitializeAlreadyCalled);
 
  template <class Matrix>
  void
  initialize(const Matrix &matrix);
 
  template <typename VectorType>
  void
  vmult(VectorType &dst, const VectorType &src) const;
 
 
  template <typename VectorType>
  void
  Tvmult(VectorType &, const VectorType &src) const;
 
  int *
  get_icntl();
 
private:
#ifdef DEAL_II_WITH_MUMPS
  mutable DMUMPS_STRUC_C id;
 
#endif // DEAL_II_WITH_MUMPS
 
  std::unique_ptr<double[]> a;
 
  mutable std::vector<double> rhs;
 
  mutable std::vector<types::mumps_index> irhs_loc;
 
  std::unique_ptr<types::mumps_index[]> irn;
 
  std::unique_ptr<types::mumps_index[]> jcn;
 
  types::global_dof_index n;
 
  types::mumps_nnz nnz;
 
  IndexSet locally_owned_rows;
 
  template <class Matrix>
  void
  initialize_matrix(const Matrix &matrix);
 
  void
  copy_solution(Vector<double> &vector) const;
 
  void
  copy_rhs_to_mumps(const Vector<double> &rhs) const;
 
  AdditionalData additional_data;
 
  const MPI_Comm mpi_communicator;
};
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // dealii_sparse_direct_h