30#ifdef DEAL_II_WITH_UMFPACK
46 template <
typename SparseMatrixType>
48 parallel_grainsize(
const SparseMatrixType &matrix)
50 const unsigned int avg_entries_per_row =
52 return std::max(1000 / avg_entries_per_row, 1u);
68#ifdef DEAL_II_WITH_UMFPACK
73 , symbolic_decomposition(nullptr)
74 , numeric_decomposition(nullptr)
75 , control(UMFPACK_CONTROL)
77 umfpack_dl_defaults(
control.data());
99 std::vector<types::suitesparse_index> tmp;
104 std::vector<types::suitesparse_index> tmp;
109 std::vector<double> tmp;
114 std::vector<double> tmp;
118 umfpack_dl_defaults(
control.data());
123template <
typename number>
138 for (size_type row = row_begin; row < row_end; ++row)
151 long int cursor = Ap[row];
152 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
154 std::swap(Ai[cursor], Ai[cursor + 1]);
156 std::swap(Ax[cursor], Ax[cursor + 1]);
157 if (numbers::NumberTraits<number>::is_complex == true)
158 std::swap(Az[cursor], Az[cursor + 1]);
164 parallel_grainsize(matrix));
169template <
typename number>
178 for (size_type row = row_begin; row < row_end; ++row)
180 long int cursor = Ap[row];
181 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
183 std::swap(Ai[cursor], Ai[cursor + 1]);
185 std::swap(Ax[cursor], Ax[cursor + 1]);
186 if (numbers::NumberTraits<number>::is_complex == true)
187 std::swap(Az[cursor], Az[cursor + 1]);
193 parallel_grainsize(matrix));
198template <
typename number>
211 for (size_type row = row_begin; row < row_end; ++row)
213 long int cursor = Ap[row];
214 for (size_type block = 0; block < matrix.n_block_cols(); ++block)
217 while ((cursor < Ap[row + 1] - 1) &&
218 (Ai[cursor] < Ai[cursor + 1]))
222 if (cursor == Ap[row + 1] - 1)
227 long int element = cursor;
228 while ((element < Ap[row + 1] - 1) &&
229 (Ai[element] > Ai[element + 1]))
231 std::swap(Ai[element], Ai[element + 1]);
233 std::swap(Ax[element], Ax[element + 1]);
234 if (numbers::NumberTraits<number>::is_complex == true)
235 std::swap(Az[element], Az[element + 1]);
242 parallel_grainsize(matrix));
247template <
class Matrix>
255 using number =
typename Matrix::value_type;
281 Ai.resize(matrix.n_nonzero_elements());
282 Ax.resize(matrix.n_nonzero_elements());
284 Az.resize(matrix.n_nonzero_elements());
289 Ap[row] =
Ap[row - 1] + matrix.get_row_length(row - 1);
299 for (size_type row = row_begin; row < row_end; ++row)
301 long int index = Ap[row];
302 for (typename Matrix::const_iterator p = matrix.begin(row);
303 p != matrix.end(row);
307 Ai[index] = p->column();
308 Ax[index] = std::real(p->value());
309 if (numbers::NumberTraits<number>::is_complex == true)
310 Az[index] = std::imag(p->value());
315 Assert(index == Ap[row + 1], ExcInternalError());
318 parallel_grainsize(matrix));
327 status = umfpack_dl_symbolic(N,
332 &symbolic_decomposition,
336 status = umfpack_zl_symbolic(N,
342 &symbolic_decomposition,
346 ExcUMFPACKError(
"umfpack_dl_symbolic", status));
349 status = umfpack_dl_numeric(Ap.data(),
352 symbolic_decomposition,
353 &numeric_decomposition,
357 status = umfpack_zl_numeric(Ap.data(),
361 symbolic_decomposition,
362 &numeric_decomposition,
366 ExcUMFPACKError(
"umfpack_dl_numeric", status));
368 umfpack_dl_free_symbolic(&symbolic_decomposition);
383 ExcMessage(
"You have previously factored a matrix using this class "
384 "that had complex-valued entries. This then requires "
385 "applying the factored matrix to a complex-valued "
386 "vector, but you are only providing a real-valued vector "
390 rhs = rhs_and_solution;
398 const int status = umfpack_dl_solve(
transpose ? UMFPACK_A : UMFPACK_At,
402 rhs_and_solution.
begin(),
417# ifdef DEAL_II_WITH_COMPLEX_VALUES
447 for (
unsigned int i = 0; i < rhs_and_solution.size(); ++i)
449 rhs_re(i) = std::real(rhs_and_solution(i));
450 rhs_im(i) = std::imag(rhs_and_solution(i));
465 const int status = umfpack_zl_solve(
transpose ? UMFPACK_A : UMFPACK_Aat,
481 for (
unsigned int i = 0; i < rhs_and_solution.size(); ++i)
482 rhs_and_solution(i) = {solution_re(i), solution_im(i)};
497 for (
unsigned int i = 0; i < rhs.
size(); ++i)
498 rhs_real_or_imag(i) = std::real(rhs(i));
502 rhs_and_solution = rhs_real_or_imag;
507 for (
unsigned int i = 0; i < rhs.
size(); ++i)
508 rhs_real_or_imag(i) = std::imag(rhs(i));
512 for (
unsigned int i = 0; i < rhs.
size(); ++i)
513 rhs_and_solution(i).imag(rhs_real_or_imag(i));
518 (void)rhs_and_solution;
522 "This function can't be called if deal.II has been configured "
523 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
536 tmp = rhs_and_solution;
538 rhs_and_solution = tmp;
547# ifdef DEAL_II_WITH_COMPLEX_VALUES
552 tmp = rhs_and_solution;
554 rhs_and_solution = tmp;
557 (void)rhs_and_solution;
561 "This function can't be called if deal.II has been configured "
562 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
568template <
class Matrix>
580template <
class Matrix>
583 Vector<std::complex<double>> &rhs_and_solution,
586# ifdef DEAL_II_WITH_COMPLEX_VALUES
593 (void)rhs_and_solution;
597 "This function can't be called if deal.II has been configured "
598 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
604template <
class Matrix>
616template <
class Matrix>
619 BlockVector<std::complex<double>> &rhs_and_solution,
622# ifdef DEAL_II_WITH_COMPLEX_VALUES
629 (void)rhs_and_solution;
633 "This function can't be called if deal.II has been configured "
634 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
645 , symbolic_decomposition(nullptr)
646 , numeric_decomposition(nullptr)
656template <
class Matrix>
663 "To call this function you need UMFPACK, but you configured deal.II "
664 "without passing the necessary switch to 'cmake'. Please consult the "
665 "installation instructions at https://dealii.org/current/readme.html"));
675 "To call this function you need UMFPACK, but you configured deal.II "
676 "without passing the necessary switch to 'cmake'. Please consult the "
677 "installation instructions at https://dealii.org/current/readme.html"));
688 "To call this function you need UMFPACK, but you configured deal.II "
689 "without passing the necessary switch to 'cmake'. Please consult the "
690 "installation instructions at https://dealii.org/current/readme.html"));
701 "To call this function you need UMFPACK, but you configured deal.II "
702 "without passing the necessary switch to 'cmake'. Please consult the "
703 "installation instructions at https://dealii.org/current/readme.html"));
715 "To call this function you need UMFPACK, but you configured deal.II "
716 "without passing the necessary switch to 'cmake'. Please consult the "
717 "installation instructions at https://dealii.org/current/readme.html"));
722template <
class Matrix>
729 "To call this function you need UMFPACK, but you configured deal.II "
730 "without passing the necessary switch to 'cmake'. Please consult the "
731 "installation instructions at https://dealii.org/current/readme.html"));
736template <
class Matrix>
739 Vector<std::complex<double>> &,
745 "To call this function you need UMFPACK, but you configured deal.II "
746 "without passing the necessary switch to 'cmake'. Please consult the "
747 "installation instructions at https://dealii.org/current/readme.html"));
752template <
class Matrix>
759 "To call this function you need UMFPACK, but you configured deal.II "
760 "without passing the necessary switch to 'cmake'. Please consult the "
761 "installation instructions at https://dealii.org/current/readme.html"));
766template <
class Matrix>
775 "To call this function you need UMFPACK, but you configured deal.II "
776 "without passing the necessary switch to 'cmake'. Please consult the "
777 "installation instructions at https://dealii.org/current/readme.html"));
783template <
class Matrix>
814 this->
solve(dst,
true);
824 this->
solve(dst,
true);
843#define InstantiateUMFPACK(MatrixType) \
844 template void SparseDirectUMFPACK::factorize(const MatrixType &); \
845 template void SparseDirectUMFPACK::solve(const MatrixType &, \
848 template void SparseDirectUMFPACK::solve(const MatrixType &, \
849 Vector<std::complex<double>> &, \
851 template void SparseDirectUMFPACK::solve(const MatrixType &, \
852 BlockVector<double> &, \
854 template void SparseDirectUMFPACK::solve( \
855 const MatrixType &, BlockVector<std::complex<double>> &, const bool); \
856 template void SparseDirectUMFPACK::initialize(const MatrixType &, \
857 const AdditionalData)
843#define InstantiateUMFPACK(MatrixType) \ …
868#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)