Reference documentation for deal.II version Git bed997f895 2020-09-22 11:49:20 -0400
\(\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 
88 {
89 public:
94 
100  {};
101 
102 
108 
112  ~SparseDirectUMFPACK() override;
113 
125  void
126  initialize(const SparsityPattern &sparsity_pattern);
127 
145  template <class Matrix>
146  void
147  factorize(const Matrix &matrix);
148 
152  template <class Matrix>
153  void
154  initialize(const Matrix & matrix,
155  const AdditionalData additional_data = AdditionalData());
156 
177  void
178  vmult(Vector<double> &dst, const Vector<double> &src) const;
179 
183  void
184  vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
185 
190  void
191  Tvmult(Vector<double> &dst, const Vector<double> &src) const;
192 
196  void
197  Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
198 
203  size_type
204  m() const;
205 
210  size_type
211  n() const;
212 
241  void
242  solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
243 
255  void
256  solve(Vector<std::complex<double>> &rhs_and_solution,
257  const bool transpose = false) const;
258 
262  void
263  solve(BlockVector<double> &rhs_and_solution,
264  const bool transpose = false) const;
265 
269  void
270  solve(BlockVector<std::complex<double>> &rhs_and_solution,
271  const bool transpose = false) const;
272 
279  template <class Matrix>
280  void
281  solve(const Matrix & matrix,
282  Vector<double> &rhs_and_solution,
283  const bool transpose = false);
284 
288  template <class Matrix>
289  void
290  solve(const Matrix & matrix,
291  Vector<std::complex<double>> &rhs_and_solution,
292  const bool transpose = false);
293 
297  template <class Matrix>
298  void
299  solve(const Matrix & matrix,
300  BlockVector<double> &rhs_and_solution,
301  const bool transpose = false);
302 
306  template <class Matrix>
307  void
308  solve(const Matrix & matrix,
309  BlockVector<std::complex<double>> &rhs_and_solution,
310  const bool transpose = false);
311 
322  ExcUMFPACKError,
323  std::string,
324  int,
325  << "UMFPACK routine " << arg1 << " returned error status " << arg2 << "."
326  << "\n\n"
327  << ("A complete list of error codes can be found in the file "
328  "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
329  "\n\n"
330  "That said, the two most common errors that can happen are "
331  "that your matrix cannot be factorized because it is "
332  "rank deficient, and that UMFPACK runs out of memory "
333  "because your problem is too large."
334  "\n\n"
335  "The first of these cases most often happens if you "
336  "forget terms in your bilinear form necessary to ensure "
337  "that the matrix has full rank, or if your equation has a "
338  "spatially variable coefficient (or nonlinearity) that is "
339  "supposed to be strictly positive but, for whatever "
340  "reasons, is negative or zero. In either case, you probably "
341  "want to check your assembly procedure. Similarly, a "
342  "matrix can be rank deficient if you forgot to apply the "
343  "appropriate boundary conditions. For example, the "
344  "Laplace equation for a problem where only Neumann boundary "
345  "conditions are posed (or where you forget to apply Dirichlet "
346  "boundary conditions) has exactly one eigenvalue equal to zero "
347  "and its rank is therefore deficient by one. Finally, the matrix "
348  "may be rank deficient because you are using a quadrature "
349  "formula with too few quadrature points."
350  "\n\n"
351  "The other common situation is that you run out of memory. "
352  "On a typical laptop or desktop, it should easily be possible "
353  "to solve problems with 100,000 unknowns in 2d. If you are "
354  "solving problems with many more unknowns than that, in "
355  "particular if you are in 3d, then you may be running out "
356  "of memory and you will need to consider iterative "
357  "solvers instead of the direct solver employed by "
358  "UMFPACK."));
359 
360 private:
365 
371 
379 
383  void
384  clear();
385 
392  template <typename number>
393  void
394  sort_arrays(const SparseMatrixEZ<number> &);
395 
396  template <typename number>
397  void
398  sort_arrays(const SparseMatrix<number> &);
399 
400  template <typename number>
401  void
402  sort_arrays(const BlockSparseMatrix<number> &);
403 
415  std::vector<types::suitesparse_index> Ap;
416  std::vector<types::suitesparse_index> Ai;
417  std::vector<double> Ax;
418  std::vector<double> Az;
419 
423  std::vector<double> control;
424 };
425 
427 
428 #endif // dealii_sparse_direct_h
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
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:363
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:362
std::vector< types::suitesparse_index > Ap