29#ifdef DEAL_II_WITH_UMFPACK
45 template <
typename SparseMatrixType>
47 parallel_grainsize(
const SparseMatrixType &matrix)
49 const unsigned int avg_entries_per_row =
51 return std::max(1000 / avg_entries_per_row, 1u);
67#ifdef DEAL_II_WITH_UMFPACK
72 , symbolic_decomposition(nullptr)
73 , numeric_decomposition(nullptr)
74 , control(UMFPACK_CONTROL)
76 umfpack_dl_defaults(
control.data());
98 std::vector<types::suitesparse_index> tmp;
103 std::vector<types::suitesparse_index> tmp;
108 std::vector<double> tmp;
113 std::vector<double> tmp;
117 umfpack_dl_defaults(
control.data());
122template <
typename number>
137 for (size_type row = row_begin; row < row_end; ++row)
150 long int cursor = Ap[row];
151 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
153 std::swap(Ai[cursor], Ai[cursor + 1]);
155 std::swap(Ax[cursor], Ax[cursor + 1]);
156 if (numbers::NumberTraits<number>::is_complex == true)
157 std::swap(Az[cursor], Az[cursor + 1]);
163 parallel_grainsize(matrix));
168template <
typename number>
177 for (size_type row = row_begin; row < row_end; ++row)
179 long int cursor = Ap[row];
180 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
182 std::swap(Ai[cursor], Ai[cursor + 1]);
184 std::swap(Ax[cursor], Ax[cursor + 1]);
185 if (numbers::NumberTraits<number>::is_complex == true)
186 std::swap(Az[cursor], Az[cursor + 1]);
192 parallel_grainsize(matrix));
197template <
typename number>
210 for (size_type row = row_begin; row < row_end; ++row)
212 long int cursor = Ap[row];
213 for (size_type block = 0; block < matrix.n_block_cols(); ++block)
216 while ((cursor < Ap[row + 1] - 1) &&
217 (Ai[cursor] < Ai[cursor + 1]))
221 if (cursor == Ap[row + 1] - 1)
226 long int element = cursor;
227 while ((element < Ap[row + 1] - 1) &&
228 (Ai[element] > Ai[element + 1]))
230 std::swap(Ai[element], Ai[element + 1]);
232 std::swap(Ax[element], Ax[element + 1]);
233 if (numbers::NumberTraits<number>::is_complex == true)
234 std::swap(Az[element], Az[element + 1]);
241 parallel_grainsize(matrix));
246template <
class Matrix>
254 using number =
typename Matrix::value_type;
280 Ai.resize(matrix.n_nonzero_elements());
281 Ax.resize(matrix.n_nonzero_elements());
283 Az.resize(matrix.n_nonzero_elements());
288 Ap[row] =
Ap[row - 1] + matrix.get_row_length(row - 1);
298 for (size_type row = row_begin; row < row_end; ++row)
300 long int index = Ap[row];
301 for (typename Matrix::const_iterator p = matrix.begin(row);
302 p != matrix.end(row);
306 Ai[index] = p->column();
307 Ax[index] = std::real(p->value());
308 if (numbers::NumberTraits<number>::is_complex == true)
309 Az[index] = std::imag(p->value());
314 Assert(index == Ap[row + 1], ExcInternalError());
317 parallel_grainsize(matrix));
326 status = umfpack_dl_symbolic(N,
331 &symbolic_decomposition,
335 status = umfpack_zl_symbolic(N,
341 &symbolic_decomposition,
345 ExcUMFPACKError(
"umfpack_dl_symbolic", status));
348 status = umfpack_dl_numeric(Ap.data(),
351 symbolic_decomposition,
352 &numeric_decomposition,
356 status = umfpack_zl_numeric(Ap.data(),
360 symbolic_decomposition,
361 &numeric_decomposition,
365 ExcUMFPACKError(
"umfpack_dl_numeric", status));
367 umfpack_dl_free_symbolic(&symbolic_decomposition);
382 ExcMessage(
"You have previously factored a matrix using this class "
383 "that had complex-valued entries. This then requires "
384 "applying the factored matrix to a complex-valued "
385 "vector, but you are only providing a real-valued vector "
389 rhs = rhs_and_solution;
397 const int status = umfpack_dl_solve(
transpose ? UMFPACK_A : UMFPACK_At,
401 rhs_and_solution.
begin(),
416# ifdef DEAL_II_WITH_COMPLEX_VALUES
446 for (
unsigned int i = 0; i < rhs_and_solution.size(); ++i)
448 rhs_re(i) = std::real(rhs_and_solution(i));
449 rhs_im(i) = std::imag(rhs_and_solution(i));
464 const int status = umfpack_zl_solve(
transpose ? UMFPACK_A : UMFPACK_Aat,
480 for (
unsigned int i = 0; i < rhs_and_solution.size(); ++i)
481 rhs_and_solution(i) = {solution_re(i), solution_im(i)};
496 for (
unsigned int i = 0; i < rhs.
size(); ++i)
497 rhs_real_or_imag(i) = std::real(rhs(i));
501 rhs_and_solution = rhs_real_or_imag;
506 for (
unsigned int i = 0; i < rhs.
size(); ++i)
507 rhs_real_or_imag(i) = std::imag(rhs(i));
511 for (
unsigned int i = 0; i < rhs.
size(); ++i)
512 rhs_and_solution(i).imag(rhs_real_or_imag(i));
517 (void)rhs_and_solution;
521 "This function can't be called if deal.II has been configured "
522 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
535 tmp = rhs_and_solution;
537 rhs_and_solution = tmp;
546# ifdef DEAL_II_WITH_COMPLEX_VALUES
551 tmp = rhs_and_solution;
553 rhs_and_solution = tmp;
556 (void)rhs_and_solution;
560 "This function can't be called if deal.II has been configured "
561 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
567template <
class Matrix>
579template <
class Matrix>
582 Vector<std::complex<double>> &rhs_and_solution,
585# ifdef DEAL_II_WITH_COMPLEX_VALUES
592 (void)rhs_and_solution;
596 "This function can't be called if deal.II has been configured "
597 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
603template <
class Matrix>
615template <
class Matrix>
618 BlockVector<std::complex<double>> &rhs_and_solution,
621# ifdef DEAL_II_WITH_COMPLEX_VALUES
628 (void)rhs_and_solution;
632 "This function can't be called if deal.II has been configured "
633 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
644 , symbolic_decomposition(nullptr)
645 , numeric_decomposition(nullptr)
655template <
class Matrix>
662 "To call this function you need UMFPACK, but you configured deal.II "
663 "without passing the necessary switch to 'cmake'. Please consult the "
664 "installation instructions at https://dealii.org/current/readme.html"));
674 "To call this function you need UMFPACK, but you configured deal.II "
675 "without passing the necessary switch to 'cmake'. Please consult the "
676 "installation instructions at https://dealii.org/current/readme.html"));
687 "To call this function you need UMFPACK, but you configured deal.II "
688 "without passing the necessary switch to 'cmake'. Please consult the "
689 "installation instructions at https://dealii.org/current/readme.html"));
700 "To call this function you need UMFPACK, but you configured deal.II "
701 "without passing the necessary switch to 'cmake'. Please consult the "
702 "installation instructions at https://dealii.org/current/readme.html"));
714 "To call this function you need UMFPACK, but you configured deal.II "
715 "without passing the necessary switch to 'cmake'. Please consult the "
716 "installation instructions at https://dealii.org/current/readme.html"));
721template <
class Matrix>
728 "To call this function you need UMFPACK, but you configured deal.II "
729 "without passing the necessary switch to 'cmake'. Please consult the "
730 "installation instructions at https://dealii.org/current/readme.html"));
735template <
class Matrix>
738 Vector<std::complex<double>> &,
744 "To call this function you need UMFPACK, but you configured deal.II "
745 "without passing the necessary switch to 'cmake'. Please consult the "
746 "installation instructions at https://dealii.org/current/readme.html"));
751template <
class Matrix>
758 "To call this function you need UMFPACK, but you configured deal.II "
759 "without passing the necessary switch to 'cmake'. Please consult the "
760 "installation instructions at https://dealii.org/current/readme.html"));
765template <
class Matrix>
774 "To call this function you need UMFPACK, but you configured deal.II "
775 "without passing the necessary switch to 'cmake'. Please consult the "
776 "installation instructions at https://dealii.org/current/readme.html"));
782template <
class Matrix>
813 this->
solve(dst,
true);
823 this->
solve(dst,
true);
842#define InstantiateUMFPACK(MatrixType) \
843 template void SparseDirectUMFPACK::factorize(const MatrixType &); \
844 template void SparseDirectUMFPACK::solve(const MatrixType &, \
847 template void SparseDirectUMFPACK::solve(const MatrixType &, \
848 Vector<std::complex<double>> &, \
850 template void SparseDirectUMFPACK::solve(const MatrixType &, \
851 BlockVector<double> &, \
853 template void SparseDirectUMFPACK::solve( \
854 const MatrixType &, BlockVector<std::complex<double>> &, const bool); \
855 template void SparseDirectUMFPACK::initialize(const MatrixType &, \
856 const AdditionalData)
842#define InstantiateUMFPACK(MatrixType) \ …
867#ifdef DEAL_II_WITH_COMPLEX_VALUES
virtual size_type size() const override
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)
@ matrix
Contents is actually a matrix.
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 > &)
#define InstantiateUMFPACK(MatrixType)