Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08: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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
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
34#ifdef DEAL_II_WITH_MUMPS
35# include <dmumps_c.h>
36#endif
37
39
40namespace types
41{
46#ifdef SuiteSparse_long
47 using suitesparse_index = SuiteSparse_long;
48#else
49 using suitesparse_index = long int;
50#endif
51} // namespace types
52
97{
98public:
103
109 {};
110
111
117
121 ~SparseDirectUMFPACK() override;
122
134 void
135 initialize(const SparsityPattern &sparsity_pattern);
136
154 template <class Matrix>
155 void
156 factorize(const Matrix &matrix);
157
161 template <class Matrix>
162 void
163 initialize(const Matrix &matrix,
164 const AdditionalData additional_data = AdditionalData());
165
186 void
187 vmult(Vector<double> &dst, const Vector<double> &src) const;
188
192 void
193 vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
194
199 void
200 Tvmult(Vector<double> &dst, const Vector<double> &src) const;
201
205 void
206 Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
207
213 m() const;
214
220 n() const;
221
250 void
251 solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
252
264 void
265 solve(Vector<std::complex<double>> &rhs_and_solution,
266 const bool transpose = false) const;
267
271 void
272 solve(BlockVector<double> &rhs_and_solution,
273 const bool transpose = false) const;
274
278 void
279 solve(BlockVector<std::complex<double>> &rhs_and_solution,
280 const bool transpose = false) const;
281
288 template <class Matrix>
289 void
290 solve(const Matrix &matrix,
291 Vector<double> &rhs_and_solution,
292 const bool transpose = false);
293
297 template <class Matrix>
298 void
299 solve(const Matrix &matrix,
300 Vector<std::complex<double>> &rhs_and_solution,
301 const bool transpose = false);
302
306 template <class Matrix>
307 void
308 solve(const Matrix &matrix,
309 BlockVector<double> &rhs_and_solution,
310 const bool transpose = false);
311
315 template <class Matrix>
316 void
317 solve(const Matrix &matrix,
318 BlockVector<std::complex<double>> &rhs_and_solution,
319 const bool transpose = false);
320
332 std::string,
333 int,
334 << "UMFPACK routine " << arg1 << " returned error status " << arg2 << '.'
335 << "\n\n"
336 << ("A complete list of error codes can be found in the file "
337 "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
338 "\n\n"
339 "That said, the two most common errors that can happen are "
340 "that your matrix cannot be factorized because it is "
341 "rank deficient, and that UMFPACK runs out of memory "
342 "because your problem is too large."
343 "\n\n"
344 "The first of these cases most often happens if you "
345 "forget terms in your bilinear form necessary to ensure "
346 "that the matrix has full rank, or if your equation has a "
347 "spatially variable coefficient (or nonlinearity) that is "
348 "supposed to be strictly positive but, for whatever "
349 "reasons, is negative or zero. In either case, you probably "
350 "want to check your assembly procedure. Similarly, a "
351 "matrix can be rank deficient if you forgot to apply the "
352 "appropriate boundary conditions. For example, the "
353 "Laplace equation for a problem where only Neumann boundary "
354 "conditions are posed (or where you forget to apply Dirichlet "
355 "boundary conditions) has exactly one eigenvalue equal to zero "
356 "and its rank is therefore deficient by one. Finally, the matrix "
357 "may be rank deficient because you are using a quadrature "
358 "formula with too few quadrature points."
359 "\n\n"
360 "The other common situation is that you run out of memory. "
361 "On a typical laptop or desktop, it should easily be possible "
362 "to solve problems with 100,000 unknowns in 2d. If you are "
363 "solving problems with many more unknowns than that, in "
364 "particular if you are in 3d, then you may be running out "
365 "of memory and you will need to consider iterative "
366 "solvers instead of the direct solver employed by "
367 "UMFPACK."));
368
369private:
374
380
388
392 void
393 clear();
394
401 template <typename number>
402 void
404
405 template <typename number>
406 void
408
409 template <typename number>
410 void
412
424 std::vector<types::suitesparse_index> Ap;
425 std::vector<types::suitesparse_index> Ai;
426 std::vector<double> Ax;
427 std::vector<double> Az;
428
432 std::vector<double> control;
433};
434
435
454{
455private:
456#ifdef DEAL_II_WITH_MUMPS
457 DMUMPS_STRUC_C id;
458#endif // DEAL_II_WITH_MUMPS
459
460 double *a;
461 std::vector<double> rhs;
462 int *irn;
463 int *jcn;
466
471 template <class Matrix>
472 void
473 initialize_matrix(const Matrix &matrix);
474
478 void
480
484 void
486
492
493public:
498
503
508
513
519 template <class Matrix>
520 void
521 initialize(const Matrix &matrix, const Vector<double> &rhs_vector);
522
527 template <class Matrix>
528 void
529 initialize(const Matrix &matrix);
530
536 void
538
544 void
546};
547
549
550#endif // dealii_sparse_direct_h
void initialize(const Matrix &matrix, const Vector< double > &rhs_vector)
void solve(Vector< double > &vector)
void initialize_matrix(const Matrix &matrix)
void initialize(const Matrix &matrix)
std::vector< double > rhs
types::global_dof_index n
types::global_dof_index nz
void copy_solution(Vector< double > &vector)
void copy_rhs_to_mumps(const Vector< double > &rhs)
void vmult(Vector< double > &dst, const Vector< double > &src)
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:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DeclException0(Exception0)
static ::ExceptionBase & ExcUMFPACKError(std::string arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
static ::ExceptionBase & ExcInitializeAlreadyCalled()
Definition types.h:32
unsigned int global_dof_index
Definition types.h:94
long int suitesparse_index