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