15#ifndef dealii_sparse_direct_h
16#define dealii_sparse_direct_h
30#ifdef DEAL_II_WITH_UMFPACK
34#ifdef DEAL_II_WITH_MUMPS
46#ifdef SuiteSparse_long
52#ifdef DEAL_II_WITH_MUMPS
164 template <
class Matrix>
171 template <
class Matrix>
275 solve(
Vector<std::complex<double>> &rhs_and_solution,
298 template <
class Matrix>
300 solve(
const Matrix &matrix,
307 template <
class Matrix>
309 solve(
const Matrix &matrix,
310 Vector<std::complex<double>> &rhs_and_solution,
316 template <
class Matrix>
318 solve(
const Matrix &matrix,
325 template <
class Matrix>
327 solve(
const Matrix &matrix,
328 BlockVector<std::complex<double>> &rhs_and_solution,
344 <<
"UMFPACK routine " << arg1 <<
" returned error status " << arg2 <<
'.'
346 << (
"A complete list of error codes can be found in the file "
347 "<bundled/umfpack/UMFPACK/Include/umfpack.h>."
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."
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."
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 "
411 template <
typename number>
415 template <
typename number>
419 template <
typename number>
434 std::vector<types::suitesparse_index>
Ap;
435 std::vector<types::suitesparse_index>
Ai;
436 std::vector<double>
Ax;
437 std::vector<double>
Az;
504 const bool posdef =
false,
553 const MPI_Comm &communicator = MPI_COMM_WORLD);
570 template <
class Matrix>
605#ifdef DEAL_II_WITH_MUMPS
606 mutable DMUMPS_STRUC_C id;
614 std::unique_ptr<double[]>
a;
619 mutable std::vector<double>
rhs;
624 std::unique_ptr<types::mumps_index[]>
irn;
629 std::unique_ptr<types::mumps_index[]>
jcn;
644 template <
class Matrix>
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
AdditionalData additional_data
void * numeric_decomposition
~SparseDirectUMFPACK() override
void initialize(const SparsityPattern &sparsity_pattern)
void Tvmult(Vector< double > &dst, const Vector< double > &src) const
void * symbolic_decomposition
void solve(Vector< double > &rhs_and_solution, const bool transpose=false) const
void sort_arrays(const SparseMatrixEZ< number > &)
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
#define DEAL_II_NAMESPACE_CLOSE
#define DeclException0(Exception0)
static ::ExceptionBase & ExcUMFPACKError(std::string arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
static ::ExceptionBase & ExcInitializeAlreadyCalled()
unsigned int global_dof_index
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())