Reference documentation for deal.II version Git a39983f358 2020-01-20 20:47:10 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
sparse_direct.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 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 
16 #include <deal.II/base/memory_consumption.h>
17 
18 #include <deal.II/lac/block_sparse_matrix.h>
19 #include <deal.II/lac/sparse_direct.h>
20 #include <deal.II/lac/sparse_matrix.h>
21 #include <deal.II/lac/vector.h>
22 
23 #include <cerrno>
24 #include <iostream>
25 #include <list>
26 #include <typeinfo>
27 #include <vector>
28 
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 
33 // include UMFPACK file.
34 #ifdef DEAL_II_WITH_UMFPACK
35 # include <umfpack.h>
36 #endif
37 
38 
39 
41 {
42  clear();
43 }
44 
45 
46 void
48 {}
49 
50 
51 #ifdef DEAL_II_WITH_UMFPACK
52 
54  : n_rows(0)
55  , n_cols(0)
56  , symbolic_decomposition(nullptr)
57  , numeric_decomposition(nullptr)
58  , control(UMFPACK_CONTROL)
59 {
60  umfpack_dl_defaults(control.data());
61 }
62 
63 
64 
65 void
67 {
68  // delete objects that haven't been deleted yet
69  if (symbolic_decomposition != nullptr)
70  {
71  umfpack_dl_free_symbolic(&symbolic_decomposition);
72  symbolic_decomposition = nullptr;
73  }
74 
75  if (numeric_decomposition != nullptr)
76  {
77  umfpack_dl_free_numeric(&numeric_decomposition);
78  numeric_decomposition = nullptr;
79  }
80 
81  {
82  std::vector<long int> tmp;
83  tmp.swap(Ap);
84  }
85 
86  {
87  std::vector<long int> tmp;
88  tmp.swap(Ai);
89  }
90 
91  {
92  std::vector<double> tmp;
93  tmp.swap(Ax);
94  }
95 
96  umfpack_dl_defaults(control.data());
97 }
98 
99 
100 
101 template <typename number>
102 void
104 {
105  // do the copying around of entries so that the diagonal entry is in the
106  // right place. note that this is easy to detect: since all entries apart
107  // from the diagonal entry are sorted, we know that the diagonal entry is
108  // in the wrong place if and only if its column index is larger than the
109  // column index of the second entry in a row
110  //
111  // ignore rows with only one or no entry
112  for (size_type row = 0; row < matrix.m(); ++row)
113  {
114  // we may have to move some elements that are left of the diagonal
115  // but presently after the diagonal entry to the left, whereas the
116  // diagonal entry has to move to the right. we could first figure out
117  // where to move everything to, but for simplicity we just make a
118  // series of swaps instead (this is kind of a single run of
119  // bubble-sort, which gives us the desired result since the array is
120  // already "almost" sorted)
121  //
122  // in the first loop, the condition in the while-header also checks
123  // that the row has at least two entries and that the diagonal entry
124  // is really in the wrong place
125  long int cursor = Ap[row];
126  while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
127  {
128  std::swap(Ai[cursor], Ai[cursor + 1]);
129  std::swap(Ax[cursor], Ax[cursor + 1]);
130  ++cursor;
131  }
132  }
133 }
134 
135 
136 
137 template <typename number>
138 void
140 {
141  // same thing for SparseMatrixEZ
142  for (size_type row = 0; row < matrix.m(); ++row)
143  {
144  long int cursor = Ap[row];
145  while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
146  {
147  std::swap(Ai[cursor], Ai[cursor + 1]);
148  std::swap(Ax[cursor], Ax[cursor + 1]);
149  ++cursor;
150  }
151  }
152 }
153 
154 
155 
156 template <typename number>
157 void
159 {
160  // the case for block matrices is a bit more difficult, since all we know
161  // is that *within each block*, the diagonal of that block may come
162  // first. however, that means that there may be as many entries per row
163  // in the wrong place as there are block columns. we can do the same
164  // thing as above, but we have to do it multiple times
165  for (size_type row = 0; row < matrix.m(); ++row)
166  {
167  long int cursor = Ap[row];
168  for (size_type block = 0; block < matrix.n_block_cols(); ++block)
169  {
170  // find the next out-of-order element
171  while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] < Ai[cursor + 1]))
172  ++cursor;
173 
174  // if there is none, then just go on
175  if (cursor == Ap[row + 1] - 1)
176  break;
177 
178  // otherwise swap this entry with successive ones as long as
179  // necessary
180  long int element = cursor;
181  while ((element < Ap[row + 1] - 1) && (Ai[element] > Ai[element + 1]))
182  {
183  std::swap(Ai[element], Ai[element + 1]);
184  std::swap(Ax[element], Ax[element + 1]);
185  ++element;
186  }
187  }
188  }
189 }
190 
191 
192 
193 template <class Matrix>
194 void
195 SparseDirectUMFPACK::factorize(const Matrix &matrix)
196 {
197  Assert(matrix.m() == matrix.n(), ExcNotQuadratic())
198 
199  clear();
200 
201  n_rows = matrix.m();
202  n_cols = matrix.n();
203 
204  const size_type N = matrix.m();
205 
206  // copy over the data from the matrix to the data structures UMFPACK
207  // wants. note two things: first, UMFPACK wants compressed column storage
208  // whereas we always do compressed row storage; we work around this by,
209  // rather than shuffling things around, copy over the data we have, but
210  // then call the umfpack_dl_solve function with the UMFPACK_At argument,
211  // meaning that we want to solve for the transpose system
212  //
213  // second: the data we have in the sparse matrices is "almost" right
214  // already; UMFPACK wants the entries in each row (i.e. really: column)
215  // to be sorted in ascending order. we almost have that, except that we
216  // usually store the diagonal first in each row to allow for some
217  // optimizations. thus, we have to resort things a little bit, but only
218  // within each row
219  //
220  // final note: if the matrix has entries in the sparsity pattern that are
221  // actually occupied by entries that have a zero numerical value, then we
222  // keep them anyway. people are supposed to provide accurate sparsity
223  // patterns.
224  Ap.resize(N + 1);
225  Ai.resize(matrix.n_nonzero_elements());
226  Ax.resize(matrix.n_nonzero_elements());
227 
228  // first fill row lengths array
229  Ap[0] = 0;
230  for (size_type row = 1; row <= N; ++row)
231  Ap[row] = Ap[row - 1] + matrix.get_row_length(row - 1);
232  Assert(static_cast<size_type>(Ap.back()) == Ai.size(), ExcInternalError());
233 
234  // then copy over matrix elements. note that for sparse matrices,
235  // iterators are sorted so that they traverse each row from start to end
236  // before moving on to the next row. however, this isn't true for block
237  // matrices, so we have to do a bit of book keeping
238  {
239  // have an array that for each row points to the first entry not yet
240  // written to
241  std::vector<long int> row_pointers = Ap;
242 
243  // loop over the elements of the matrix row by row, as suggested in the
244  // documentation of the sparse matrix iterator class
245  for (size_type row = 0; row < matrix.m(); ++row)
246  {
247  for (typename Matrix::const_iterator p = matrix.begin(row);
248  p != matrix.end(row);
249  ++p)
250  {
251  // write entry into the first free one for this row
252  Ai[row_pointers[row]] = p->column();
253  Ax[row_pointers[row]] = p->value();
254 
255  // then move pointer ahead
256  ++row_pointers[row];
257  }
258  }
259 
260  // at the end, we should have written all rows completely
261  for (size_type i = 0; i < Ap.size() - 1; ++i)
262  Assert(row_pointers[i] == Ap[i + 1], ExcInternalError());
263  }
264 
265  // make sure that the elements in each row are sorted. we have to be more
266  // careful for block sparse matrices, so ship this task out to a
267  // different function
268  sort_arrays(matrix);
269 
270  int status;
271  status = umfpack_dl_symbolic(N,
272  N,
273  Ap.data(),
274  Ai.data(),
275  Ax.data(),
277  control.data(),
278  nullptr);
279  AssertThrow(status == UMFPACK_OK,
280  ExcUMFPACKError("umfpack_dl_symbolic", status));
281 
282  status = umfpack_dl_numeric(Ap.data(),
283  Ai.data(),
284  Ax.data(),
286  &numeric_decomposition,
287  control.data(),
288  nullptr);
289  AssertThrow(status == UMFPACK_OK,
290  ExcUMFPACKError("umfpack_dl_numeric", status));
291 
292  umfpack_dl_free_symbolic(&symbolic_decomposition);
293 }
294 
295 
296 
297 void
299  const bool transpose /*=false*/) const
300 {
301  // make sure that some kind of factorize() call has happened before
302  Assert(Ap.size() != 0, ExcNotInitialized());
303  Assert(Ai.size() != 0, ExcNotInitialized());
304  Assert(Ai.size() == Ax.size(), ExcNotInitialized());
305 
306  Vector<double> rhs(rhs_and_solution.size());
307  rhs = rhs_and_solution;
308 
309  // solve the system. note that since UMFPACK wants compressed column
310  // storage instead of the compressed row storage format we use in
311  // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead
312 
313  // Conversely, if we solve for the transpose, we have to use UMFPACK_A
314  // instead.
315  const int status = umfpack_dl_solve(transpose ? UMFPACK_A : UMFPACK_At,
316  Ap.data(),
317  Ai.data(),
318  Ax.data(),
319  rhs_and_solution.begin(),
320  rhs.begin(),
321  numeric_decomposition,
322  control.data(),
323  nullptr);
324  AssertThrow(status == UMFPACK_OK,
325  ExcUMFPACKError("umfpack_dl_solve", status));
326 }
327 
328 
329 void
331  const bool transpose /*=false*/) const
332 {
333  // the UMFPACK functions want a contiguous array of elements, so
334  // there is no way around copying data around. thus, just copy the
335  // data into a regular vector and back
336  Vector<double> tmp(rhs_and_solution.size());
337  tmp = rhs_and_solution;
338  solve(tmp, transpose);
339  rhs_and_solution = tmp;
340 }
341 
342 
343 
344 template <class Matrix>
345 void
346 SparseDirectUMFPACK::solve(const Matrix & matrix,
347  Vector<double> &rhs_and_solution,
348  const bool transpose /*=false*/)
349 {
350  factorize(matrix);
351  solve(rhs_and_solution, transpose);
352 }
353 
354 
355 template <class Matrix>
356 void
357 SparseDirectUMFPACK::solve(const Matrix & matrix,
358  BlockVector<double> &rhs_and_solution,
359  const bool transpose /*=false*/)
360 {
361  factorize(matrix);
362  solve(rhs_and_solution, transpose);
363 }
364 
365 
366 #else
367 
368 
370  : n_rows(0)
371  , n_cols(0)
372  , symbolic_decomposition(nullptr)
373  , numeric_decomposition(nullptr)
374  , control(0)
375 {}
376 
377 
378 void
380 {}
381 
382 
383 template <class Matrix>
384 void
385 SparseDirectUMFPACK::factorize(const Matrix &)
386 {
387  AssertThrow(
388  false,
389  ExcMessage(
390  "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
391 }
392 
393 
394 void
395 SparseDirectUMFPACK::solve(Vector<double> &, const bool) const
396 {
397  AssertThrow(
398  false,
399  ExcMessage(
400  "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
401 }
402 
403 
404 
405 void
407 {
408  AssertThrow(
409  false,
410  ExcMessage(
411  "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
412 }
413 
414 
415 template <class Matrix>
416 void
417 SparseDirectUMFPACK::solve(const Matrix &, Vector<double> &, const bool)
418 {
419  AssertThrow(
420  false,
421  ExcMessage(
422  "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
423 }
424 
425 
426 
427 template <class Matrix>
428 void
429 SparseDirectUMFPACK::solve(const Matrix &, BlockVector<double> &, const bool)
430 {
431  AssertThrow(
432  false,
433  ExcMessage(
434  "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
435 }
436 
437 #endif
438 
439 
440 template <class Matrix>
441 void
443 {
444  this->factorize(M);
445 }
446 
447 
448 void
450 {
451  dst = src;
452  this->solve(dst);
453 }
454 
455 
456 
457 void
459  const BlockVector<double> &src) const
460 {
461  dst = src;
462  this->solve(dst);
463 }
464 
465 
466 void
468  const Vector<double> &src) const
469 {
470  dst = src;
471  this->solve(dst, /*transpose=*/true);
472 }
473 
474 
475 
476 void
478  const BlockVector<double> &src) const
479 {
480  dst = src;
481  this->solve(dst, /*transpose=*/true);
482 }
483 
486 {
488  return n_rows;
489 }
490 
493 {
495  return n_cols;
496 }
497 
498 
499 // explicit instantiations for SparseMatrixUMFPACK
500 #define InstantiateUMFPACK(MatrixType) \
501  template void SparseDirectUMFPACK::factorize(const MatrixType &); \
502  template void SparseDirectUMFPACK::solve(const MatrixType &, \
503  Vector<double> &, \
504  bool); \
505  template void SparseDirectUMFPACK::solve(const MatrixType &, \
506  BlockVector<double> &, \
507  bool); \
508  template void SparseDirectUMFPACK::initialize(const MatrixType &, \
509  const AdditionalData)
510 
511 InstantiateUMFPACK(SparseMatrix<double>);
512 InstantiateUMFPACK(SparseMatrix<float>);
513 InstantiateUMFPACK(SparseMatrixEZ<double>);
514 InstantiateUMFPACK(SparseMatrixEZ<float>);
515 InstantiateUMFPACK(BlockSparseMatrix<double>);
516 InstantiateUMFPACK(BlockSparseMatrix<float>);
517 
518 
519 DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcUMFPACKError(std::string arg1, int arg2)
size_type n() const
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
void factorize(const Matrix &matrix)
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
types::global_dof_index size_type
Definition: sparse_direct.h:96
iterator begin()
std::vector< double > control
static ::ExceptionBase & ExcNotQuadratic()
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
void initialize(const SparsityPattern &sparsity_pattern)
size_type m() const
std::vector< types::suitesparse_index > Ap
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()