Reference documentation for deal.II version Git 1c31021568 2020-01-21 10:20:35 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
sparse_direct.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_sparse_direct_h
17 #define dealii_sparse_direct_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/subscriptor.h>
24 
25 #include <deal.II/lac/block_sparse_matrix.h>
26 #include <deal.II/lac/sparse_matrix.h>
27 #include <deal.II/lac/sparse_matrix_ez.h>
28 #include <deal.II/lac/vector.h>
29 
30 #ifdef DEAL_II_WITH_UMFPACK
31 # include <umfpack.h>
32 #endif
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 namespace types
37 {
42 #ifdef SuiteSparse_long
43  using suitesparse_index = SuiteSparse_long;
44 #else
45  using suitesparse_index = long int;
46 #endif
47 } // namespace types
48 
91 {
92 public:
97 
103  {};
104 
105 
111 
115  ~SparseDirectUMFPACK() override;
116 
128  void
129  initialize(const SparsityPattern &sparsity_pattern);
130 
148  template <class Matrix>
149  void
150  factorize(const Matrix &matrix);
151 
155  template <class Matrix>
156  void
157  initialize(const Matrix & matrix,
158  const AdditionalData additional_data = AdditionalData());
159 
180  void
181  vmult(Vector<double> &dst, const Vector<double> &src) const;
182 
186  void
187  vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
188 
193  void
194  Tvmult(Vector<double> &dst, const Vector<double> &src) const;
195 
199  void
200  Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
201 
206  size_type
207  m() const;
208 
213  size_type
214  n() const;
215 
244  void
245  solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
246 
250  void
251  solve(BlockVector<double> &rhs_and_solution,
252  const bool transpose = false) const;
253 
260  template <class Matrix>
261  void
262  solve(const Matrix & matrix,
263  Vector<double> &rhs_and_solution,
264  const bool transpose = false);
265 
269  template <class Matrix>
270  void
271  solve(const Matrix & matrix,
272  BlockVector<double> &rhs_and_solution,
273  const bool transpose = false);
274 
285  ExcUMFPACKError,
286  std::string,
287  int,
288  << "UMFPACK routine " << arg1 << " returned error status " << arg2 << "."
289  << "\n\n"
290  << ("A complete list of error codes can be found in the file "
291  "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
292  "\n\n"
293  "That said, the two most common errors that can happen are "
294  "that your matrix cannot be factorized because it is "
295  "rank deficient, and that UMFPACK runs out of memory "
296  "because your problem is too large."
297  "\n\n"
298  "The first of these cases most often happens if you "
299  "forget terms in your bilinear form necessary to ensure "
300  "that the matrix has full rank, or if your equation has a "
301  "spatially variable coefficient (or nonlinearity) that is "
302  "supposed to be strictly positive but, for whatever "
303  "reasons, is negative or zero. In either case, you probably "
304  "want to check your assembly procedure. Similarly, a "
305  "matrix can be rank deficient if you forgot to apply the "
306  "appropriate boundary conditions. For example, the "
307  "Laplace equation for a problem where only Neumann boundary "
308  "conditions are posed (or where you forget to apply Dirichlet "
309  "boundary conditions) has exactly one eigenvalue equal to zero "
310  "and its rank is therefore deficient by one. Finally, the matrix "
311  "may be rank deficient because you are using a quadrature "
312  "formula with too few quadrature points."
313  "\n\n"
314  "The other common situation is that you run out of memory. "
315  "On a typical laptop or desktop, it should easily be possible "
316  "to solve problems with 100,000 unknowns in 2d. If you are "
317  "solving problems with many more unknowns than that, in "
318  "particular if you are in 3d, then you may be running out "
319  "of memory and you will need to consider iterative "
320  "solvers instead of the direct solver employed by "
321  "UMFPACK."));
322 
323 private:
328 
334 
341  void *numeric_decomposition;
342 
346  void
347  clear();
348 
355  template <typename number>
356  void
357  sort_arrays(const SparseMatrixEZ<number> &);
358 
359  template <typename number>
360  void
361  sort_arrays(const SparseMatrix<number> &);
362 
363  template <typename number>
364  void
365  sort_arrays(const BlockSparseMatrix<number> &);
366 
370  std::vector<types::suitesparse_index> Ap;
371  std::vector<types::suitesparse_index> Ai;
372  std::vector<double> Ax;
373 
377  std::vector<double> control;
378 };
379 
380 DEAL_II_NAMESPACE_CLOSE
381 
382 #endif // dealii_sparse_direct_h
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
long int suitesparse_index
Definition: sparse_direct.h:45
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
Definition: types.h:31
types::global_dof_index size_type
Definition: sparse_direct.h:96
std::vector< double > control
unsigned int global_dof_index
Definition: types.h:89
std::vector< types::suitesparse_index > Ap