Reference documentation for deal.II version Git 28ce845196 2019-06-26 10:31:03 +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 #include <deal.II/base/thread_management.h>
25 
26 #include <deal.II/lac/block_sparse_matrix.h>
27 #include <deal.II/lac/sparse_matrix.h>
28 #include <deal.II/lac/sparse_matrix_ez.h>
29 #include <deal.II/lac/vector.h>
30 
31 #ifdef DEAL_II_WITH_UMFPACK
32 # include <umfpack.h>
33 #endif
34 
35 DEAL_II_NAMESPACE_OPEN
36 
37 namespace types
38 {
43 #ifdef SuiteSparse_long
44  using suitesparse_index = SuiteSparse_long;
45 #else
46  using suitesparse_index = long int;
47 #endif
48 } // namespace types
49 
95 {
96 public:
101 
107  {};
108 
109 
115 
119  ~SparseDirectUMFPACK() override;
120 
132  void
133  initialize(const SparsityPattern &sparsity_pattern);
134 
152  template <class Matrix>
153  void
154  factorize(const Matrix &matrix);
155 
159  template <class Matrix>
160  void
161  initialize(const Matrix & matrix,
162  const AdditionalData additional_data = AdditionalData());
163 
184  void
185  vmult(Vector<double> &dst, const Vector<double> &src) const;
186 
190  void
191  vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
192 
197  void
198  Tvmult(Vector<double> &dst, const Vector<double> &src) const;
199 
203  void
204  Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
205 
210  size_type
211  m() const;
212 
217  size_type
218  n() const;
219 
246  void
247  solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
248 
252  void
253  solve(BlockVector<double> &rhs_and_solution,
254  const bool transpose = false) const;
255 
262  template <class Matrix>
263  void
264  solve(const Matrix & matrix,
265  Vector<double> &rhs_and_solution,
266  const bool transpose = false);
267 
271  template <class Matrix>
272  void
273  solve(const Matrix & matrix,
274  BlockVector<double> &rhs_and_solution,
275  const bool transpose = false);
276 
287  ExcUMFPACKError,
288  std::string,
289  int,
290  << "UMFPACK routine " << arg1 << " returned error status " << arg2 << "."
291  << "\n\n"
292  << ("A complete list of error codes can be found in the file "
293  "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
294  "\n\n"
295  "That said, the two most common errors that can happen are "
296  "that your matrix cannot be factorized because it is "
297  "rank deficient, and that UMFPACK runs out of memory "
298  "because your problem is too large."
299  "\n\n"
300  "The first of these cases most often happens if you "
301  "forget terms in your bilinear form necessary to ensure "
302  "that the matrix has full rank, or if your equation has a "
303  "spatially variable coefficient (or nonlinearity) that is "
304  "supposed to be strictly positive but, for whatever "
305  "reasons, is negative or zero. In either case, you probably "
306  "want to check your assembly procedure. Similarly, a "
307  "matrix can be rank deficient if you forgot to apply the "
308  "appropriate boundary conditions. For example, the "
309  "Laplace equation without boundary conditions has a "
310  "single zero eigenvalue and its rank is therefore "
311  "deficient by one."
312  "\n\n"
313  "The other common situation is that you run out of memory. "
314  "On a typical laptop or desktop, it should easily be possible "
315  "to solve problems with 100,000 unknowns in 2d. If you are "
316  "solving problems with many more unknowns than that, in "
317  "particular if you are in 3d, then you may be running out "
318  "of memory and you will need to consider iterative "
319  "solvers instead of the direct solver employed by "
320  "UMFPACK."));
321 
322 private:
327 
332 
339  void *numeric_decomposition;
340 
344  void
345  clear();
346 
353  template <typename number>
354  void
355  sort_arrays(const SparseMatrixEZ<number> &);
356 
357  template <typename number>
358  void
359  sort_arrays(const SparseMatrix<number> &);
360 
361  template <typename number>
362  void
363  sort_arrays(const BlockSparseMatrix<number> &);
364 
368  std::vector<types::suitesparse_index> Ap;
369  std::vector<types::suitesparse_index> Ai;
370  std::vector<double> Ax;
371 
375  std::vector<double> control;
376 };
377 
378 DEAL_II_NAMESPACE_CLOSE
379 
380 #endif // dealii_sparse_direct_h
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
long int suitesparse_index
Definition: sparse_direct.h:46
Definition: types.h:31
types::global_dof_index size_type
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
std::vector< double > control
unsigned int global_dof_index
Definition: types.h:89
std::vector< types::suitesparse_index > Ap