deal.II version GIT relicensing-3512-g0a98d4ed9f 2025-06-14 14: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
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
52#ifdef DEAL_II_WITH_MUMPS
53 using mumps_index = MUMPS_INT;
54 using mumps_nnz = MUMPS_INT8;
55#else
57 using mumps_nnz = std::size_t;
58#endif
59
60
61} // namespace types
62
107{
108public:
113
119 {};
120
121
127
131 ~SparseDirectUMFPACK() override;
132
144 void
145 initialize(const SparsityPattern &sparsity_pattern);
146
164 template <class Matrix>
165 void
166 factorize(const Matrix &matrix);
167
171 template <class Matrix>
172 void
173 initialize(const Matrix &matrix,
174 const AdditionalData additional_data = AdditionalData());
175
196 void
197 vmult(Vector<double> &dst, const Vector<double> &src) const;
198
202 void
203 vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
204
209 void
210 Tvmult(Vector<double> &dst, const Vector<double> &src) const;
211
215 void
216 Tvmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
217
223 m() const;
224
230 n() const;
231
260 void
261 solve(Vector<double> &rhs_and_solution, const bool transpose = false) const;
262
274 void
275 solve(Vector<std::complex<double>> &rhs_and_solution,
276 const bool transpose = false) const;
277
281 void
282 solve(BlockVector<double> &rhs_and_solution,
283 const bool transpose = false) const;
284
288 void
289 solve(BlockVector<std::complex<double>> &rhs_and_solution,
290 const bool transpose = false) const;
291
298 template <class Matrix>
299 void
300 solve(const Matrix &matrix,
301 Vector<double> &rhs_and_solution,
302 const bool transpose = false);
303
307 template <class Matrix>
308 void
309 solve(const Matrix &matrix,
310 Vector<std::complex<double>> &rhs_and_solution,
311 const bool transpose = false);
312
316 template <class Matrix>
317 void
318 solve(const Matrix &matrix,
319 BlockVector<double> &rhs_and_solution,
320 const bool transpose = false);
321
325 template <class Matrix>
326 void
327 solve(const Matrix &matrix,
328 BlockVector<std::complex<double>> &rhs_and_solution,
329 const bool transpose = false);
330
342 std::string,
343 int,
344 << "UMFPACK routine " << arg1 << " returned error status " << arg2 << '.'
345 << "\n\n"
346 << ("A complete list of error codes can be found in the file "
347 "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
348 "\n\n"
349 "That said, the two most common errors that can happen are "
350 "that your matrix cannot be factorized because it is "
351 "rank deficient, and that UMFPACK runs out of memory "
352 "because your problem is too large."
353 "\n\n"
354 "The first of these cases most often happens if you "
355 "forget terms in your bilinear form necessary to ensure "
356 "that the matrix has full rank, or if your equation has a "
357 "spatially variable coefficient (or nonlinearity) that is "
358 "supposed to be strictly positive but, for whatever "
359 "reasons, is negative or zero. In either case, you probably "
360 "want to check your assembly procedure. Similarly, a "
361 "matrix can be rank deficient if you forgot to apply the "
362 "appropriate boundary conditions. For example, the "
363 "Laplace equation for a problem where only Neumann boundary "
364 "conditions are posed (or where you forget to apply Dirichlet "
365 "boundary conditions) has exactly one eigenvalue equal to zero "
366 "and its rank is therefore deficient by one. Finally, the matrix "
367 "may be rank deficient because you are using a quadrature "
368 "formula with too few quadrature points."
369 "\n\n"
370 "The other common situation is that you run out of memory. "
371 "On a typical laptop or desktop, it should easily be possible "
372 "to solve problems with 100,000 unknowns in 2d. If you are "
373 "solving problems with many more unknowns than that, in "
374 "particular if you are in 3d, then you may be running out "
375 "of memory and you will need to consider iterative "
376 "solvers instead of the direct solver employed by "
377 "UMFPACK."));
378
379private:
384
390
398
402 void
403 clear();
404
411 template <typename number>
412 void
414
415 template <typename number>
416 void
418
419 template <typename number>
420 void
422
434 std::vector<types::suitesparse_index> Ap;
435 std::vector<types::suitesparse_index> Ai;
436 std::vector<double> Ax;
437 std::vector<double> Az;
438
442 std::vector<double> control;
443};
444
445
446
465{
466public:
471
472
478 {
483 {
484 BlockLowRank(const bool blr_ucfs = false,
485 const double lowrank_threshold = 1e-8)
488 {}
494
499 };
500
501 AdditionalData(const bool output_details = false,
502 const bool error_statistics = false,
503 const bool symmetric = false,
504 const bool posdef = false,
505 const bool blr_factorization = false,
506 const BlockLowRank &blr = BlockLowRank())
510 , posdef(posdef)
512 , blr(blr)
513 {}
514
515 /*
516 * If true, the MUMPS solver will print out details of the execution.
517 */
519
520 /*
521 * If true, the MUMPS solver will print out error statistics.
522 */
524
525 /*
526 * If true, the MUMPS solver will use the symmetric factorization. This is
527 * only possible if the matrix is symmetric.
528 */
530
531 /*
532 * If true, the MUMPS solver will use the positive definite factorization.
533 * This is only possible if the matrix is symmetric and positive definite.
534 */
535 bool posdef;
536
537 /*
538 * If true, the MUMPS solver will use the Block Low-Rank factorization.
539 */
541
542 /*
543 * Stores Block Low-Rank approximation settings to be used in MUMPS
544 * factorization, if enabled.
545 */
547 };
548
553 const MPI_Comm &communicator = MPI_COMM_WORLD);
554
559
564
570 template <class Matrix>
571 void
572 initialize(const Matrix &matrix);
573
579 void
580 vmult(Vector<double> &dst, const Vector<double> &src) const;
581
587 void
588 Tvmult(Vector<double> &dst, const Vector<double> &src) const;
589
601 int *
603
604private:
605#ifdef DEAL_II_WITH_MUMPS
606 mutable DMUMPS_STRUC_C id;
607
608#endif // DEAL_II_WITH_MUMPS
609
614 std::unique_ptr<double[]> a;
615
619 mutable std::vector<double> rhs;
620
624 std::unique_ptr<types::mumps_index[]> irn;
625
629 std::unique_ptr<types::mumps_index[]> jcn;
630
635
640
644 template <class Matrix>
645 void
646 initialize_matrix(const Matrix &matrix);
647
651 void
653
657 void
659
664
669};
670
672
673#endif // dealii_sparse_direct_h
void copy_solution(Vector< double > &vector) const
void vmult(Vector< double > &dst, const Vector< double > &src) const
void initialize_matrix(const Matrix &matrix)
void initialize(const Matrix &matrix)
std::vector< double > rhs
std::unique_ptr< double[]> a
void copy_rhs_to_mumps(const Vector< double > &rhs) const
types::global_dof_index n
const MPI_Comm mpi_communicator
std::unique_ptr< types::mumps_index[]> irn
void Tvmult(Vector< double > &dst, const Vector< double > &src) const
std::unique_ptr< types::mumps_index[]> jcn
types::mumps_nnz nnz
AdditionalData additional_data
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
std::size_t mumps_nnz
long int suitesparse_index
BlockLowRank(const bool blr_ucfs=false, const double lowrank_threshold=1e-8)
AdditionalData(const bool output_details=false, const bool error_statistics=false, const bool symmetric=false, const bool posdef=false, const bool blr_factorization=false, const BlockLowRank &blr=BlockLowRank())