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