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