sparse_direct.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2020 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.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/exceptions.h>
#include <deal.II/base/subscriptor.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 Subscriptor
{
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