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