Reference documentation for deal.II version GIT relicensing-43-g3fc2fae7c3 2024-03-03 21:10:02+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\}}\)
Loading...
Searching...
No Matches
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 <cerrno>
24#include <complex>
25#include <iostream>
26#include <list>
27#include <typeinfo>
28#include <vector>
29
30
32
33
34// include UMFPACK file.
35#ifdef DEAL_II_WITH_UMFPACK
36# include <umfpack.h>
37#endif
38
39
40namespace
41{
49 template <typename SparseMatrixType>
50 unsigned int
51 parallel_grainsize(const SparseMatrixType &matrix)
52 {
53 const unsigned int avg_entries_per_row =
54 matrix.n_nonzero_elements() / matrix.m();
55 return std::max(1000 / avg_entries_per_row, 1u);
56 }
57} // namespace
58
59
64
65
66void
69
70
71#ifdef DEAL_II_WITH_UMFPACK
72
74 : n_rows(0)
75 , n_cols(0)
76 , symbolic_decomposition(nullptr)
77 , numeric_decomposition(nullptr)
78 , control(UMFPACK_CONTROL)
79{
80 umfpack_dl_defaults(control.data());
81}
82
83
84
85void
87{
88 // delete objects that haven't been deleted yet
89 if (symbolic_decomposition != nullptr)
90 {
91 umfpack_dl_free_symbolic(&symbolic_decomposition);
92 symbolic_decomposition = nullptr;
93 }
94
95 if (numeric_decomposition != nullptr)
96 {
97 umfpack_dl_free_numeric(&numeric_decomposition);
98 numeric_decomposition = nullptr;
99 }
100
101 {
102 std::vector<types::suitesparse_index> tmp;
103 tmp.swap(Ap);
104 }
105
106 {
107 std::vector<types::suitesparse_index> tmp;
108 tmp.swap(Ai);
109 }
110
111 {
112 std::vector<double> tmp;
113 tmp.swap(Ax);
114 }
115
116 {
117 std::vector<double> tmp;
118 tmp.swap(Az);
119 }
120
121 umfpack_dl_defaults(control.data());
122}
123
124
125
126template <typename number>
127void
129{
130 // do the copying around of entries so that the diagonal entry is in the
131 // right place. note that this is easy to detect: since all entries apart
132 // from the diagonal entry are sorted, we know that the diagonal entry is
133 // in the wrong place if and only if its column index is larger than the
134 // column index of the second entry in a row
135 //
136 // ignore rows with only one or no entry
138 0,
139 matrix.m(),
140 [this](const size_type row_begin, const size_type row_end) {
141 for (size_type row = row_begin; row < row_end; ++row)
142 {
143 // we may have to move some elements that are left of the diagonal
144 // but presently after the diagonal entry to the left, whereas the
145 // diagonal entry has to move to the right. we could first figure out
146 // where to move everything to, but for simplicity we just make a
147 // series of swaps instead (this is kind of a single run of
148 // bubble-sort, which gives us the desired result since the array is
149 // already "almost" sorted)
150 //
151 // in the first loop, the condition in the while-header also checks
152 // that the row has at least two entries and that the diagonal entry
153 // is really in the wrong place
154 long int cursor = Ap[row];
155 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
156 {
157 std::swap(Ai[cursor], Ai[cursor + 1]);
158
159 std::swap(Ax[cursor], Ax[cursor + 1]);
160 if (numbers::NumberTraits<number>::is_complex == true)
161 std::swap(Az[cursor], Az[cursor + 1]);
162
163 ++cursor;
164 }
165 }
166 },
167 parallel_grainsize(matrix));
168}
169
170
171
172template <typename number>
173void
175{
176 // same thing for SparseMatrixEZ
178 0,
179 matrix.m(),
180 [this](const size_type row_begin, const size_type row_end) {
181 for (size_type row = row_begin; row < row_end; ++row)
182 {
183 long int cursor = Ap[row];
184 while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
185 {
186 std::swap(Ai[cursor], Ai[cursor + 1]);
187
188 std::swap(Ax[cursor], Ax[cursor + 1]);
189 if (numbers::NumberTraits<number>::is_complex == true)
190 std::swap(Az[cursor], Az[cursor + 1]);
191
192 ++cursor;
193 }
194 }
195 },
196 parallel_grainsize(matrix));
197}
198
199
200
201template <typename number>
202void
204{
205 // the case for block matrices is a bit more difficult, since all we know
206 // is that *within each block*, the diagonal of that block may come
207 // first. however, that means that there may be as many entries per row
208 // in the wrong place as there are block columns. we can do the same
209 // thing as above, but we have to do it multiple times
211 0,
212 matrix.m(),
213 [this, &matrix](const size_type row_begin, const size_type row_end) {
214 for (size_type row = row_begin; row < row_end; ++row)
215 {
216 long int cursor = Ap[row];
217 for (size_type block = 0; block < matrix.n_block_cols(); ++block)
218 {
219 // find the next out-of-order element
220 while ((cursor < Ap[row + 1] - 1) &&
221 (Ai[cursor] < Ai[cursor + 1]))
222 ++cursor;
223
224 // if there is none, then just go on
225 if (cursor == Ap[row + 1] - 1)
226 break;
227
228 // otherwise swap this entry with successive ones as long as
229 // necessary
230 long int element = cursor;
231 while ((element < Ap[row + 1] - 1) &&
232 (Ai[element] > Ai[element + 1]))
233 {
234 std::swap(Ai[element], Ai[element + 1]);
235
236 std::swap(Ax[element], Ax[element + 1]);
237 if (numbers::NumberTraits<number>::is_complex == true)
238 std::swap(Az[element], Az[element + 1]);
239
240 ++element;
241 }
242 }
243 }
244 },
245 parallel_grainsize(matrix));
246}
247
248
249
250template <class Matrix>
251void
253{
254 Assert(matrix.m() == matrix.n(), ExcNotQuadratic());
255
256 clear();
257
258 using number = typename Matrix::value_type;
259
260 n_rows = matrix.m();
261 n_cols = matrix.n();
262
263 const size_type N = matrix.m();
264
265 // copy over the data from the matrix to the data structures UMFPACK
266 // wants. note two things: first, UMFPACK wants compressed column storage
267 // whereas we always do compressed row storage; we work around this by,
268 // rather than shuffling things around, copy over the data we have, but
269 // then call the umfpack_dl_solve function with the UMFPACK_At argument,
270 // meaning that we want to solve for the transpose system
271 //
272 // second: the data we have in the sparse matrices is "almost" right
273 // already; UMFPACK wants the entries in each row (i.e. really: column)
274 // to be sorted in ascending order. we almost have that, except that we
275 // usually store the diagonal first in each row to allow for some
276 // optimizations. thus, we have to resort things a little bit, but only
277 // within each row
278 //
279 // final note: if the matrix has entries in the sparsity pattern that are
280 // actually occupied by entries that have a zero numerical value, then we
281 // keep them anyway. people are supposed to provide accurate sparsity
282 // patterns.
283 Ap.resize(N + 1);
284 Ai.resize(matrix.n_nonzero_elements());
285 Ax.resize(matrix.n_nonzero_elements());
287 Az.resize(matrix.n_nonzero_elements());
288
289 // first fill row lengths array
290 Ap[0] = 0;
291 for (size_type row = 1; row <= N; ++row)
292 Ap[row] = Ap[row - 1] + matrix.get_row_length(row - 1);
293 Assert(static_cast<size_type>(Ap.back()) == Ai.size(), ExcInternalError());
294
295 // then copy over matrix elements. note that for sparse matrices,
296 // iterators are sorted so that they traverse each row from start to end
297 // before moving on to the next row.
299 0,
300 matrix.m(),
301 [this, &matrix](const size_type row_begin, const size_type row_end) {
302 for (size_type row = row_begin; row < row_end; ++row)
303 {
304 long int index = Ap[row];
305 for (typename Matrix::const_iterator p = matrix.begin(row);
306 p != matrix.end(row);
307 ++p)
308 {
309 // write entry into the first free one for this row
310 Ai[index] = p->column();
311 Ax[index] = std::real(p->value());
312 if (numbers::NumberTraits<number>::is_complex == true)
313 Az[index] = std::imag(p->value());
314
315 // then move pointer ahead
316 ++index;
317 }
318 Assert(index == Ap[row + 1], ExcInternalError());
319 }
320 },
321 parallel_grainsize(matrix));
322
323 // make sure that the elements in each row are sorted. we have to be more
324 // careful for block sparse matrices, so ship this task out to a
325 // different function
326 sort_arrays(matrix);
327
328 int status;
330 status = umfpack_dl_symbolic(N,
331 N,
332 Ap.data(),
333 Ai.data(),
334 Ax.data(),
335 &symbolic_decomposition,
336 control.data(),
337 nullptr);
338 else
339 status = umfpack_zl_symbolic(N,
340 N,
341 Ap.data(),
342 Ai.data(),
343 Ax.data(),
344 Az.data(),
345 &symbolic_decomposition,
346 control.data(),
347 nullptr);
348 AssertThrow(status == UMFPACK_OK,
349 ExcUMFPACKError("umfpack_dl_symbolic", status));
350
352 status = umfpack_dl_numeric(Ap.data(),
353 Ai.data(),
354 Ax.data(),
355 symbolic_decomposition,
356 &numeric_decomposition,
357 control.data(),
358 nullptr);
359 else
360 status = umfpack_zl_numeric(Ap.data(),
361 Ai.data(),
362 Ax.data(),
363 Az.data(),
364 symbolic_decomposition,
365 &numeric_decomposition,
366 control.data(),
367 nullptr);
368 AssertThrow(status == UMFPACK_OK,
369 ExcUMFPACKError("umfpack_dl_numeric", status));
370
371 umfpack_dl_free_symbolic(&symbolic_decomposition);
372}
373
374
375
376void
378 const bool transpose /*=false*/) const
379{
380 // make sure that some kind of factorize() call has happened before
381 Assert(Ap.size() != 0, ExcNotInitialized());
382 Assert(Ai.size() != 0, ExcNotInitialized());
383 Assert(Ai.size() == Ax.size(), ExcNotInitialized());
384
385 Assert(Az.empty(),
386 ExcMessage("You have previously factored a matrix using this class "
387 "that had complex-valued entries. This then requires "
388 "applying the factored matrix to a complex-valued "
389 "vector, but you are only providing a real-valued vector "
390 "here."));
391
392 Vector<double> rhs(rhs_and_solution.size());
393 rhs = rhs_and_solution;
394
395 // solve the system. note that since UMFPACK wants compressed column
396 // storage instead of the compressed row storage format we use in
397 // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
398
399 // Conversely, if we solve for the transpose, we have to use UMFPACK_A
400 // instead.
401 const int status = umfpack_dl_solve(transpose ? UMFPACK_A : UMFPACK_At,
402 Ap.data(),
403 Ai.data(),
404 Ax.data(),
405 rhs_and_solution.begin(),
406 rhs.begin(),
408 control.data(),
409 nullptr);
410 AssertThrow(status == UMFPACK_OK,
411 ExcUMFPACKError("umfpack_dl_solve", status));
412}
413
414
415
416void
417SparseDirectUMFPACK::solve(Vector<std::complex<double>> &rhs_and_solution,
418 const bool transpose /*=false*/) const
419{
420# ifdef DEAL_II_WITH_COMPLEX_VALUES
421 // make sure that some kind of factorize() call has happened before
422 Assert(Ap.size() != 0, ExcNotInitialized());
423 Assert(Ai.size() != 0, ExcNotInitialized());
424 Assert(Ai.size() == Ax.size(), ExcNotInitialized());
425
426 // First see whether the matrix that was factorized was complex-valued.
427 // If so, just apply the complex factorization to the vector.
428 if (Az.size() != 0)
429 {
430 Assert(Ax.size() == Az.size(), ExcInternalError());
431
432 // It would be nice if we could just present a pointer to the
433 // first element of the complex-valued solution vector and let
434 // UMFPACK fill both the real and imaginary parts of the solution
435 // vector at that address. UMFPACK calls this the 'packed' format,
436 // and in those cases it only takes one pointer to the entire
437 // vector, rather than a pointer to the real and one pointer to
438 // the imaginary parts of the vector. The problem is that if we
439 // want to pack, then we also need to pack the matrix, and the
440 // functions above have already decided that we don't want to pack
441 // the matrix but instead deal with split format for the matrix,
442 // and then UMFPACK decides that it can't deal with a split matrix
443 // and a packed vector. We have to choose one or the other, not
444 // mix.
445 //
446 // So create four vectors, one each for the real and imaginary parts
447 // of the right hand side and solution.
448 Vector<double> rhs_re(rhs_and_solution.size());
449 Vector<double> rhs_im(rhs_and_solution.size());
450 for (unsigned int i = 0; i < rhs_and_solution.size(); ++i)
451 {
452 rhs_re(i) = std::real(rhs_and_solution(i));
453 rhs_im(i) = std::imag(rhs_and_solution(i));
454 }
455
456 Vector<double> solution_re(rhs_and_solution.size());
457 Vector<double> solution_im(rhs_and_solution.size());
458
459 // Solve the system. note that since UMFPACK wants compressed column
460 // storage instead of the compressed row storage format we use in
461 // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
462 //
463 // Conversely, if we solve for the transpose, we have to use UMFPACK_A
464 // instead.
465 //
466 // Note that for the complex case here, the transpose is selected using
467 // UMFPACK_Aat, not UMFPACK_At.
468 const int status = umfpack_zl_solve(transpose ? UMFPACK_A : UMFPACK_Aat,
469 Ap.data(),
470 Ai.data(),
471 Ax.data(),
472 Az.data(),
473 solution_re.data(),
474 solution_im.data(),
475 rhs_re.data(),
476 rhs_im.data(),
478 control.data(),
479 nullptr);
480 AssertThrow(status == UMFPACK_OK,
481 ExcUMFPACKError("umfpack_dl_solve", status));
482
483 // Now put things back together into the output vector
484 for (unsigned int i = 0; i < rhs_and_solution.size(); ++i)
485 rhs_and_solution(i) = {solution_re(i), solution_im(i)};
486 }
487 else
488 {
489 // We have factorized a real-valued matrix, but the rhs and solution
490 // vectors are complex-valued. UMFPACK does not natively support this
491 // case, but we can just apply the factorization to real and imaginary
492 // parts of the right hand side separately
493 const Vector<std::complex<double>> rhs = rhs_and_solution;
494
495 // Get the real part of the right hand side, solve with it, and copy the
496 // results into the result vector by just copying the real output
497 // into the complex-valued result vector (which implies setting the
498 // imaginary parts to zero):
499 Vector<double> rhs_real_or_imag(rhs_and_solution.size());
500 for (unsigned int i = 0; i < rhs.size(); ++i)
501 rhs_real_or_imag(i) = std::real(rhs(i));
502
503 solve(rhs_real_or_imag, transpose);
504
505 rhs_and_solution = rhs_real_or_imag;
506
507 // Then repeat the whole thing with the imaginary part. The copying step
508 // is more complicated now because we can only touch the imaginary
509 // component of the output vector:
510 for (unsigned int i = 0; i < rhs.size(); ++i)
511 rhs_real_or_imag(i) = std::imag(rhs(i));
512
513 solve(rhs_real_or_imag, transpose);
514
515 for (unsigned int i = 0; i < rhs.size(); ++i)
516 rhs_and_solution(i).imag(rhs_real_or_imag(i));
517 }
518
519# else
520
521 (void)rhs_and_solution;
522 (void)transpose;
523 Assert(false,
525 "This function can't be called if deal.II has been configured "
526 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
527# endif
528}
529
530
531void
533 const bool transpose /*=false*/) const
534{
535 // the UMFPACK functions want a contiguous array of elements, so
536 // there is no way around copying data around. thus, just copy the
537 // data into a regular vector and back
538 Vector<double> tmp(rhs_and_solution.size());
539 tmp = rhs_and_solution;
540 solve(tmp, transpose);
541 rhs_and_solution = tmp;
542}
543
544
545
546void
547SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &rhs_and_solution,
548 const bool transpose /*=false*/) const
549{
550# ifdef DEAL_II_WITH_COMPLEX_VALUES
551 // the UMFPACK functions want a contiguous array of elements, so
552 // there is no way around copying data around. thus, just copy the
553 // data into a regular vector and back
554 Vector<std::complex<double>> tmp(rhs_and_solution.size());
555 tmp = rhs_and_solution;
556 solve(tmp, transpose);
557 rhs_and_solution = tmp;
558
559# else
560 (void)rhs_and_solution;
561 (void)transpose;
562 Assert(false,
564 "This function can't be called if deal.II has been configured "
565 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
566# endif
567}
568
569
570
571template <class Matrix>
572void
573SparseDirectUMFPACK::solve(const Matrix &matrix,
574 Vector<double> &rhs_and_solution,
575 const bool transpose /*=false*/)
576{
577 factorize(matrix);
578 solve(rhs_and_solution, transpose);
579}
580
581
582
583template <class Matrix>
584void
585SparseDirectUMFPACK::solve(const Matrix &matrix,
586 Vector<std::complex<double>> &rhs_and_solution,
587 const bool transpose /*=false*/)
588{
589# ifdef DEAL_II_WITH_COMPLEX_VALUES
590 factorize(matrix);
591 solve(rhs_and_solution, transpose);
592
593# else
594
595 (void)matrix;
596 (void)rhs_and_solution;
597 (void)transpose;
598 Assert(false,
600 "This function can't be called if deal.II has been configured "
601 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
602# endif
603}
604
605
606
607template <class Matrix>
608void
609SparseDirectUMFPACK::solve(const Matrix &matrix,
610 BlockVector<double> &rhs_and_solution,
611 const bool transpose /*=false*/)
612{
613 factorize(matrix);
614 solve(rhs_and_solution, transpose);
615}
616
617
618
619template <class Matrix>
620void
621SparseDirectUMFPACK::solve(const Matrix &matrix,
622 BlockVector<std::complex<double>> &rhs_and_solution,
623 const bool transpose /*=false*/)
624{
625# ifdef DEAL_II_WITH_COMPLEX_VALUES
626 factorize(matrix);
627 solve(rhs_and_solution, transpose);
628
629# else
630
631 (void)matrix;
632 (void)rhs_and_solution;
633 (void)transpose;
634 Assert(false,
636 "This function can't be called if deal.II has been configured "
637 "with DEAL_II_WITH_COMPLEX_VALUES=FALSE."));
638# endif
639}
640
641
642#else
643
644
646 : n_rows(0)
647 , n_cols(0)
648 , symbolic_decomposition(nullptr)
649 , numeric_decomposition(nullptr)
650 , control(0)
651{}
652
653
654void
656{}
657
658
659template <class Matrix>
660void
661SparseDirectUMFPACK::factorize(const Matrix &)
662{
664 false,
666 "To call this function you need UMFPACK, but you configured deal.II "
667 "without passing the necessary switch to 'cmake'. Please consult the "
668 "installation instructions at https://dealii.org/current/readme.html"));
669}
670
671
672void
674{
676 false,
678 "To call this function you need UMFPACK, but you configured deal.II "
679 "without passing the necessary switch to 'cmake'. Please consult the "
680 "installation instructions at https://dealii.org/current/readme.html"));
681}
682
683
684
685void
686SparseDirectUMFPACK::solve(Vector<std::complex<double>> &, const bool) const
687{
689 false,
691 "To call this function you need UMFPACK, but you configured deal.II "
692 "without passing the necessary switch to 'cmake'. Please consult the "
693 "installation instructions at https://dealii.org/current/readme.html"));
694}
695
696
697
698void
700{
702 false,
704 "To call this function you need UMFPACK, but you configured deal.II "
705 "without passing the necessary switch to 'cmake'. Please consult the "
706 "installation instructions at https://dealii.org/current/readme.html"));
707}
708
709
710
711void
712SparseDirectUMFPACK::solve(BlockVector<std::complex<double>> &,
713 const bool) const
714{
716 false,
718 "To call this function you need UMFPACK, but you configured deal.II "
719 "without passing the necessary switch to 'cmake'. Please consult the "
720 "installation instructions at https://dealii.org/current/readme.html"));
721}
722
723
724
725template <class Matrix>
726void
727SparseDirectUMFPACK::solve(const Matrix &, Vector<double> &, const bool)
728{
730 false,
732 "To call this function you need UMFPACK, but you configured deal.II "
733 "without passing the necessary switch to 'cmake'. Please consult the "
734 "installation instructions at https://dealii.org/current/readme.html"));
735}
736
737
738
739template <class Matrix>
740void
741SparseDirectUMFPACK::solve(const Matrix &,
742 Vector<std::complex<double>> &,
743 const bool)
744{
746 false,
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"));
751}
752
753
754
755template <class Matrix>
756void
757SparseDirectUMFPACK::solve(const Matrix &, BlockVector<double> &, const bool)
758{
760 false,
762 "To call this function you need UMFPACK, but you configured deal.II "
763 "without passing the necessary switch to 'cmake'. Please consult the "
764 "installation instructions at https://dealii.org/current/readme.html"));
765}
766
767
768
769template <class Matrix>
770void
771SparseDirectUMFPACK::solve(const Matrix &,
772 BlockVector<std::complex<double>> &,
773 const bool)
774{
776 false,
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"));
781}
782
783#endif
784
785
786template <class Matrix>
787void
789{
790 this->factorize(M);
791}
792
793
794void
796{
797 dst = src;
798 this->solve(dst);
799}
800
801
802
803void
805 const BlockVector<double> &src) const
806{
807 dst = src;
808 this->solve(dst);
809}
810
811
812void
814 const Vector<double> &src) const
815{
816 dst = src;
817 this->solve(dst, /*transpose=*/true);
818}
819
820
821
822void
824 const BlockVector<double> &src) const
825{
826 dst = src;
827 this->solve(dst, /*transpose=*/true);
828}
829
832{
834 return n_rows;
835}
836
839{
841 return n_cols;
842}
843
844
845// explicit instantiations for SparseMatrixUMFPACK
846#define InstantiateUMFPACK(MatrixType) \
847 template void SparseDirectUMFPACK::factorize(const MatrixType &); \
848 template void SparseDirectUMFPACK::solve(const MatrixType &, \
849 Vector<double> &, \
850 const bool); \
851 template void SparseDirectUMFPACK::solve(const MatrixType &, \
852 Vector<std::complex<double>> &, \
853 const bool); \
854 template void SparseDirectUMFPACK::solve(const MatrixType &, \
855 BlockVector<double> &, \
856 const bool); \
857 template void SparseDirectUMFPACK::solve( \
858 const MatrixType &, BlockVector<std::complex<double>> &, const bool); \
859 template void SparseDirectUMFPACK::initialize(const MatrixType &, \
860 const AdditionalData)
861
862// Instantiate everything for real-valued matrices
869
870// Now also for complex-valued matrices
871#ifdef DEAL_II_WITH_COMPLEX_VALUES
872InstantiateUMFPACK(SparseMatrix<std::complex<double>>);
873InstantiateUMFPACK(SparseMatrix<std::complex<float>>);
876#endif
877
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:452
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
#define InstantiateUMFPACK(MatrixType)