Reference documentation for deal.II version 9.2.0
\(\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\}}\)
petsc_matrix_base.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 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 
18 #ifdef DEAL_II_WITH_PETSC
19 
20 # include <deal.II/lac/exceptions.h>
25 
27 
28 namespace PETScWrappers
29 {
30  namespace MatrixIterators
31  {
32 # ifndef DOXYGEN
33  void
35  {
36  // if we are asked to visit the past-the-end line (or a line that is not
37  // stored on the current processor), then simply release all our caches
38  // and go on with life
39  if (matrix->in_local_range(this->a_row) == false)
40  {
41  colnum_cache.reset();
42  value_cache.reset();
43 
44  return;
45  }
46 
47  // get a representation of the present row
48  PetscInt ncols;
49  const PetscInt * colnums;
50  const PetscScalar *values;
51 
52  PetscErrorCode ierr =
53  MatGetRow(*matrix, this->a_row, &ncols, &colnums, &values);
54  AssertThrow(ierr == 0, ExcPETScError(ierr));
55 
56  // copy it into our caches if the line
57  // isn't empty. if it is, then we've
58  // done something wrong, since we
59  // shouldn't have initialized an
60  // iterator for an empty line (what
61  // would it point to?)
62  Assert(ncols != 0, ExcInternalError());
63  colnum_cache =
64  std::make_shared<std::vector<size_type>>(colnums, colnums + ncols);
65  value_cache =
66  std::make_shared<std::vector<PetscScalar>>(values, values + ncols);
67 
68  // and finally restore the matrix
69  ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values);
70  AssertThrow(ierr == 0, ExcPETScError(ierr));
71  }
72 # endif
73  } // namespace MatrixIterators
74 
75 
76 
78  : matrix(nullptr)
79  , last_action(VectorOperation::unknown)
80  {}
81 
82 
83 
85  {
87  }
88 
89 
90 
91  void
93  {
94  // destroy the matrix...
95  {
96  const PetscErrorCode ierr = destroy_matrix(matrix);
97  AssertThrow(ierr == 0, ExcPETScError(ierr));
98  }
99 
100  // ...and replace it by an empty
101  // sequential matrix
102  const int m = 0, n = 0, n_nonzero_per_row = 0;
103  const PetscErrorCode ierr = MatCreateSeqAIJ(
104  PETSC_COMM_SELF, m, n, n_nonzero_per_row, nullptr, &matrix);
105  AssertThrow(ierr == 0, ExcPETScError(ierr));
106  }
107 
108 
109 
110  MatrixBase &
112  {
113  (void)d;
115 
117 
118  const PetscErrorCode ierr = MatZeroEntries(matrix);
119  AssertThrow(ierr == 0, ExcPETScError(ierr));
120 
121  return *this;
122  }
123 
124 
125 
126  void
127  MatrixBase::clear_row(const size_type row, const PetscScalar new_diag_value)
128  {
129  std::vector<size_type> rows(1, row);
130  clear_rows(rows, new_diag_value);
131  }
132 
133 
134 
135  void
136  MatrixBase::clear_rows(const std::vector<size_type> &rows,
137  const PetscScalar new_diag_value)
138  {
140 
141  // now set all the entries of these rows
142  // to zero
143  const std::vector<PetscInt> petsc_rows(rows.begin(), rows.end());
144 
145  // call the functions. note that we have
146  // to call them even if #rows is empty,
147  // since this is a collective operation
148  IS index_set;
149 
150  ISCreateGeneral(get_mpi_communicator(),
151  rows.size(),
152  petsc_rows.data(),
153  PETSC_COPY_VALUES,
154  &index_set);
155 
156  const PetscErrorCode ierr =
157  MatZeroRowsIS(matrix, index_set, new_diag_value, nullptr, nullptr);
158  AssertThrow(ierr == 0, ExcPETScError(ierr));
159  ISDestroy(&index_set);
160  }
161 
162 
163 
164  PetscScalar
165  MatrixBase::el(const size_type i, const size_type j) const
166  {
167  PetscInt petsc_i = i, petsc_j = j;
168 
169  PetscScalar value;
170 
171  const PetscErrorCode ierr =
172  MatGetValues(matrix, 1, &petsc_i, 1, &petsc_j, &value);
173  AssertThrow(ierr == 0, ExcPETScError(ierr));
174 
175  return value;
176  }
177 
178 
179 
180  PetscScalar
182  {
183  Assert(m() == n(), ExcNotQuadratic());
184 
185  // this doesn't seem to work any
186  // different than any other element
187  return el(i, i);
188  }
189 
190 
191 
192  void
194  {
195  {
196 # ifdef DEBUG
197 # ifdef DEAL_II_WITH_MPI
198  // Check that all processors agree that last_action is the same (or none!)
199 
200  int my_int_last_action = last_action;
201  int all_int_last_action;
202 
203  const int ierr = MPI_Allreduce(&my_int_last_action,
204  &all_int_last_action,
205  1,
206  MPI_INT,
207  MPI_BOR,
209  AssertThrowMPI(ierr);
210 
211  AssertThrow(all_int_last_action !=
213  ExcMessage("Error: not all processors agree on the last "
214  "VectorOperation before this compress() call."));
215 # endif
216 # endif
217  }
218 
219  AssertThrow(
221  ExcMessage(
222  "Missing compress() or calling with wrong VectorOperation argument."));
223 
224  // flush buffers
225  PetscErrorCode ierr = MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY);
226  AssertThrow(ierr == 0, ExcPETScError(ierr));
227 
228  ierr = MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY);
229  AssertThrow(ierr == 0, ExcPETScError(ierr));
230 
232  }
233 
234 
235 
238  {
239  PetscInt n_rows, n_cols;
240 
241  const PetscErrorCode ierr = MatGetSize(matrix, &n_rows, &n_cols);
242  AssertThrow(ierr == 0, ExcPETScError(ierr));
243 
244  return n_rows;
245  }
246 
247 
248 
251  {
252  PetscInt n_rows, n_cols;
253 
254  const PetscErrorCode ierr = MatGetSize(matrix, &n_rows, &n_cols);
255  AssertThrow(ierr == 0, ExcPETScError(ierr));
256 
257  return n_cols;
258  }
259 
260 
261 
264  {
265  PetscInt n_rows, n_cols;
266 
267  const PetscErrorCode ierr = MatGetLocalSize(matrix, &n_rows, &n_cols);
268  AssertThrow(ierr == 0, ExcPETScError(ierr));
269 
270  return n_rows;
271  }
272 
273 
274 
275  std::pair<MatrixBase::size_type, MatrixBase::size_type>
277  {
278  PetscInt begin, end;
279 
280  const PetscErrorCode ierr =
281  MatGetOwnershipRange(static_cast<const Mat &>(matrix), &begin, &end);
282  AssertThrow(ierr == 0, ExcPETScError(ierr));
283 
284  return std::make_pair(begin, end);
285  }
286 
287 
288 
291  {
292  MatInfo mat_info;
293  const PetscErrorCode ierr = MatGetInfo(matrix, MAT_GLOBAL_SUM, &mat_info);
294  AssertThrow(ierr == 0, ExcPETScError(ierr));
295 
296  return static_cast<size_type>(mat_info.nz_used);
297  }
298 
299 
300 
303  {
304  // TODO: this function will probably only work if compress() was called on
305  // the matrix previously. however, we can't do this here, since it would
306  // impose global communication and one would have to make sure that this
307  // function is called the same number of times from all processors,
308  // something that is unreasonable. there should simply be a way in PETSc to
309  // query the number of entries in a row bypassing the call to compress(),
310  // but I can't find one
311  Assert(row < m(), ExcInternalError());
312 
313  // get a representation of the present
314  // row
315  PetscInt ncols;
316  const PetscInt * colnums;
317  const PetscScalar *values;
318 
319  // TODO: this is probably horribly inefficient; we should lobby for a way to
320  // query this information from PETSc
321  PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values);
322  AssertThrow(ierr == 0, ExcPETScError(ierr));
323 
324  // then restore the matrix and return the number of columns in this row as
325  // queried previously. Starting with PETSc 3.4, MatRestoreRow actually
326  // resets the last three arguments to nullptr, to avoid abuse of pointers
327  // now dangling. as a consequence, we need to save the size of the array
328  // and return the saved value.
329  const PetscInt ncols_saved = ncols;
330  ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values);
331  AssertThrow(ierr == 0, ExcPETScError(ierr));
332 
333  return ncols_saved;
334  }
335 
336 
337  PetscReal
339  {
340  PetscReal result;
341 
342  const PetscErrorCode ierr = MatNorm(matrix, NORM_1, &result);
343  AssertThrow(ierr == 0, ExcPETScError(ierr));
344 
345  return result;
346  }
347 
348 
349 
350  PetscReal
352  {
353  PetscReal result;
354 
355  const PetscErrorCode ierr = MatNorm(matrix, NORM_INFINITY, &result);
356  AssertThrow(ierr == 0, ExcPETScError(ierr));
357 
358  return result;
359  }
360 
361 
362 
363  PetscReal
365  {
366  PetscReal result;
367 
368  const PetscErrorCode ierr = MatNorm(matrix, NORM_FROBENIUS, &result);
369  AssertThrow(ierr == 0, ExcPETScError(ierr));
370 
371  return result;
372  }
373 
374 
375  PetscScalar
377  {
378  VectorBase tmp(v);
379  vmult(tmp, v);
380  return tmp * v;
381  }
382 
383 
384  PetscScalar
386  const VectorBase &v) const
387  {
388  VectorBase tmp(u);
389  vmult(tmp, v);
390  return u * tmp;
391  }
392 
393 
394  PetscScalar
396  {
397  PetscScalar result;
398 
399  const PetscErrorCode ierr = MatGetTrace(matrix, &result);
400  AssertThrow(ierr == 0, ExcPETScError(ierr));
401 
402  return result;
403  }
404 
405 
406 
407  MatrixBase &
408  MatrixBase::operator*=(const PetscScalar a)
409  {
410  const PetscErrorCode ierr = MatScale(matrix, a);
411  AssertThrow(ierr == 0, ExcPETScError(ierr));
412 
413  return *this;
414  }
415 
416 
417 
418  MatrixBase &
419  MatrixBase::operator/=(const PetscScalar a)
420  {
421  const PetscScalar factor = 1. / a;
422  const PetscErrorCode ierr = MatScale(matrix, factor);
423  AssertThrow(ierr == 0, ExcPETScError(ierr));
424 
425  return *this;
426  }
427 
428 
429  MatrixBase &
430  MatrixBase::add(const PetscScalar factor, const MatrixBase &other)
431  {
432  const PetscErrorCode ierr =
433  MatAXPY(matrix, factor, other, DIFFERENT_NONZERO_PATTERN);
434  AssertThrow(ierr == 0, ExcPETScError(ierr));
435 
436  return *this;
437  }
438 
439 
440 
441  MatrixBase &
442  MatrixBase::add(const MatrixBase &other, const PetscScalar factor)
443  {
444  return add(factor, other);
445  }
446 
447 
448  void
449  MatrixBase::vmult(VectorBase &dst, const VectorBase &src) const
450  {
451  Assert(&src != &dst, ExcSourceEqualsDestination());
452 
453  const PetscErrorCode ierr = MatMult(matrix, src, dst);
454  AssertThrow(ierr == 0, ExcPETScError(ierr));
455  }
456 
457 
458 
459  void
460  MatrixBase::Tvmult(VectorBase &dst, const VectorBase &src) const
461  {
462  Assert(&src != &dst, ExcSourceEqualsDestination());
463 
464  const PetscErrorCode ierr = MatMultTranspose(matrix, src, dst);
465  AssertThrow(ierr == 0, ExcPETScError(ierr));
466  }
467 
468 
469 
470  void
472  {
473  Assert(&src != &dst, ExcSourceEqualsDestination());
474 
475  const PetscErrorCode ierr = MatMultAdd(matrix, src, dst, dst);
476  AssertThrow(ierr == 0, ExcPETScError(ierr));
477  }
478 
479 
480 
481  void
483  {
484  Assert(&src != &dst, ExcSourceEqualsDestination());
485 
486  const PetscErrorCode ierr = MatMultTransposeAdd(matrix, src, dst, dst);
487  AssertThrow(ierr == 0, ExcPETScError(ierr));
488  }
489 
490 
491  namespace internals
492  {
493  void
494  perform_mmult(const MatrixBase &inputleft,
495  const MatrixBase &inputright,
496  MatrixBase & result,
497  const VectorBase &V,
498  const bool transpose_left)
499  {
500  const bool use_vector = (V.size() == inputright.m() ? true : false);
501  if (transpose_left == false)
502  {
503  Assert(inputleft.n() == inputright.m(),
504  ExcDimensionMismatch(inputleft.n(), inputright.m()));
505  }
506  else
507  {
508  Assert(inputleft.m() == inputright.m(),
509  ExcDimensionMismatch(inputleft.m(), inputright.m()));
510  }
511 
512  result.clear();
513 
514  PetscErrorCode ierr;
515 
516  if (use_vector == false)
517  {
518  if (transpose_left)
519  {
520  ierr = MatTransposeMatMult(inputleft,
521  inputright,
522  MAT_INITIAL_MATRIX,
523  PETSC_DEFAULT,
524  &result.petsc_matrix());
525  AssertThrow(ierr == 0, ExcPETScError(ierr));
526  }
527  else
528  {
529  ierr = MatMatMult(inputleft,
530  inputright,
531  MAT_INITIAL_MATRIX,
532  PETSC_DEFAULT,
533  &result.petsc_matrix());
534  AssertThrow(ierr == 0, ExcPETScError(ierr));
535  }
536  }
537  else
538  {
539  Mat tmp;
540  ierr = MatDuplicate(inputleft, MAT_COPY_VALUES, &tmp);
541  AssertThrow(ierr == 0, ExcPETScError(ierr));
542  if (transpose_left)
543  {
544 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
545  ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp);
546 # else
547  ierr = MatTranspose(tmp, MAT_INPLACE_MATRIX, &tmp);
548 # endif
549  AssertThrow(ierr == 0, ExcPETScError(ierr));
550  }
551  ierr = MatDiagonalScale(tmp, nullptr, V);
552  AssertThrow(ierr == 0, ExcPETScError(ierr));
553  ierr = MatMatMult(tmp,
554  inputright,
555  MAT_INITIAL_MATRIX,
556  PETSC_DEFAULT,
557  &result.petsc_matrix());
558  AssertThrow(ierr == 0, ExcPETScError(ierr));
559  }
560  }
561  } // namespace internals
562 
563  void
565  const MatrixBase &B,
566  const VectorBase &V) const
567  {
568  internals::perform_mmult(*this, B, C, V, false);
569  }
570 
571  void
573  const MatrixBase &B,
574  const VectorBase &V) const
575  {
576  internals::perform_mmult(*this, B, C, V, true);
577  }
578 
579  PetscScalar
581  const VectorBase &x,
582  const VectorBase &b) const
583  {
584  // avoid the use of a temporary, and
585  // rather do one negation pass more than
586  // necessary
587  vmult(dst, x);
588  dst -= b;
589  dst *= -1;
590 
591  return dst.l2_norm();
592  }
593 
594 
595 
596  MatrixBase::operator Mat() const
597  {
598  return matrix;
599  }
600 
601  Mat &
603  {
604  return matrix;
605  }
606 
607  void
609  {
610 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
611  const PetscErrorCode ierr = MatTranspose(matrix, MAT_REUSE_MATRIX, &matrix);
612 # else
613  const PetscErrorCode ierr =
614  MatTranspose(matrix, MAT_INPLACE_MATRIX, &matrix);
615 # endif
616  AssertThrow(ierr == 0, ExcPETScError(ierr));
617  }
618 
619  PetscBool
620  MatrixBase::is_symmetric(const double tolerance)
621  {
622  PetscBool truth;
624  const PetscErrorCode ierr = MatIsSymmetric(matrix, tolerance, &truth);
625  AssertThrow(ierr == 0, ExcPETScError(ierr));
626  return truth;
627  }
628 
629  PetscBool
630  MatrixBase::is_hermitian(const double tolerance)
631  {
632  PetscBool truth;
633 
635  const PetscErrorCode ierr = MatIsHermitian(matrix, tolerance, &truth);
636  AssertThrow(ierr == 0, ExcPETScError(ierr));
637 
638  return truth;
639  }
640 
641  void
642  MatrixBase::write_ascii(const PetscViewerFormat format)
643  {
645 
646  // Set options
647  PetscErrorCode ierr =
648  PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, format);
649  AssertThrow(ierr == 0, ExcPETScError(ierr));
650 
651  // Write to screen
652  ierr = MatView(matrix, PETSC_VIEWER_STDOUT_WORLD);
653  AssertThrow(ierr == 0, ExcPETScError(ierr));
654  }
655 
656  void
657  MatrixBase::print(std::ostream &out, const bool /*alternative_output*/) const
658  {
659  std::pair<MatrixBase::size_type, MatrixBase::size_type> loc_range =
660  local_range();
661 
662  PetscInt ncols;
663  const PetscInt * colnums;
664  const PetscScalar *values;
665 
667  for (row = loc_range.first; row < loc_range.second; ++row)
668  {
669  PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values);
670  AssertThrow(ierr == 0, ExcPETScError(ierr));
671 
672  for (PetscInt col = 0; col < ncols; ++col)
673  {
674  out << "(" << row << "," << colnums[col] << ") " << values[col]
675  << std::endl;
676  }
677 
678  ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values);
679  AssertThrow(ierr == 0, ExcPETScError(ierr));
680  }
681 
682  AssertThrow(out, ExcIO());
683  }
684 
685 
686 
687  std::size_t
689  {
690  MatInfo info;
691  const PetscErrorCode ierr = MatGetInfo(matrix, MAT_LOCAL, &info);
692  AssertThrow(ierr == 0, ExcPETScError(ierr));
693 
694  return sizeof(*this) + static_cast<size_type>(info.memory);
695  }
696 
697 } // namespace PETScWrappers
698 
700 
701 #endif // DEAL_II_WITH_PETSC
PetscErrorCode destroy_matrix(Mat &matrix)
void Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
Contents is actually a matrix.
PetscScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
void add(const size_type i, const size_type j, const PetscScalar value)
static ::ExceptionBase & ExcIO()
void vmult(VectorBase &dst, const VectorBase &src) const
const_iterator begin() const
PetscReal frobenius_norm() const
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
static const char V
size_type row_length(const size_type row) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
PetscBool is_symmetric(const double tolerance=1.e-12)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
PetscBool is_hermitian(const double tolerance=1.e-12)
void clear_rows(const std::vector< size_type > &rows, const PetscScalar new_diag_value=0)
void compress(const VectorOperation::values operation)
size_type n_nonzero_elements() const
MatrixBase & operator=(const MatrixBase &)=delete
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1419
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void vmult_add(VectorBase &dst, const VectorBase &src) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
PetscScalar residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const
void mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
PetscScalar diag_element(const size_type i) const
void Tvmult(VectorBase &dst, const VectorBase &src) const
PetscScalar el(const size_type i, const size_type j) const
void Tvmult_add(VectorBase &dst, const VectorBase &src) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const_iterator end() const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void print(std::ostream &out, const bool alternative_output=false) const
static ::ExceptionBase & ExcNotQuadratic()
virtual ~MatrixBase() override
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
VectorOperation::values last_action
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1707
MatrixBase & operator*=(const PetscScalar factor)
std::pair< size_type, size_type > local_range() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
PetscScalar matrix_norm_square(const VectorBase &v) const
MatrixBase & operator/=(const PetscScalar factor)
static ::ExceptionBase & ExcSourceEqualsDestination()
void perform_mmult(const MatrixBase &inputleft, const MatrixBase &inputright, MatrixBase &result, const VectorBase &V, const bool transpose_left)
static const bool value
virtual const MPI_Comm & get_mpi_communicator() const =0
void clear_row(const size_type row, const PetscScalar new_diag_value=0)
std::size_t memory_consumption() const
static ::ExceptionBase & ExcInternalError()