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
866 : additional_data(
data)
867 , mpi_communicator(communicator)
873 Assert(!(additional_data.symmetric ==
false &&
874 additional_data.posdef ==
true),
876 "You can't have a positive definite matrix that is not symmetric."));
878 if (additional_data.symmetric ==
true)
880 if (additional_data.posdef ==
true)
888 id.comm_fortran = (MUMPS_INT)MPI_Comm_c2f(mpi_communicator);
891 if (additional_data.output_details ==
false)
900 if (additional_data.error_statistics ==
true)
903 if (additional_data.blr_factorization ==
true)
906 if (additional_data.blr.blr_ucfs ==
true)
908 Assert(additional_data.blr.lowrank_threshold > 0,
909 ExcMessage(
"Lowrank threshold must be positive."));
910 id.cntl[6] = additional_data.blr.lowrank_threshold;
925template <
class Matrix>
938 nnz =
matrix.n_actually_nonzero_elements();
941 a = std::make_unique<double[]>(
nnz);
946 irn = std::make_unique<MUMPS_INT[]>(
nnz);
947 jcn = std::make_unique<MUMPS_INT[]>(
nnz);
957 for (
typename Matrix::const_iterator ptr =
matrix.begin(row);
960 if (
std::abs(ptr->value()) > 0.0 && ptr->column() >= row)
962 a[n_non_zero_elements] = ptr->value();
963 irn[n_non_zero_elements] = row + 1;
964 jcn[n_non_zero_elements] = ptr->column() + 1;
966 ++n_non_zero_elements;
974 for (
typename Matrix::const_iterator ptr =
matrix.begin(row);
979 a[n_non_zero_elements] = ptr->value();
980 irn[n_non_zero_elements] = row + 1;
981 jcn[n_non_zero_elements] = ptr->column() + 1;
982 ++n_non_zero_elements;
987 id.nnz = n_non_zero_elements;
1000 ExcMessage(
"Matrix size and rhs length must be equal."));
1006 rhs[i] = new_rhs(i);
1018 ExcMessage(
"Matrix size and solution vector length must be equal."));
1020 ExcMessage(
"Class not initialized with a rhs vector."));
1034template <
class Matrix>
1099 : mpi_communicator(MPI_COMM_SELF)
1104 "To call this function you need MUMPS, but you configured deal.II "
1105 "without passing the necessary switch to 'cmake'. Please consult the "
1106 "installation instructions at https://dealii.org/current/readme.html"));
1120#define InstantiateUMFPACK(MatrixType) \
1121 template void SparseDirectUMFPACK::factorize(const MatrixType &); \
1122 template void SparseDirectUMFPACK::solve(const MatrixType &, \
1125 template void SparseDirectUMFPACK::solve(const MatrixType &, \
1126 Vector<std::complex<double>> &, \
1128 template void SparseDirectUMFPACK::solve(const MatrixType &, \
1129 BlockVector<double> &, \
1131 template void SparseDirectUMFPACK::solve( \
1132 const MatrixType &, BlockVector<std::complex<double>> &, const bool); \
1133 template void SparseDirectUMFPACK::initialize(const MatrixType &, \
1134 const AdditionalData)
1145#ifdef DEAL_II_WITH_COMPLEX_VALUES
1153#ifdef DEAL_II_WITH_MUMPS
1154# define InstantiateMUMPS(MATRIX) \
1155 template void SparseDirectMUMPS::initialize(const MATRIX &);
virtual size_type size() const override
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
SparseDirectMUMPS(const AdditionalData &additional_data=AdditionalData(), const MPI_Comm &communicator=MPI_COMM_WORLD)
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::global_dof_index size_type
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
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 & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
@ 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)