Reference documentation for deal.II version GIT 3a483e2de5 2022-11-28 08:40: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 
432  MatrixBase &
433  MatrixBase::add(const PetscScalar factor, const MatrixBase &other)
434  {
435  const PetscErrorCode ierr =
436  MatAXPY(matrix, factor, other, DIFFERENT_NONZERO_PATTERN);
437  AssertThrow(ierr == 0, ExcPETScError(ierr));
438 
439  return *this;
440  }
441 
442 
443  void
444  MatrixBase::vmult(VectorBase &dst, const VectorBase &src) const
445  {
446  Assert(&src != &dst, ExcSourceEqualsDestination());
447 
448  const PetscErrorCode ierr = MatMult(matrix, src, dst);
449  AssertThrow(ierr == 0, ExcPETScError(ierr));
450  }
451 
452 
453 
454  void
455  MatrixBase::Tvmult(VectorBase &dst, const VectorBase &src) const
456  {
457  Assert(&src != &dst, ExcSourceEqualsDestination());
458 
459  const PetscErrorCode ierr = MatMultTranspose(matrix, src, dst);
460  AssertThrow(ierr == 0, ExcPETScError(ierr));
461  }
462 
463 
464 
465  void
467  {
468  Assert(&src != &dst, ExcSourceEqualsDestination());
469 
470  const PetscErrorCode ierr = MatMultAdd(matrix, src, dst, dst);
471  AssertThrow(ierr == 0, ExcPETScError(ierr));
472  }
473 
474 
475 
476  void
478  {
479  Assert(&src != &dst, ExcSourceEqualsDestination());
480 
481  const PetscErrorCode ierr = MatMultTransposeAdd(matrix, src, dst, dst);
482  AssertThrow(ierr == 0, ExcPETScError(ierr));
483  }
484 
485 
486  namespace internals
487  {
488  void
489  perform_mmult(const MatrixBase &inputleft,
490  const MatrixBase &inputright,
491  MatrixBase & result,
492  const VectorBase &V,
493  const bool transpose_left)
494  {
495  const bool use_vector = (V.size() == inputright.m() ? true : false);
496  if (transpose_left == false)
497  {
498  Assert(inputleft.n() == inputright.m(),
499  ExcDimensionMismatch(inputleft.n(), inputright.m()));
500  }
501  else
502  {
503  Assert(inputleft.m() == inputright.m(),
504  ExcDimensionMismatch(inputleft.m(), inputright.m()));
505  }
506 
507  result.clear();
508 
509  PetscErrorCode ierr;
510 
511  if (use_vector == false)
512  {
513  if (transpose_left)
514  {
515  ierr = MatTransposeMatMult(inputleft,
516  inputright,
517  MAT_INITIAL_MATRIX,
518  PETSC_DEFAULT,
519  &result.petsc_matrix());
520  AssertThrow(ierr == 0, ExcPETScError(ierr));
521  }
522  else
523  {
524  ierr = MatMatMult(inputleft,
525  inputright,
526  MAT_INITIAL_MATRIX,
527  PETSC_DEFAULT,
528  &result.petsc_matrix());
529  AssertThrow(ierr == 0, ExcPETScError(ierr));
530  }
531  }
532  else
533  {
534  Mat tmp;
535  ierr = MatDuplicate(inputleft, MAT_COPY_VALUES, &tmp);
536  AssertThrow(ierr == 0, ExcPETScError(ierr));
537  if (transpose_left)
538  {
539 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
540  ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp);
541 # else
542  ierr = MatTranspose(tmp, MAT_INPLACE_MATRIX, &tmp);
543 # endif
544  AssertThrow(ierr == 0, ExcPETScError(ierr));
545  }
546  ierr = MatDiagonalScale(tmp, nullptr, V);
547  AssertThrow(ierr == 0, ExcPETScError(ierr));
548  ierr = MatMatMult(tmp,
549  inputright,
550  MAT_INITIAL_MATRIX,
551  PETSC_DEFAULT,
552  &result.petsc_matrix());
553  AssertThrow(ierr == 0, ExcPETScError(ierr));
554  ierr = PETScWrappers::destroy_matrix(tmp);
555  AssertThrow(ierr == 0, ExcPETScError(ierr));
556  }
557  }
558  } // namespace internals
559 
560  void
562  const MatrixBase &B,
563  const VectorBase &V) const
564  {
565  internals::perform_mmult(*this, B, C, V, false);
566  }
567 
568  void
570  const MatrixBase &B,
571  const VectorBase &V) const
572  {
573  internals::perform_mmult(*this, B, C, V, true);
574  }
575 
576  PetscScalar
578  const VectorBase &x,
579  const VectorBase &b) const
580  {
581  // avoid the use of a temporary, and
582  // rather do one negation pass more than
583  // necessary
584  vmult(dst, x);
585  dst -= b;
586  dst *= -1;
587 
588  return dst.l2_norm();
589  }
590 
591 
592 
593  MatrixBase::operator Mat() const
594  {
595  return matrix;
596  }
597 
598  Mat &
600  {
601  return matrix;
602  }
603 
604  void
606  {
607 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
608  const PetscErrorCode ierr = MatTranspose(matrix, MAT_REUSE_MATRIX, &matrix);
609 # else
610  const PetscErrorCode ierr =
611  MatTranspose(matrix, MAT_INPLACE_MATRIX, &matrix);
612 # endif
613  AssertThrow(ierr == 0, ExcPETScError(ierr));
614  }
615 
616  PetscBool
617  MatrixBase::is_symmetric(const double tolerance)
618  {
619  PetscBool truth;
621  const PetscErrorCode ierr = MatIsSymmetric(matrix, tolerance, &truth);
622  AssertThrow(ierr == 0, ExcPETScError(ierr));
623  return truth;
624  }
625 
626  PetscBool
627  MatrixBase::is_hermitian(const double tolerance)
628  {
629  PetscBool truth;
630 
632  const PetscErrorCode ierr = MatIsHermitian(matrix, tolerance, &truth);
633  AssertThrow(ierr == 0, ExcPETScError(ierr));
634 
635  return truth;
636  }
637 
638  void
639  MatrixBase::write_ascii(const PetscViewerFormat format)
640  {
642 
643  // Set options
644  PetscErrorCode ierr =
645  PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, format);
646  AssertThrow(ierr == 0, ExcPETScError(ierr));
647 
648  // Write to screen
649  ierr = MatView(matrix, PETSC_VIEWER_STDOUT_WORLD);
650  AssertThrow(ierr == 0, ExcPETScError(ierr));
651  }
652 
653  void
654  MatrixBase::print(std::ostream &out, const bool /*alternative_output*/) const
655  {
656  std::pair<MatrixBase::size_type, MatrixBase::size_type> loc_range =
657  local_range();
658 
659  PetscInt ncols;
660  const PetscInt * colnums;
661  const PetscScalar *values;
662 
664  for (row = loc_range.first; row < loc_range.second; ++row)
665  {
666  PetscErrorCode ierr = MatGetRow(*this, row, &ncols, &colnums, &values);
667  AssertThrow(ierr == 0, ExcPETScError(ierr));
668 
669  for (PetscInt col = 0; col < ncols; ++col)
670  {
671  out << "(" << row << "," << colnums[col] << ") " << values[col]
672  << std::endl;
673  }
674 
675  ierr = MatRestoreRow(*this, row, &ncols, &colnums, &values);
676  AssertThrow(ierr == 0, ExcPETScError(ierr));
677  }
678 
679  AssertThrow(out.fail() == false, ExcIO());
680  }
681 
682 
683 
684  std::size_t
686  {
687  MatInfo info;
688  const PetscErrorCode ierr = MatGetInfo(matrix, MAT_LOCAL, &info);
689  AssertThrow(ierr == 0, ExcPETScError(ierr));
690 
691  return sizeof(*this) + static_cast<size_type>(info.memory);
692  }
693 
694 } // namespace PETScWrappers
695 
697 
698 #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:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1818
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcSourceEqualsDestination()
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 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)