Reference documentation for deal.II version GIT 5dcc62ab57 2022-07-04 21:05:02+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\}}\)
petsc_matrix_base.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 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 
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 
289  std::uint64_t
291  {
292  MatInfo mat_info;
293  const PetscErrorCode ierr = MatGetInfo(matrix, MAT_GLOBAL_SUM, &mat_info);
294  AssertThrow(ierr == 0, ExcPETScError(ierr));
295 
296  // MatInfo logs quantities as PetscLogDouble. So we need to cast it to match
297  // our interface.
298  return static_cast<std::uint64_t>(mat_info.nz_used);
299  }
300 
301 
302 
305  {
306  // TODO: this function will probably only work if compress() was called on
307  // the matrix previously. however, we can't do this here, since it would
308  // impose global communication and one would have to make sure that this
309  // function is called the same number of times from all processors,
310  // something that is unreasonable. there should simply be a way in PETSc to
311  // query the number of entries in a row bypassing the call to compress(),
312  // but I can't find one
313  Assert(row < m(), ExcInternalError());
314 
315  // get a representation of the present
316  // row
317  PetscInt ncols;
318  const PetscInt * colnums;
319  const PetscScalar *values;
320 
321  // TODO: this is probably horribly inefficient; we should lobby for a way to
322  // query this information from PETSc
323  PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values);
324  AssertThrow(ierr == 0, ExcPETScError(ierr));
325 
326  // then restore the matrix and return the number of columns in this row as
327  // queried previously. Starting with PETSc 3.4, MatRestoreRow actually
328  // resets the last three arguments to nullptr, to avoid abuse of pointers
329  // now dangling. as a consequence, we need to save the size of the array
330  // and return the saved value.
331  const PetscInt ncols_saved = ncols;
332  ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values);
333  AssertThrow(ierr == 0, ExcPETScError(ierr));
334 
335  return ncols_saved;
336  }
337 
338 
339  PetscReal
341  {
342  PetscReal result;
343 
344  const PetscErrorCode ierr = MatNorm(matrix, NORM_1, &result);
345  AssertThrow(ierr == 0, ExcPETScError(ierr));
346 
347  return result;
348  }
349 
350 
351 
352  PetscReal
354  {
355  PetscReal result;
356 
357  const PetscErrorCode ierr = MatNorm(matrix, NORM_INFINITY, &result);
358  AssertThrow(ierr == 0, ExcPETScError(ierr));
359 
360  return result;
361  }
362 
363 
364 
365  PetscReal
367  {
368  PetscReal result;
369 
370  const PetscErrorCode ierr = MatNorm(matrix, NORM_FROBENIUS, &result);
371  AssertThrow(ierr == 0, ExcPETScError(ierr));
372 
373  return result;
374  }
375 
376 
377  PetscScalar
379  {
380  VectorBase tmp(v);
381  vmult(tmp, v);
382  return tmp * v;
383  }
384 
385 
386  PetscScalar
388  const VectorBase &v) const
389  {
390  VectorBase tmp(u);
391  vmult(tmp, v);
392  return u * tmp;
393  }
394 
395 
396  PetscScalar
398  {
399  PetscScalar result;
400 
401  const PetscErrorCode ierr = MatGetTrace(matrix, &result);
402  AssertThrow(ierr == 0, ExcPETScError(ierr));
403 
404  return result;
405  }
406 
407 
408 
409  MatrixBase &
410  MatrixBase::operator*=(const PetscScalar a)
411  {
412  const PetscErrorCode ierr = MatScale(matrix, a);
413  AssertThrow(ierr == 0, ExcPETScError(ierr));
414 
415  return *this;
416  }
417 
418 
419 
420  MatrixBase &
421  MatrixBase::operator/=(const PetscScalar a)
422  {
423  const PetscScalar factor = 1. / a;
424  const PetscErrorCode ierr = MatScale(matrix, factor);
425  AssertThrow(ierr == 0, ExcPETScError(ierr));
426 
427  return *this;
428  }
429 
430 
431  MatrixBase &
432  MatrixBase::add(const PetscScalar factor, const MatrixBase &other)
433  {
434  const PetscErrorCode ierr =
435  MatAXPY(matrix, factor, other, DIFFERENT_NONZERO_PATTERN);
436  AssertThrow(ierr == 0, ExcPETScError(ierr));
437 
438  return *this;
439  }
440 
441 
442  void
443  MatrixBase::vmult(VectorBase &dst, const VectorBase &src) const
444  {
445  Assert(&src != &dst, ExcSourceEqualsDestination());
446 
447  const PetscErrorCode ierr = MatMult(matrix, src, dst);
448  AssertThrow(ierr == 0, ExcPETScError(ierr));
449  }
450 
451 
452 
453  void
454  MatrixBase::Tvmult(VectorBase &dst, const VectorBase &src) const
455  {
456  Assert(&src != &dst, ExcSourceEqualsDestination());
457 
458  const PetscErrorCode ierr = MatMultTranspose(matrix, src, dst);
459  AssertThrow(ierr == 0, ExcPETScError(ierr));
460  }
461 
462 
463 
464  void
466  {
467  Assert(&src != &dst, ExcSourceEqualsDestination());
468 
469  const PetscErrorCode ierr = MatMultAdd(matrix, src, dst, dst);
470  AssertThrow(ierr == 0, ExcPETScError(ierr));
471  }
472 
473 
474 
475  void
477  {
478  Assert(&src != &dst, ExcSourceEqualsDestination());
479 
480  const PetscErrorCode ierr = MatMultTransposeAdd(matrix, src, dst, dst);
481  AssertThrow(ierr == 0, ExcPETScError(ierr));
482  }
483 
484 
485  namespace internals
486  {
487  void
488  perform_mmult(const MatrixBase &inputleft,
489  const MatrixBase &inputright,
490  MatrixBase & result,
491  const VectorBase &V,
492  const bool transpose_left)
493  {
494  const bool use_vector = (V.size() == inputright.m() ? true : false);
495  if (transpose_left == false)
496  {
497  Assert(inputleft.n() == inputright.m(),
498  ExcDimensionMismatch(inputleft.n(), inputright.m()));
499  }
500  else
501  {
502  Assert(inputleft.m() == inputright.m(),
503  ExcDimensionMismatch(inputleft.m(), inputright.m()));
504  }
505 
506  result.clear();
507 
508  PetscErrorCode ierr;
509 
510  if (use_vector == false)
511  {
512  if (transpose_left)
513  {
514  ierr = MatTransposeMatMult(inputleft,
515  inputright,
516  MAT_INITIAL_MATRIX,
517  PETSC_DEFAULT,
518  &result.petsc_matrix());
519  AssertThrow(ierr == 0, ExcPETScError(ierr));
520  }
521  else
522  {
523  ierr = MatMatMult(inputleft,
524  inputright,
525  MAT_INITIAL_MATRIX,
526  PETSC_DEFAULT,
527  &result.petsc_matrix());
528  AssertThrow(ierr == 0, ExcPETScError(ierr));
529  }
530  }
531  else
532  {
533  Mat tmp;
534  ierr = MatDuplicate(inputleft, MAT_COPY_VALUES, &tmp);
535  AssertThrow(ierr == 0, ExcPETScError(ierr));
536  if (transpose_left)
537  {
538 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
539  ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp);
540 # else
541  ierr = MatTranspose(tmp, MAT_INPLACE_MATRIX, &tmp);
542 # endif
543  AssertThrow(ierr == 0, ExcPETScError(ierr));
544  }
545  ierr = MatDiagonalScale(tmp, nullptr, V);
546  AssertThrow(ierr == 0, ExcPETScError(ierr));
547  ierr = MatMatMult(tmp,
548  inputright,
549  MAT_INITIAL_MATRIX,
550  PETSC_DEFAULT,
551  &result.petsc_matrix());
552  AssertThrow(ierr == 0, ExcPETScError(ierr));
553  ierr = PETScWrappers::destroy_matrix(tmp);
554  AssertThrow(ierr == 0, ExcPETScError(ierr));
555  }
556  }
557  } // namespace internals
558 
559  void
561  const MatrixBase &B,
562  const VectorBase &V) const
563  {
564  internals::perform_mmult(*this, B, C, V, false);
565  }
566 
567  void
569  const MatrixBase &B,
570  const VectorBase &V) const
571  {
572  internals::perform_mmult(*this, B, C, V, true);
573  }
574 
575  PetscScalar
577  const VectorBase &x,
578  const VectorBase &b) const
579  {
580  // avoid the use of a temporary, and
581  // rather do one negation pass more than
582  // necessary
583  vmult(dst, x);
584  dst -= b;
585  dst *= -1;
586 
587  return dst.l2_norm();
588  }
589 
590 
591 
592  MatrixBase::operator Mat() const
593  {
594  return matrix;
595  }
596 
597  Mat &
599  {
600  return matrix;
601  }
602 
603  void
605  {
606 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
607  const PetscErrorCode ierr = MatTranspose(matrix, MAT_REUSE_MATRIX, &matrix);
608 # else
609  const PetscErrorCode ierr =
610  MatTranspose(matrix, MAT_INPLACE_MATRIX, &matrix);
611 # endif
612  AssertThrow(ierr == 0, ExcPETScError(ierr));
613  }
614 
615  PetscBool
616  MatrixBase::is_symmetric(const double tolerance)
617  {
618  PetscBool truth;
620  const PetscErrorCode ierr = MatIsSymmetric(matrix, tolerance, &truth);
621  AssertThrow(ierr == 0, ExcPETScError(ierr));
622  return truth;
623  }
624 
625  PetscBool
626  MatrixBase::is_hermitian(const double tolerance)
627  {
628  PetscBool truth;
629 
631  const PetscErrorCode ierr = MatIsHermitian(matrix, tolerance, &truth);
632  AssertThrow(ierr == 0, ExcPETScError(ierr));
633 
634  return truth;
635  }
636 
637  void
638  MatrixBase::write_ascii(const PetscViewerFormat format)
639  {
641 
642  // Set options
643  PetscErrorCode ierr =
644  PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, format);
645  AssertThrow(ierr == 0, ExcPETScError(ierr));
646 
647  // Write to screen
648  ierr = MatView(matrix, PETSC_VIEWER_STDOUT_WORLD);
649  AssertThrow(ierr == 0, ExcPETScError(ierr));
650  }
651 
652  void
653  MatrixBase::print(std::ostream &out, const bool /*alternative_output*/) const
654  {
655  std::pair<MatrixBase::size_type, MatrixBase::size_type> loc_range =
656  local_range();
657 
658  PetscInt ncols;
659  const PetscInt * colnums;
660  const PetscScalar *values;
661 
663  for (row = loc_range.first; row < loc_range.second; ++row)
664  {
665  PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values);
666  AssertThrow(ierr == 0, ExcPETScError(ierr));
667 
668  for (PetscInt col = 0; col < ncols; ++col)
669  {
670  out << "(" << row << "," << colnums[col] << ") " << values[col]
671  << std::endl;
672  }
673 
674  ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values);
675  AssertThrow(ierr == 0, ExcPETScError(ierr));
676  }
677 
678  AssertThrow(out.fail() == false, ExcIO());
679  }
680 
681 
682 
683  std::size_t
685  {
686  MatInfo info;
687  const PetscErrorCode ierr = MatGetInfo(matrix, MAT_LOCAL, &info);
688  AssertThrow(ierr == 0, ExcPETScError(ierr));
689 
690  return sizeof(*this) + static_cast<size_type>(info.memory);
691  }
692 
693 } // namespace PETScWrappers
694 
696 
697 #endif // DEAL_II_WITH_PETSC
void add(const size_type i, const size_type j, const PetscScalar value)
size_type row_length(const size_type row) const
std::size_t memory_consumption() const
VectorOperation::values last_action
void vmult(VectorBase &dst, const VectorBase &src) const
void clear_rows(const std::vector< size_type > &rows, const PetscScalar new_diag_value=0)
PetscScalar diag_element(const size_type i) const
const_iterator begin() const
MatrixBase & operator/=(const PetscScalar factor)
void mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
PetscBool is_symmetric(const double tolerance=1.e-12)
types::global_dof_index size_type
PetscReal frobenius_norm() const
MatrixBase & operator=(const MatrixBase &)=delete
virtual ~MatrixBase() override
const_iterator end() const
void Tvmult_add(VectorBase &dst, const VectorBase &src) const
void print(std::ostream &out, const bool alternative_output=false) const
PetscScalar el(const size_type i, const size_type j) const
PetscBool is_hermitian(const double tolerance=1.e-12)
PetscScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
void Tvmult(VectorBase &dst, const VectorBase &src) const
MatrixBase & operator*=(const PetscScalar factor)
PetscScalar residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const
void Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
void vmult_add(VectorBase &dst, const VectorBase &src) const
std::pair< size_type, size_type > local_range() const
virtual const MPI_Comm & get_mpi_communicator() const =0
void compress(const VectorOperation::values operation)
PetscScalar matrix_norm_square(const VectorBase &v) const
std::uint64_t n_nonzero_elements() const
void clear_row(const size_type row, const PetscScalar new_diag_value=0)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
@ matrix
Contents is actually a matrix.
static const char V
void perform_mmult(const MatrixBase &inputleft, const MatrixBase &inputright, MatrixBase &result, const VectorBase &V, const bool transpose_left)
PetscErrorCode destroy_matrix(Mat &matrix)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)