sparse_direct.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2001 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#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
 
DEAL_II_NAMESPACE_OPEN
 
namespace types
{
#ifdef SuiteSparse_long
  using suitesparse_index = SuiteSparse_long;
#else
  using suitesparse_index = long int;
#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;
};
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // dealii_sparse_direct_h