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