Reference documentation for deal.II version Git 082d75bebd 2019-10-16 19:44:02 +0200
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 
94 {
95 public:
100 
106  {};
107 
108 
114 
118  ~SparseDirectUMFPACK() override;
119 
131  void
132  initialize(const SparsityPattern &sparsity_pattern);
133 
151  template <class Matrix>
152  void
153  factorize(const Matrix &matrix);
154 
158  template <class Matrix>
159  void
160  initialize(const Matrix & matrix,
161  const AdditionalData additional_data = AdditionalData());
162 
183  void
184  vmult(Vector<double> &dst, const Vector<double> &src) const;
185 
189  void
190  vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
191 
196  void
197  Tvmult(Vector<double> &dst, const Vector<double> &src) const;
198 
202  void
203  Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
204 
209  size_type
210  m() const;
211 
216  size_type
217  n() const;
218 
245  void
246  solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
247 
251  void
252  solve(BlockVector<double> &rhs_and_solution,
253  const bool transpose = false) const;
254 
261  template <class Matrix>
262  void
263  solve(const Matrix & matrix,
264  Vector<double> &rhs_and_solution,
265  const bool transpose = false);
266 
270  template <class Matrix>
271  void
272  solve(const Matrix & matrix,
273  BlockVector<double> &rhs_and_solution,
274  const bool transpose = false);
275 
286  ExcUMFPACKError,
287  std::string,
288  int,
289  << "UMFPACK routine " << arg1 << " returned error status " << arg2 << "."
290  << "\n\n"
291  << ("A complete list of error codes can be found in the file "
292  "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
293  "\n\n"
294  "That said, the two most common errors that can happen are "
295  "that your matrix cannot be factorized because it is "
296  "rank deficient, and that UMFPACK runs out of memory "
297  "because your problem is too large."
298  "\n\n"
299  "The first of these cases most often happens if you "
300  "forget terms in your bilinear form necessary to ensure "
301  "that the matrix has full rank, or if your equation has a "
302  "spatially variable coefficient (or nonlinearity) that is "
303  "supposed to be strictly positive but, for whatever "
304  "reasons, is negative or zero. In either case, you probably "
305  "want to check your assembly procedure. Similarly, a "
306  "matrix can be rank deficient if you forgot to apply the "
307  "appropriate boundary conditions. For example, the "
308  "Laplace equation without boundary conditions has a "
309  "single zero eigenvalue and its rank is therefore "
310  "deficient by one."
311  "\n\n"
312  "The other common situation is that you run out of memory. "
313  "On a typical laptop or desktop, it should easily be possible "
314  "to solve problems with 100,000 unknowns in 2d. If you are "
315  "solving problems with many more unknowns than that, in "
316  "particular if you are in 3d, then you may be running out "
317  "of memory and you will need to consider iterative "
318  "solvers instead of the direct solver employed by "
319  "UMFPACK."));
320 
321 private:
326 
331 
338  void *numeric_decomposition;
339 
343  void
344  clear();
345 
352  template <typename number>
353  void
354  sort_arrays(const SparseMatrixEZ<number> &);
355 
356  template <typename number>
357  void
358  sort_arrays(const SparseMatrix<number> &);
359 
360  template <typename number>
361  void
362  sort_arrays(const BlockSparseMatrix<number> &);
363 
367  std::vector<types::suitesparse_index> Ap;
368  std::vector<types::suitesparse_index> Ai;
369  std::vector<double> Ax;
370 
374  std::vector<double> control;
375 };
376 
377 DEAL_II_NAMESPACE_CLOSE
378 
379 #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:99
std::vector< double > control
unsigned int global_dof_index
Definition: types.h:89
std::vector< types::suitesparse_index > Ap