deal.II version GIT relicensing-2084-gab787ef0dd 2024-11-04 03:10:00+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
sparse_direct.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_sparse_direct_h
16#define dealii_sparse_direct_h
17
18
19#include <deal.II/base/config.h>
20
23
27#include <deal.II/lac/vector.h>
28
29#ifdef DEAL_II_WITH_UMFPACK
30# include <umfpack.h>
31#endif
32
34
35namespace types
36{
41#ifdef SuiteSparse_long
42 using suitesparse_index = SuiteSparse_long;
43#else
44 using suitesparse_index = long int;
45#endif
46} // namespace types
47
92{
93public:
98
104 {};
105
106
112
116 ~SparseDirectUMFPACK() override;
117
129 void
130 initialize(const SparsityPattern &sparsity_pattern);
131
149 template <class Matrix>
150 void
151 factorize(const Matrix &matrix);
152
156 template <class Matrix>
157 void
158 initialize(const Matrix &matrix,
159 const AdditionalData additional_data = AdditionalData());
160
181 void
182 vmult(Vector<double> &dst, const Vector<double> &src) const;
183
187 void
188 vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
189
194 void
195 Tvmult(Vector<double> &dst, const Vector<double> &src) const;
196
200 void
201 Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
202
208 m() const;
209
215 n() const;
216
245 void
246 solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
247
259 void
260 solve(Vector<std::complex<double>> &rhs_and_solution,
261 const bool transpose = false) const;
262
266 void
267 solve(BlockVector<double> &rhs_and_solution,
268 const bool transpose = false) const;
269
273 void
274 solve(BlockVector<std::complex<double>> &rhs_and_solution,
275 const bool transpose = false) const;
276
283 template <class Matrix>
284 void
285 solve(const Matrix &matrix,
286 Vector<double> &rhs_and_solution,
287 const bool transpose = false);
288
292 template <class Matrix>
293 void
294 solve(const Matrix &matrix,
295 Vector<std::complex<double>> &rhs_and_solution,
296 const bool transpose = false);
297
301 template <class Matrix>
302 void
303 solve(const Matrix &matrix,
304 BlockVector<double> &rhs_and_solution,
305 const bool transpose = false);
306
310 template <class Matrix>
311 void
312 solve(const Matrix &matrix,
313 BlockVector<std::complex<double>> &rhs_and_solution,
314 const bool transpose = false);
315
327 std::string,
328 int,
329 << "UMFPACK routine " << arg1 << " returned error status " << arg2 << '.'
330 << "\n\n"
331 << ("A complete list of error codes can be found in the file "
332 "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
333 "\n\n"
334 "That said, the two most common errors that can happen are "
335 "that your matrix cannot be factorized because it is "
336 "rank deficient, and that UMFPACK runs out of memory "
337 "because your problem is too large."
338 "\n\n"
339 "The first of these cases most often happens if you "
340 "forget terms in your bilinear form necessary to ensure "
341 "that the matrix has full rank, or if your equation has a "
342 "spatially variable coefficient (or nonlinearity) that is "
343 "supposed to be strictly positive but, for whatever "
344 "reasons, is negative or zero. In either case, you probably "
345 "want to check your assembly procedure. Similarly, a "
346 "matrix can be rank deficient if you forgot to apply the "
347 "appropriate boundary conditions. For example, the "
348 "Laplace equation for a problem where only Neumann boundary "
349 "conditions are posed (or where you forget to apply Dirichlet "
350 "boundary conditions) has exactly one eigenvalue equal to zero "
351 "and its rank is therefore deficient by one. Finally, the matrix "
352 "may be rank deficient because you are using a quadrature "
353 "formula with too few quadrature points."
354 "\n\n"
355 "The other common situation is that you run out of memory. "
356 "On a typical laptop or desktop, it should easily be possible "
357 "to solve problems with 100,000 unknowns in 2d. If you are "
358 "solving problems with many more unknowns than that, in "
359 "particular if you are in 3d, then you may be running out "
360 "of memory and you will need to consider iterative "
361 "solvers instead of the direct solver employed by "
362 "UMFPACK."));
363
364private:
369
375
383
387 void
388 clear();
389
396 template <typename number>
397 void
399
400 template <typename number>
401 void
403
404 template <typename number>
405 void
407
419 std::vector<types::suitesparse_index> Ap;
420 std::vector<types::suitesparse_index> Ai;
421 std::vector<double> Ax;
422 std::vector<double> Az;
423
427 std::vector<double> control;
428};
429
431
432#endif // dealii_sparse_direct_h
std::vector< double > Az
~SparseDirectUMFPACK() override
void initialize(const SparsityPattern &sparsity_pattern)
void Tvmult(Vector< double > &dst, const Vector< double > &src) const
size_type m() const
void solve(Vector< double > &rhs_and_solution, const bool transpose=false) const
std::vector< double > Ax
void sort_arrays(const SparseMatrixEZ< number > &)
size_type n() const
void factorize(const Matrix &matrix)
std::vector< double > control
void vmult(Vector< double > &dst, const Vector< double > &src) const
std::vector< types::suitesparse_index > Ap
std::vector< types::suitesparse_index > Ai
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
static ::ExceptionBase & ExcUMFPACKError(std::string arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:534
Definition types.h:32
unsigned int global_dof_index
Definition types.h:81
long int suitesparse_index