sparse_direct.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2016 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/base/thread_management.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_matrix_ez.h>
#include <deal.II/lac/block_sparse_matrix.h>

#ifdef DEAL_II_WITH_UMFPACK
#  include <umfpack.h>
#endif
#ifndef SuiteSparse_long
#  define SuiteSparse_long long int
#endif

DEAL_II_NAMESPACE_OPEN

class SparseDirectUMFPACK : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;

  class AdditionalData
  {};


  SparseDirectUMFPACK ();

  ~SparseDirectUMFPACK ();

  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, bool transpose = false) const;

  void solve (BlockVector<double> &rhs_and_solution, bool transpose = false) const;

  template <class Matrix>
  void solve (const Matrix   &matrix,
              Vector<double> &rhs_and_solution,
              bool            transpose = false);

  template <class Matrix>
  void solve (const Matrix        &matrix,
              BlockVector<double> &rhs_and_solution,
              bool                 transpose = false);

  DeclException2 (ExcUMFPACKError, char *, 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 without boundary conditions has a "
                      "single zero eigenvalue and its rank is therefore "
                      "deficient by one."
                      "\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 _m;

  size_type _n;

  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<SuiteSparse_long> Ap;
  std::vector<SuiteSparse_long> Ai;
  std::vector<double> Ax;

  std::vector<double> control;
};

DEAL_II_NAMESPACE_CLOSE

#endif // dealii__sparse_direct_h