29#ifdef DEAL_II_WITH_UMFPACK
33#ifdef DEAL_II_WITH_MUMPS
48 template <
typename SparseMatrixType>
50 parallel_grainsize(
const SparseMatrixType &matrix)
52 const unsigned int avg_entries_per_row =
54 return std::max(1000 / avg_entries_per_row, 1u);
70#ifdef DEAL_II_WITH_UMFPACK
75 , symbolic_decomposition(nullptr)
76 , numeric_decomposition(nullptr)
77 , control(UMFPACK_CONTROL)
79 umfpack_dl_defaults(
control.data());
101 std::vector<types::suitesparse_index> tmp;
106 std::vector<types::suitesparse_index> tmp;
111 std::vector<double> tmp;
116 std::vector<double> tmp;
120 umfpack_dl_defaults(
control.data());
125template <
typename number>
140 for (size_type row = row_begin; row < row_end; ++row)
153 long int cursor = Ap[row];
154 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
156 std::swap(Ai[cursor], Ai[cursor + 1]);
158 std::swap(Ax[cursor], Ax[cursor + 1]);
159 if (numbers::NumberTraits<number>::is_complex == true)
160 std::swap(Az[cursor], Az[cursor + 1]);
166 parallel_grainsize(matrix));
171template <
typename number>
180 for (size_type row = row_begin; row < row_end; ++row)
182 long int cursor = Ap[row];
183 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
185 std::swap(Ai[cursor], Ai[cursor + 1]);
187 std::swap(Ax[cursor], Ax[cursor + 1]);
188 if (numbers::NumberTraits<number>::is_complex == true)
189 std::swap(Az[cursor], Az[cursor + 1]);
195 parallel_grainsize(matrix));
200template <
typename number>
213 for (size_type row = row_begin; row < row_end; ++row)
215 long int cursor = Ap[row];
216 for (size_type block = 0; block < matrix.n_block_cols(); ++block)
219 while ((cursor < Ap[row + 1] - 1) &&
220 (Ai[cursor] < Ai[cursor + 1]))
224 if (cursor == Ap[row + 1] - 1)
229 long int element = cursor;
230 while ((element < Ap[row + 1] - 1) &&
231 (Ai[element] > Ai[element + 1]))
233 std::swap(Ai[element], Ai[element + 1]);
235 std::swap(Ax[element], Ax[element + 1]);
236 if (numbers::NumberTraits<number>::is_complex == true)
237 std::swap(Az[element], Az[element + 1]);
244 parallel_grainsize(matrix));
249template <
class Matrix>
257 using number =
typename Matrix::value_type;
283 Ai.resize(matrix.n_nonzero_elements());
284 Ax.resize(matrix.n_nonzero_elements());
286 Az.resize(matrix.n_nonzero_elements());
291 Ap[row] =
Ap[row - 1] + matrix.get_row_length(row - 1);
301 for (size_type row = row_begin; row < row_end; ++row)
303 long int index = Ap[row];
304 for (typename Matrix::const_iterator p = matrix.begin(row);
305 p != matrix.end(row);
309 Ai[index] = p->column();
310 Ax[index] = std::real(p->value());
311 if (numbers::NumberTraits<number>::is_complex == true)
312 Az[index] = std::imag(p->value());
317 Assert(index == Ap[row + 1], ExcInternalError());
320 parallel_grainsize(matrix));
329 status = umfpack_dl_symbolic(N,
334 &symbolic_decomposition,
338 status = umfpack_zl_symbolic(N,
344 &symbolic_decomposition,
348 ExcUMFPACKError(
"umfpack_dl_symbolic", status));
351 status = umfpack_dl_numeric(Ap.data(),
354 symbolic_decomposition,
355 &numeric_decomposition,
359 status = umfpack_zl_numeric(Ap.data(),
363 symbolic_decomposition,
364 &numeric_decomposition,
369 umfpack_dl_free_symbolic(&symbolic_decomposition);
370 if (status == UMFPACK_WARNING_singular_matrix)
377 "UMFPACK reports that the matrix is singular, "
378 "but that the factorization was successful anyway. "
379 "You can try and see whether you can still "
380 "solve a linear system with such a factorization "
381 "by catching and ignoring this exception, "
382 "though in practice this will typically not "
387 ExcUMFPACKError(
"umfpack_dl_numeric", status));
402 ExcMessage(
"You have previously factored a matrix using this class "
403 "that had complex-valued entries. This then requires "
404 "applying the factored matrix to a complex-valued "
405 "vector, but you are only providing a real-valued vector "
409 rhs = rhs_and_solution;
417 const int status = umfpack_dl_solve(
transpose ? UMFPACK_A : UMFPACK_At,
421 rhs_and_solution.
begin(),
436# ifdef DEAL_II_WITH_COMPLEX_VALUES
466 for (
unsigned int i = 0; i < rhs_and_solution.size(); ++i)
468 rhs_re(i) = std::real(rhs_and_solution(i));
469 rhs_im(i) = std::imag(rhs_and_solution(i));
484 const int status = umfpack_zl_solve(
transpose ? UMFPACK_A : UMFPACK_Aat,
500 for (
unsigned int i = 0; i < rhs_and_solution.size(); ++i)
501 rhs_and_solution(i) = {solution_re(i), solution_im(i)};
516 for (
unsigned int i = 0; i < rhs.
size(); ++i)
517 rhs_real_or_imag(i) = std::real(rhs(i));
521 rhs_and_solution = rhs_real_or_imag;
526 for (
unsigned int i = 0; i < rhs.
size(); ++i)
527 rhs_real_or_imag(i) = std::imag(rhs(i));
531 for (
unsigned int i = 0; i < rhs.
size(); ++i)
532 rhs_and_solution(i).imag(rhs_real_or_imag(i));
537 (void)rhs_and_solution;
541 "This function can't be called if deal.II has been configured "
542 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
555 tmp = rhs_and_solution;
557 rhs_and_solution = tmp;
566# ifdef DEAL_II_WITH_COMPLEX_VALUES
571 tmp = rhs_and_solution;
573 rhs_and_solution = tmp;
576 (void)rhs_and_solution;
580 "This function can't be called if deal.II has been configured "
581 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
587template <
class Matrix>
599template <
class Matrix>
602 Vector<std::complex<double>> &rhs_and_solution,
605# ifdef DEAL_II_WITH_COMPLEX_VALUES
612 (void)rhs_and_solution;
616 "This function can't be called if deal.II has been configured "
617 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
623template <
class Matrix>
635template <
class Matrix>
638 BlockVector<std::complex<double>> &rhs_and_solution,
641# ifdef DEAL_II_WITH_COMPLEX_VALUES
648 (void)rhs_and_solution;
652 "This function can't be called if deal.II has been configured "
653 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
664 , symbolic_decomposition(nullptr)
665 , numeric_decomposition(nullptr)
675template <
class Matrix>
682 "To call this function you need UMFPACK, but you configured deal.II "
683 "without passing the necessary switch to 'cmake'. Please consult the "
684 "installation instructions at https://dealii.org/current/readme.html"));
694 "To call this function you need UMFPACK, but you configured deal.II "
695 "without passing the necessary switch to 'cmake'. Please consult the "
696 "installation instructions at https://dealii.org/current/readme.html"));
707 "To call this function you need UMFPACK, but you configured deal.II "
708 "without passing the necessary switch to 'cmake'. Please consult the "
709 "installation instructions at https://dealii.org/current/readme.html"));
720 "To call this function you need UMFPACK, but you configured deal.II "
721 "without passing the necessary switch to 'cmake'. Please consult the "
722 "installation instructions at https://dealii.org/current/readme.html"));
734 "To call this function you need UMFPACK, but you configured deal.II "
735 "without passing the necessary switch to 'cmake'. Please consult the "
736 "installation instructions at https://dealii.org/current/readme.html"));
741template <
class Matrix>
748 "To call this function you need UMFPACK, but you configured deal.II "
749 "without passing the necessary switch to 'cmake'. Please consult the "
750 "installation instructions at https://dealii.org/current/readme.html"));
755template <
class Matrix>
758 Vector<std::complex<double>> &,
764 "To call this function you need UMFPACK, but you configured deal.II "
765 "without passing the necessary switch to 'cmake'. Please consult the "
766 "installation instructions at https://dealii.org/current/readme.html"));
771template <
class Matrix>
778 "To call this function you need UMFPACK, but you configured deal.II "
779 "without passing the necessary switch to 'cmake'. Please consult the "
780 "installation instructions at https://dealii.org/current/readme.html"));
785template <
class Matrix>
794 "To call this function you need UMFPACK, but you configured deal.II "
795 "without passing the necessary switch to 'cmake'. Please consult the "
796 "installation instructions at https://dealii.org/current/readme.html"));
802template <
class Matrix>
833 this->
solve(dst,
true);
843 this->
solve(dst,
true);
862#ifdef DEAL_II_WITH_MUMPS
865 : initialize_called(false)
882template <
class Matrix>
898 id.comm_fortran = -987654;
910 nz =
matrix.n_actually_nonzero_elements();
927 for (
typename Matrix::const_iterator ptr =
matrix.begin(row);
956template <
class Matrix>
971 ExcMessage(
"Matrix size and rhs length must be equal."));
987 ExcMessage(
"Matrix size and solution vector length must be equal."));
989 ExcMessage(
"Class not initialized with a rhs vector."));
1001template <
class Matrix>
1056#define InstantiateUMFPACK(MatrixType) \
1057 template void SparseDirectUMFPACK::factorize(const MatrixType &); \
1058 template void SparseDirectUMFPACK::solve(const MatrixType &, \
1061 template void SparseDirectUMFPACK::solve(const MatrixType &, \
1062 Vector<std::complex<double>> &, \
1064 template void SparseDirectUMFPACK::solve(const MatrixType &, \
1065 BlockVector<double> &, \
1067 template void SparseDirectUMFPACK::solve( \
1068 const MatrixType &, BlockVector<std::complex<double>> &, const bool); \
1069 template void SparseDirectUMFPACK::initialize(const MatrixType &, \
1070 const AdditionalData)
1081#ifdef DEAL_II_WITH_COMPLEX_VALUES
1089#ifdef DEAL_II_WITH_MUMPS
1090# define InstantiateMUMPS(MATRIX) \
1091 template void SparseDirectMUMPS::initialize(const MATRIX &, \
1092 const Vector<double> &); \
1093 template void SparseDirectMUMPS::initialize(const MATRIX &);
virtual size_type size() const override
void initialize(const Matrix &matrix, const Vector< double > &rhs_vector)
void solve(Vector< double > &vector)
void initialize_matrix(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)
types::global_dof_index size_type
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
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcUMFPACKError(std::string arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInitializeAlreadyCalled()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
#define InstantiateUMFPACK(MatrixType)