Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
sparse_direct.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2020 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 
24 
28 #include <deal.II/lac/vector.h>
29 
30 #ifdef DEAL_II_WITH_UMFPACK
31 # include <umfpack.h>
32 #endif
33 
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 
258  void
259  solve(Vector<std::complex<double>> &rhs_and_solution,
260  const bool transpose = false) const;
261 
265  void
266  solve(BlockVector<double> &rhs_and_solution,
267  const bool transpose = false) const;
268 
272  void
273  solve(BlockVector<std::complex<double>> &rhs_and_solution,
274  const bool transpose = false) const;
275 
282  template <class Matrix>
283  void
284  solve(const Matrix & matrix,
285  Vector<double> &rhs_and_solution,
286  const bool transpose = false);
287 
291  template <class Matrix>
292  void
293  solve(const Matrix & matrix,
294  Vector<std::complex<double>> &rhs_and_solution,
295  const bool transpose = false);
296 
300  template <class Matrix>
301  void
302  solve(const Matrix & matrix,
303  BlockVector<double> &rhs_and_solution,
304  const bool transpose = false);
305 
309  template <class Matrix>
310  void
311  solve(const Matrix & matrix,
312  BlockVector<std::complex<double>> &rhs_and_solution,
313  const bool transpose = false);
314 
325  ExcUMFPACKError,
326  std::string,
327  int,
328  << "UMFPACK routine " << arg1 << " returned error status " << arg2 << "."
329  << "\n\n"
330  << ("A complete list of error codes can be found in the file "
331  "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
332  "\n\n"
333  "That said, the two most common errors that can happen are "
334  "that your matrix cannot be factorized because it is "
335  "rank deficient, and that UMFPACK runs out of memory "
336  "because your problem is too large."
337  "\n\n"
338  "The first of these cases most often happens if you "
339  "forget terms in your bilinear form necessary to ensure "
340  "that the matrix has full rank, or if your equation has a "
341  "spatially variable coefficient (or nonlinearity) that is "
342  "supposed to be strictly positive but, for whatever "
343  "reasons, is negative or zero. In either case, you probably "
344  "want to check your assembly procedure. Similarly, a "
345  "matrix can be rank deficient if you forgot to apply the "
346  "appropriate boundary conditions. For example, the "
347  "Laplace equation for a problem where only Neumann boundary "
348  "conditions are posed (or where you forget to apply Dirichlet "
349  "boundary conditions) has exactly one eigenvalue equal to zero "
350  "and its rank is therefore deficient by one. Finally, the matrix "
351  "may be rank deficient because you are using a quadrature "
352  "formula with too few quadrature points."
353  "\n\n"
354  "The other common situation is that you run out of memory. "
355  "On a typical laptop or desktop, it should easily be possible "
356  "to solve problems with 100,000 unknowns in 2d. If you are "
357  "solving problems with many more unknowns than that, in "
358  "particular if you are in 3d, then you may be running out "
359  "of memory and you will need to consider iterative "
360  "solvers instead of the direct solver employed by "
361  "UMFPACK."));
362 
363 private:
368 
374 
382 
386  void
387  clear();
388 
395  template <typename number>
396  void
397  sort_arrays(const SparseMatrixEZ<number> &);
398 
399  template <typename number>
400  void
401  sort_arrays(const SparseMatrix<number> &);
402 
403  template <typename number>
404  void
405  sort_arrays(const BlockSparseMatrix<number> &);
406 
418  std::vector<types::suitesparse_index> Ap;
419  std::vector<types::suitesparse_index> Ai;
420  std::vector<double> Ax;
421  std::vector<double> Az;
422 
426  std::vector<double> control;
427 };
428 
430 
431 #endif // dealii_sparse_direct_h
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
long int suitesparse_index
Definition: sparse_direct.h:45
std::vector< double > Az
Contents is actually a matrix.
std::vector< types::suitesparse_index > Ai
Definition: types.h:31
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
std::vector< double > Ax
std::vector< double > control
unsigned int global_dof_index
Definition: types.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
std::vector< types::suitesparse_index > Ap