Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2850-g1646a780ac 2025-03-17 15:10:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
sparse_direct.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
21#include <deal.II/lac/vector.h>
22
23#include <complex>
24#include <iostream>
25#include <list>
26#include <typeinfo>
27#include <vector>
28
29#ifdef DEAL_II_WITH_UMFPACK
30# include <umfpack.h>
31#endif
32
33
35
36namespace
37{
45 template <typename SparseMatrixType>
46 unsigned int
47 parallel_grainsize(const SparseMatrixType &matrix)
48 {
49 const unsigned int avg_entries_per_row =
50 matrix.n_nonzero_elements() / matrix.m();
51 return std::max(1000 / avg_entries_per_row, 1u);
52 }
53} // namespace
54
55
60
61
62void
65
66
67#ifdef DEAL_II_WITH_UMFPACK
68
70 : n_rows(0)
71 , n_cols(0)
72 , symbolic_decomposition(nullptr)
73 , numeric_decomposition(nullptr)
74 , control(UMFPACK_CONTROL)
75{
76 umfpack_dl_defaults(control.data());
77}
78
79
80
81void
83{
84 // delete objects that haven't been deleted yet
85 if (symbolic_decomposition != nullptr)
86 {
87 umfpack_dl_free_symbolic(&symbolic_decomposition);
88 symbolic_decomposition = nullptr;
89 }
90
91 if (numeric_decomposition != nullptr)
92 {
93 umfpack_dl_free_numeric(&numeric_decomposition);
94 numeric_decomposition = nullptr;
95 }
96
97 {
98 std::vector<types::suitesparse_index> tmp;
99 tmp.swap(Ap);
100 }
101
102 {
103 std::vector<types::suitesparse_index> tmp;
104 tmp.swap(Ai);
105 }
106
107 {
108 std::vector<double> tmp;
109 tmp.swap(Ax);
110 }
111
112 {
113 std::vector<double> tmp;
114 tmp.swap(Az);
115 }
116
117 umfpack_dl_defaults(control.data());
118}
119
120
121
122template <typename number>
123void
125{
126 // do the copying around of entries so that the diagonal entry is in the
127 // right place. note that this is easy to detect: since all entries apart
128 // from the diagonal entry are sorted, we know that the diagonal entry is
129 // in the wrong place if and only if its column index is larger than the
130 // column index of the second entry in a row
131 //
132 // ignore rows with only one or no entry
134 0,
135 matrix.m(),
136 [this](const size_type row_begin, const size_type row_end) {
137 for (size_type row = row_begin; row < row_end; ++row)
138 {
139 // we may have to move some elements that are left of the diagonal
140 // but presently after the diagonal entry to the left, whereas the
141 // diagonal entry has to move to the right. we could first figure out
142 // where to move everything to, but for simplicity we just make a
143 // series of swaps instead (this is kind of a single run of
144 // bubble-sort, which gives us the desired result since the array is
145 // already "almost" sorted)
146 //
147 // in the first loop, the condition in the while-header also checks
148 // that the row has at least two entries and that the diagonal entry
149 // is really in the wrong place
150 long int cursor = Ap[row];
151 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
152 {
153 std::swap(Ai[cursor], Ai[cursor + 1]);
154
155 std::swap(Ax[cursor], Ax[cursor + 1]);
156 if (numbers::NumberTraits<number>::is_complex == true)
157 std::swap(Az[cursor], Az[cursor + 1]);
158
159 ++cursor;
160 }
161 }
162 },
163 parallel_grainsize(matrix));
164}
165
166
167
168template <typename number>
169void
171{
172 // same thing for SparseMatrixEZ
174 0,
175 matrix.m(),
176 [this](const size_type row_begin, const size_type row_end) {
177 for (size_type row = row_begin; row < row_end; ++row)
178 {
179 long int cursor = Ap[row];
180 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
181 {
182 std::swap(Ai[cursor], Ai[cursor + 1]);
183
184 std::swap(Ax[cursor], Ax[cursor + 1]);
185 if (numbers::NumberTraits<number>::is_complex == true)
186 std::swap(Az[cursor], Az[cursor + 1]);
187
188 ++cursor;
189 }
190 }
191 },
192 parallel_grainsize(matrix));
193}
194
195
196
197template <typename number>
198void
200{
201 // the case for block matrices is a bit more difficult, since all we know
202 // is that *within each block*, the diagonal of that block may come
203 // first. however, that means that there may be as many entries per row
204 // in the wrong place as there are block columns. we can do the same
205 // thing as above, but we have to do it multiple times
207 0,
208 matrix.m(),
209 [this, &matrix](const size_type row_begin, const size_type row_end) {
210 for (size_type row = row_begin; row < row_end; ++row)
211 {
212 long int cursor = Ap[row];
213 for (size_type block = 0; block < matrix.n_block_cols(); ++block)
214 {
215 // find the next out-of-order element
216 while ((cursor < Ap[row + 1] - 1) &&
217 (Ai[cursor] < Ai[cursor + 1]))
218 ++cursor;
219
220 // if there is none, then just go on
221 if (cursor == Ap[row + 1] - 1)
222 break;
223
224 // otherwise swap this entry with successive ones as long as
225 // necessary
226 long int element = cursor;
227 while ((element < Ap[row + 1] - 1) &&
228 (Ai[element] > Ai[element + 1]))
229 {
230 std::swap(Ai[element], Ai[element + 1]);
231
232 std::swap(Ax[element], Ax[element + 1]);
233 if (numbers::NumberTraits<number>::is_complex == true)
234 std::swap(Az[element], Az[element + 1]);
235
236 ++element;
237 }
238 }
239 }
240 },
241 parallel_grainsize(matrix));
242}
243
244
245
246template <class Matrix>
247void
249{
250 Assert(matrix.m() == matrix.n(), ExcNotQuadratic());
251
252 clear();
253
254 using number = typename Matrix::value_type;
255
256 n_rows = matrix.m();
257 n_cols = matrix.n();
258
259 const size_type N = matrix.m();
260
261 // copy over the data from the matrix to the data structures UMFPACK
262 // wants. note two things: first, UMFPACK wants compressed column storage
263 // whereas we always do compressed row storage; we work around this by,
264 // rather than shuffling things around, copy over the data we have, but
265 // then call the umfpack_dl_solve function with the UMFPACK_At argument,
266 // meaning that we want to solve for the transpose system
267 //
268 // second: the data we have in the sparse matrices is "almost" right
269 // already; UMFPACK wants the entries in each row (i.e. really: column)
270 // to be sorted in ascending order. we almost have that, except that we
271 // usually store the diagonal first in each row to allow for some
272 // optimizations. thus, we have to resort things a little bit, but only
273 // within each row
274 //
275 // final note: if the matrix has entries in the sparsity pattern that are
276 // actually occupied by entries that have a zero numerical value, then we
277 // keep them anyway. people are supposed to provide accurate sparsity
278 // patterns.
279 Ap.resize(N + 1);
280 Ai.resize(matrix.n_nonzero_elements());
281 Ax.resize(matrix.n_nonzero_elements());
283 Az.resize(matrix.n_nonzero_elements());
284
285 // first fill row lengths array
286 Ap[0] = 0;
287 for (size_type row = 1; row <= N; ++row)
288 Ap[row] = Ap[row - 1] + matrix.get_row_length(row - 1);
289 Assert(static_cast<size_type>(Ap.back()) == Ai.size(), ExcInternalError());
290
291 // then copy over matrix elements. note that for sparse matrices,
292 // iterators are sorted so that they traverse each row from start to end
293 // before moving on to the next row.
295 0,
296 matrix.m(),
297 [this, &matrix](const size_type row_begin, const size_type row_end) {
298 for (size_type row = row_begin; row < row_end; ++row)
299 {
300 long int index = Ap[row];
301 for (typename Matrix::const_iterator p = matrix.begin(row);
302 p != matrix.end(row);
303 ++p)
304 {
305 // write entry into the first free one for this 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());
310
311 // then move pointer ahead
312 ++index;
313 }
314 Assert(index == Ap[row + 1], ExcInternalError());
315 }
316 },
317 parallel_grainsize(matrix));
318
319 // make sure that the elements in each row are sorted. we have to be more
320 // careful for block sparse matrices, so ship this task out to a
321 // different function
322 sort_arrays(matrix);
323
324 int status;
326 status = umfpack_dl_symbolic(N,
327 N,
328 Ap.data(),
329 Ai.data(),
330 Ax.data(),
331 &symbolic_decomposition,
332 control.data(),
333 nullptr);
334 else
335 status = umfpack_zl_symbolic(N,
336 N,
337 Ap.data(),
338 Ai.data(),
339 Ax.data(),
340 Az.data(),
341 &symbolic_decomposition,
342 control.data(),
343 nullptr);
344 AssertThrow(status == UMFPACK_OK,
345 ExcUMFPACKError("umfpack_dl_symbolic", status));
346
348 status = umfpack_dl_numeric(Ap.data(),
349 Ai.data(),
350 Ax.data(),
351 symbolic_decomposition,
352 &numeric_decomposition,
353 control.data(),
354 nullptr);
355 else
356 status = umfpack_zl_numeric(Ap.data(),
357 Ai.data(),
358 Ax.data(),
359 Az.data(),
360 symbolic_decomposition,
361 &numeric_decomposition,
362 control.data(),
363 nullptr);
364 AssertThrow(status == UMFPACK_OK,
365 ExcUMFPACKError("umfpack_dl_numeric", status));
366
367 umfpack_dl_free_symbolic(&symbolic_decomposition);
368}
369
370
371
372void
374 const bool transpose /*=false*/) const
375{
376 // make sure that some kind of factorize() call has happened before
377 Assert(Ap.size() != 0, ExcNotInitialized());
378 Assert(Ai.size() != 0, ExcNotInitialized());
379 Assert(Ai.size() == Ax.size(), ExcNotInitialized());
380
381 Assert(Az.empty(),
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 "
386 "here."));
387
388 Vector<double> rhs(rhs_and_solution.size());
389 rhs = rhs_and_solution;
390
391 // solve the system. note that since UMFPACK wants compressed column
392 // storage instead of the compressed row storage format we use in
393 // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
394
395 // Conversely, if we solve for the transpose, we have to use UMFPACK_A
396 // instead.
397 const int status = umfpack_dl_solve(transpose ? UMFPACK_A : UMFPACK_At,
398 Ap.data(),
399 Ai.data(),
400 Ax.data(),
401 rhs_and_solution.begin(),
402 rhs.begin(),
404 control.data(),
405 nullptr);
406 AssertThrow(status == UMFPACK_OK,
407 ExcUMFPACKError("umfpack_dl_solve", status));
408}
409
410
411
412void
413SparseDirectUMFPACK::solve(Vector<std::complex<double>> &rhs_and_solution,
414 const bool transpose /*=false*/) const
415{
416# ifdef DEAL_II_WITH_COMPLEX_VALUES
417 // make sure that some kind of factorize() call has happened before
418 Assert(Ap.size() != 0, ExcNotInitialized());
419 Assert(Ai.size() != 0, ExcNotInitialized());
420 Assert(Ai.size() == Ax.size(), ExcNotInitialized());
421
422 // First see whether the matrix that was factorized was complex-valued.
423 // If so, just apply the complex factorization to the vector.
424 if (Az.size() != 0)
425 {
426 Assert(Ax.size() == Az.size(), ExcInternalError());
427
428 // It would be nice if we could just present a pointer to the
429 // first element of the complex-valued solution vector and let
430 // UMFPACK fill both the real and imaginary parts of the solution
431 // vector at that address. UMFPACK calls this the 'packed' format,
432 // and in those cases it only takes one pointer to the entire
433 // vector, rather than a pointer to the real and one pointer to
434 // the imaginary parts of the vector. The problem is that if we
435 // want to pack, then we also need to pack the matrix, and the
436 // functions above have already decided that we don't want to pack
437 // the matrix but instead deal with split format for the matrix,
438 // and then UMFPACK decides that it can't deal with a split matrix
439 // and a packed vector. We have to choose one or the other, not
440 // mix.
441 //
442 // So create four vectors, one each for the real and imaginary parts
443 // of the right hand side and solution.
444 Vector<double> rhs_re(rhs_and_solution.size());
445 Vector<double> rhs_im(rhs_and_solution.size());
446 for (unsigned int i = 0; i < rhs_and_solution.size(); ++i)
447 {
448 rhs_re(i) = std::real(rhs_and_solution(i));
449 rhs_im(i) = std::imag(rhs_and_solution(i));
450 }
451
452 Vector<double> solution_re(rhs_and_solution.size());
453 Vector<double> solution_im(rhs_and_solution.size());
454
455 // Solve the system. note that since UMFPACK wants compressed column
456 // storage instead of the compressed row storage format we use in
457 // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
458 //
459 // Conversely, if we solve for the transpose, we have to use UMFPACK_A
460 // instead.
461 //
462 // Note that for the complex case here, the transpose is selected using
463 // UMFPACK_Aat, not UMFPACK_At.
464 const int status = umfpack_zl_solve(transpose ? UMFPACK_A : UMFPACK_Aat,
465 Ap.data(),
466 Ai.data(),
467 Ax.data(),
468 Az.data(),
469 solution_re.data(),
470 solution_im.data(),
471 rhs_re.data(),
472 rhs_im.data(),
474 control.data(),
475 nullptr);
476 AssertThrow(status == UMFPACK_OK,
477 ExcUMFPACKError("umfpack_dl_solve", status));
478
479 // Now put things back together into the output vector
480 for (unsigned int i = 0; i < rhs_and_solution.size(); ++i)
481 rhs_and_solution(i) = {solution_re(i), solution_im(i)};
482 }
483 else
484 {
485 // We have factorized a real-valued matrix, but the rhs and solution
486 // vectors are complex-valued. UMFPACK does not natively support this
487 // case, but we can just apply the factorization to real and imaginary
488 // parts of the right hand side separately
489 const Vector<std::complex<double>> rhs = rhs_and_solution;
490
491 // Get the real part of the right hand side, solve with it, and copy the
492 // results into the result vector by just copying the real output
493 // into the complex-valued result vector (which implies setting the
494 // imaginary parts to zero):
495 Vector<double> rhs_real_or_imag(rhs_and_solution.size());
496 for (unsigned int i = 0; i < rhs.size(); ++i)
497 rhs_real_or_imag(i) = std::real(rhs(i));
498
499 solve(rhs_real_or_imag, transpose);
500
501 rhs_and_solution = rhs_real_or_imag;
502
503 // Then repeat the whole thing with the imaginary part. The copying step
504 // is more complicated now because we can only touch the imaginary
505 // component of the output vector:
506 for (unsigned int i = 0; i < rhs.size(); ++i)
507 rhs_real_or_imag(i) = std::imag(rhs(i));
508
509 solve(rhs_real_or_imag, transpose);
510
511 for (unsigned int i = 0; i < rhs.size(); ++i)
512 rhs_and_solution(i).imag(rhs_real_or_imag(i));
513 }
514
515# else
516
517 (void)rhs_and_solution;
518 (void)transpose;
519 Assert(false,
521 "This function can't be called if deal.II has been configured "
522 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
523# endif
524}
525
526
527void
529 const bool transpose /*=false*/) const
530{
531 // the UMFPACK functions want a contiguous array of elements, so
532 // there is no way around copying data around. thus, just copy the
533 // data into a regular vector and back
534 Vector<double> tmp(rhs_and_solution.size());
535 tmp = rhs_and_solution;
536 solve(tmp, transpose);
537 rhs_and_solution = tmp;
538}
539
540
541
542void
543SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &rhs_and_solution,
544 const bool transpose /*=false*/) const
545{
546# ifdef DEAL_II_WITH_COMPLEX_VALUES
547 // the UMFPACK functions want a contiguous array of elements, so
548 // there is no way around copying data around. thus, just copy the
549 // data into a regular vector and back
550 Vector<std::complex<double>> tmp(rhs_and_solution.size());
551 tmp = rhs_and_solution;
552 solve(tmp, transpose);
553 rhs_and_solution = tmp;
554
555# else
556 (void)rhs_and_solution;
557 (void)transpose;
558 Assert(false,
560 "This function can't be called if deal.II has been configured "
561 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
562# endif
563}
564
565
566
567template <class Matrix>
568void
569SparseDirectUMFPACK::solve(const Matrix &matrix,
570 Vector<double> &rhs_and_solution,
571 const bool transpose /*=false*/)
572{
573 factorize(matrix);
574 solve(rhs_and_solution, transpose);
575}
576
577
578
579template <class Matrix>
580void
581SparseDirectUMFPACK::solve(const Matrix &matrix,
582 Vector<std::complex<double>> &rhs_and_solution,
583 const bool transpose /*=false*/)
584{
585# ifdef DEAL_II_WITH_COMPLEX_VALUES
586 factorize(matrix);
587 solve(rhs_and_solution, transpose);
588
589# else
590
591 (void)matrix;
592 (void)rhs_and_solution;
593 (void)transpose;
594 Assert(false,
596 "This function can't be called if deal.II has been configured "
597 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
598# endif
599}
600
601
602
603template <class Matrix>
604void
605SparseDirectUMFPACK::solve(const Matrix &matrix,
606 BlockVector<double> &rhs_and_solution,
607 const bool transpose /*=false*/)
608{
609 factorize(matrix);
610 solve(rhs_and_solution, transpose);
611}
612
613
614
615template <class Matrix>
616void
617SparseDirectUMFPACK::solve(const Matrix &matrix,
618 BlockVector<std::complex<double>> &rhs_and_solution,
619 const bool transpose /*=false*/)
620{
621# ifdef DEAL_II_WITH_COMPLEX_VALUES
622 factorize(matrix);
623 solve(rhs_and_solution, transpose);
624
625# else
626
627 (void)matrix;
628 (void)rhs_and_solution;
629 (void)transpose;
630 Assert(false,
632 "This function can't be called if deal.II has been configured "
633 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
634# endif
635}
636
637
638#else
639
640
642 : n_rows(0)
643 , n_cols(0)
644 , symbolic_decomposition(nullptr)
645 , numeric_decomposition(nullptr)
646 , control(0)
647{}
648
649
650void
652{}
653
654
655template <class Matrix>
656void
657SparseDirectUMFPACK::factorize(const Matrix &)
658{
660 false,
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"));
665}
666
667
668void
670{
672 false,
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"));
677}
678
679
680
681void
682SparseDirectUMFPACK::solve(Vector<std::complex<double>> &, const bool) const
683{
685 false,
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"));
690}
691
692
693
694void
696{
698 false,
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"));
703}
704
705
706
707void
708SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &,
709 const bool) const
710{
712 false,
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"));
717}
718
719
720
721template <class Matrix>
722void
723SparseDirectUMFPACK::solve(const Matrix &, Vector<double> &, const bool)
724{
726 false,
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"));
731}
732
733
734
735template <class Matrix>
736void
737SparseDirectUMFPACK::solve(const Matrix &,
738 Vector<std::complex<double>> &,
739 const bool)
740{
742 false,
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"));
747}
748
749
750
751template <class Matrix>
752void
753SparseDirectUMFPACK::solve(const Matrix &, BlockVector<double> &, const bool)
754{
756 false,
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"));
761}
762
763
764
765template <class Matrix>
766void
767SparseDirectUMFPACK::solve(const Matrix &,
768 BlockVector<std::complex<double>> &,
769 const bool)
770{
772 false,
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"));
777}
778
779#endif
780
781
782template <class Matrix>
783void
785{
786 this->factorize(M);
787}
788
789
790void
792{
793 dst = src;
794 this->solve(dst);
795}
796
797
798
799void
801 const BlockVector<double> &src) const
802{
803 dst = src;
804 this->solve(dst);
805}
806
807
808void
810 const Vector<double> &src) const
811{
812 dst = src;
813 this->solve(dst, /*transpose=*/true);
814}
815
816
817
818void
820 const BlockVector<double> &src) const
821{
822 dst = src;
823 this->solve(dst, /*transpose=*/true);
824}
825
828{
830 return n_rows;
831}
832
835{
837 return n_cols;
838}
839
840
841// explicit instantiations for SparseMatrixUMFPACK
842#define InstantiateUMFPACK(MatrixType) \
843 template void SparseDirectUMFPACK::factorize(const MatrixType &); \
844 template void SparseDirectUMFPACK::solve(const MatrixType &, \
845 Vector<double> &, \
846 const bool); \
847 template void SparseDirectUMFPACK::solve(const MatrixType &, \
848 Vector<std::complex<double>> &, \
849 const bool); \
850 template void SparseDirectUMFPACK::solve(const MatrixType &, \
851 BlockVector<double> &, \
852 const bool); \
853 template void SparseDirectUMFPACK::solve( \
854 const MatrixType &, BlockVector<std::complex<double>> &, const bool); \
855 template void SparseDirectUMFPACK::initialize(const MatrixType &, \
856 const AdditionalData)
857
858// Instantiate everything for real-valued matrices
865
866// Now also for complex-valued matrices
867#ifdef DEAL_II_WITH_COMPLEX_VALUES
868InstantiateUMFPACK(SparseMatrix<std::complex<double>>);
869InstantiateUMFPACK(SparseMatrix<std::complex<float>>);
872#endif
873
virtual size_type size() const override
std::vector< double > Az
~SparseDirectUMFPACK() override
void initialize(const SparsityPattern &sparsity_pattern)
void Tvmult(Vector< double > &dst, const Vector< double > &src) const
size_type m() const
void solve(Vector< double > &rhs_and_solution, const bool transpose=false) const
std::vector< double > Ax
void sort_arrays(const SparseMatrixEZ< number > &)
size_type n() const
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
pointer data()
virtual size_type size() const override
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
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)
Definition parallel.h:540
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
#define InstantiateUMFPACK(MatrixType)