Reference documentation for deal.II version GIT b8135fa6eb 2023-03-25 15:55: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.h
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 
16 #ifndef dealii_petsc_matrix_base_h
17 #define dealii_petsc_matrix_base_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_WITH_PETSC
23 
25 
26 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/full_matrix.h>
31 
32 # include <petscmat.h>
33 
34 # include <cmath>
35 # include <memory>
36 # include <vector>
37 
39 
40 // Forward declarations
41 # ifndef DOXYGEN
42 template <typename Matrix>
43 class BlockMatrixBase;
44 # endif
45 
46 
47 namespace PETScWrappers
48 {
49  // forward declarations
50  class MatrixBase;
51 
52  namespace MatrixIterators
53  {
68  {
69  private:
73  class Accessor
74  {
75  public:
80 
86  const size_type row,
87  const size_type index);
88 
92  size_type
93  row() const;
94 
98  size_type
99  index() const;
100 
104  size_type
105  column() const;
106 
110  PetscScalar
111  value() const;
112 
121  int,
122  int,
123  int,
124  << "You tried to access row " << arg1
125  << " of a distributed matrix, but only rows " << arg2
126  << " through " << arg3
127  << " are stored locally and can be accessed.");
128 
129  private:
133  mutable MatrixBase *matrix;
134 
139 
144 
157  std::shared_ptr<const std::vector<size_type>> colnum_cache;
158 
162  std::shared_ptr<const std::vector<PetscScalar>> value_cache;
163 
169  void
171 
172  // Make enclosing class a friend.
173  friend class const_iterator;
174  };
175 
176  public:
181 
187  const size_type row,
188  const size_type index);
189 
195 
201 
205  const Accessor &
206  operator*() const;
207 
211  const Accessor *
212  operator->() const;
213 
218  bool
219  operator==(const const_iterator &) const;
223  bool
224  operator!=(const const_iterator &) const;
225 
231  bool
232  operator<(const const_iterator &) const;
233 
238  int,
239  int,
240  << "Attempt to access element " << arg2 << " of row "
241  << arg1 << " which doesn't have that many elements.");
242 
243  private:
248  };
249 
250  } // namespace MatrixIterators
251 
252 
284  class MatrixBase : public Subscriptor
285  {
286  public:
291 
296 
300  using value_type = PetscScalar;
301 
305  MatrixBase();
306 
325  explicit MatrixBase(const Mat &);
326 
332  MatrixBase(const MatrixBase &) = delete;
333 
339  MatrixBase &
340  operator=(const MatrixBase &) = delete;
341 
345  virtual ~MatrixBase() override;
346 
353  void
354  reinit(Mat A);
355 
365  MatrixBase &
366  operator=(const value_type d);
367 
372  void
373  clear();
374 
384  void
385  set(const size_type i, const size_type j, const PetscScalar value);
386 
406  void
407  set(const std::vector<size_type> & indices,
408  const FullMatrix<PetscScalar> &full_matrix,
409  const bool elide_zero_values = false);
410 
416  void
417  set(const std::vector<size_type> & row_indices,
418  const std::vector<size_type> & col_indices,
419  const FullMatrix<PetscScalar> &full_matrix,
420  const bool elide_zero_values = false);
421 
436  void
437  set(const size_type row,
438  const std::vector<size_type> & col_indices,
439  const std::vector<PetscScalar> &values,
440  const bool elide_zero_values = false);
441 
456  void
457  set(const size_type row,
458  const size_type n_cols,
459  const size_type * col_indices,
460  const PetscScalar *values,
461  const bool elide_zero_values = false);
462 
472  void
473  add(const size_type i, const size_type j, const PetscScalar value);
474 
494  void
495  add(const std::vector<size_type> & indices,
496  const FullMatrix<PetscScalar> &full_matrix,
497  const bool elide_zero_values = true);
498 
504  void
505  add(const std::vector<size_type> & row_indices,
506  const std::vector<size_type> & col_indices,
507  const FullMatrix<PetscScalar> &full_matrix,
508  const bool elide_zero_values = true);
509 
524  void
525  add(const size_type row,
526  const std::vector<size_type> & col_indices,
527  const std::vector<PetscScalar> &values,
528  const bool elide_zero_values = true);
529 
544  void
545  add(const size_type row,
546  const size_type n_cols,
547  const size_type * col_indices,
548  const PetscScalar *values,
549  const bool elide_zero_values = true,
550  const bool col_indices_are_sorted = false);
551 
568  void
569  clear_row(const size_type row, const PetscScalar new_diag_value = 0);
570 
579  void
580  clear_rows(const std::vector<size_type> &rows,
581  const PetscScalar new_diag_value = 0);
582 
586  void
587  clear_rows_columns(const std::vector<size_type> &row_and_column_indices,
588  const PetscScalar new_diag_value = 0);
589 
601  void
602  compress(const VectorOperation::values operation);
603 
615  PetscScalar
616  operator()(const size_type i, const size_type j) const;
617 
625  PetscScalar
626  el(const size_type i, const size_type j) const;
627 
637  PetscScalar
638  diag_element(const size_type i) const;
639 
643  size_type
644  m() const;
645 
649  size_type
650  n() const;
651 
660  size_type
661  local_size() const;
662 
671  std::pair<size_type, size_type>
672  local_range() const;
673 
678  bool
679  in_local_range(const size_type index) const;
680 
685  const MPI_Comm &
687 
693  std::uint64_t
694  n_nonzero_elements() const;
695 
699  size_type
700  row_length(const size_type row) const;
701 
709  PetscReal
710  l1_norm() const;
711 
719  PetscReal
720  linfty_norm() const;
721 
726  PetscReal
727  frobenius_norm() const;
728 
729 
749  PetscScalar
750  matrix_norm_square(const VectorBase &v) const;
751 
752 
766  PetscScalar
767  matrix_scalar_product(const VectorBase &u, const VectorBase &v) const;
768 
773  PetscScalar
774  trace() const;
775 
779  MatrixBase &
780  operator*=(const PetscScalar factor);
781 
785  MatrixBase &
786  operator/=(const PetscScalar factor);
787 
788 
793  MatrixBase &
794  add(const PetscScalar factor, const MatrixBase &other);
795 
807  void
808  vmult(VectorBase &dst, const VectorBase &src) const;
809 
822  void
823  Tvmult(VectorBase &dst, const VectorBase &src) const;
824 
836  void
837  vmult_add(VectorBase &dst, const VectorBase &src) const;
838 
851  void
852  Tvmult_add(VectorBase &dst, const VectorBase &src) const;
853 
866  PetscScalar
867  residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const;
868 
875  begin() const;
876 
883  end() const;
884 
894  begin(const size_type r) const;
895 
905  end(const size_type r) const;
906 
914  operator Mat() const;
915 
921  Mat &
922  petsc_matrix();
923 
927  void
928  transpose();
929 
934  PetscBool
935  is_symmetric(const double tolerance = 1.e-12);
936 
942  PetscBool
943  is_hermitian(const double tolerance = 1.e-12);
944 
951  void
952  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
953 
961  void
962  print(std::ostream &out, const bool alternative_output = false) const;
963 
967  std::size_t
968  memory_consumption() const;
969 
974  "You are attempting an operation on two vectors that "
975  "are the same object, but the operation requires that the "
976  "two objects are in fact different.");
977 
982  int,
983  int,
984  << "You tried to do a "
985  << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
986  << " operation but the matrix is currently in "
987  << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
988  << " mode. You first have to call 'compress()'.");
989 
990  protected:
995  Mat matrix;
996 
1001 
1007  void
1009 
1015  void
1017 
1028  void
1035  void
1037 
1053  void
1054  mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
1055 
1072  void
1073  Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
1074 
1075  private:
1090  mutable std::vector<PetscInt> column_indices;
1091 
1100  mutable std::vector<PetscScalar> column_values;
1101 
1102 
1103  // To allow calling protected prepare_add() and prepare_set().
1104  template <class>
1105  friend class ::BlockMatrixBase;
1106  };
1107 
1108 
1109 
1110 # ifndef DOXYGEN
1111  // ---------------------- inline and template functions ---------------------
1112 
1113 
1114  namespace MatrixIterators
1115  {
1117  const size_type row,
1118  const size_type index)
1119  : matrix(const_cast<MatrixBase *>(matrix))
1120  , a_row(row)
1121  , a_index(index)
1122  {
1123  visit_present_row();
1124  }
1125 
1126 
1127 
1130  {
1131  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1132  return a_row;
1133  }
1134 
1135 
1138  {
1139  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1140  return (*colnum_cache)[a_index];
1141  }
1142 
1143 
1146  {
1147  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1148  return a_index;
1149  }
1150 
1151 
1152  inline PetscScalar
1154  {
1155  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1156  return (*value_cache)[a_index];
1157  }
1158 
1159 
1160  inline const_iterator::const_iterator(const MatrixBase *matrix,
1161  const size_type row,
1162  const size_type index)
1163  : accessor(matrix, row, index)
1164  {}
1165 
1166 
1167 
1168  inline const_iterator &
1170  {
1171  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
1172 
1173  ++accessor.a_index;
1174 
1175  // if at end of line: do one step, then cycle until we find a
1176  // row with a nonzero number of entries
1177  if (accessor.a_index >= accessor.colnum_cache->size())
1178  {
1179  accessor.a_index = 0;
1180  ++accessor.a_row;
1181 
1182  while ((accessor.a_row < accessor.matrix->m()) &&
1183  (accessor.a_row < accessor.matrix->local_range().second) &&
1184  (accessor.matrix->row_length(accessor.a_row) == 0))
1185  ++accessor.a_row;
1186 
1187  accessor.visit_present_row();
1188  }
1189  return *this;
1190  }
1191 
1192 
1193  inline const_iterator
1195  {
1196  const const_iterator old_state = *this;
1197  ++(*this);
1198  return old_state;
1199  }
1200 
1201 
1202  inline const const_iterator::Accessor &
1204  {
1205  return accessor;
1206  }
1207 
1208 
1209  inline const const_iterator::Accessor *
1210  const_iterator::operator->() const
1211  {
1212  return &accessor;
1213  }
1214 
1215 
1216  inline bool
1217  const_iterator::operator==(const const_iterator &other) const
1218  {
1219  return (accessor.a_row == other.accessor.a_row &&
1220  accessor.a_index == other.accessor.a_index);
1221  }
1222 
1223 
1224  inline bool
1225  const_iterator::operator!=(const const_iterator &other) const
1226  {
1227  return !(*this == other);
1228  }
1229 
1230 
1231  inline bool
1232  const_iterator::operator<(const const_iterator &other) const
1233  {
1234  return (accessor.row() < other.accessor.row() ||
1235  (accessor.row() == other.accessor.row() &&
1236  accessor.index() < other.accessor.index()));
1237  }
1238 
1239  } // namespace MatrixIterators
1240 
1241 
1242 
1243  // Inline the set() and add()
1244  // functions, since they will be
1245  // called frequently, and the
1246  // compiler can optimize away
1247  // some unnecessary loops when
1248  // the sizes are given at
1249  // compile time.
1250  inline void
1251  MatrixBase::set(const size_type i, const size_type j, const PetscScalar value)
1252  {
1253  AssertIsFinite(value);
1254 
1255  set(i, 1, &j, &value, false);
1256  }
1257 
1258 
1259 
1260  inline void
1261  MatrixBase::set(const std::vector<size_type> & indices,
1263  const bool elide_zero_values)
1264  {
1265  Assert(indices.size() == values.m(),
1266  ExcDimensionMismatch(indices.size(), values.m()));
1267  Assert(values.m() == values.n(), ExcNotQuadratic());
1268 
1269  for (size_type i = 0; i < indices.size(); ++i)
1270  set(indices[i],
1271  indices.size(),
1272  indices.data(),
1273  &values(i, 0),
1274  elide_zero_values);
1275  }
1276 
1277 
1278 
1279  inline void
1280  MatrixBase::set(const std::vector<size_type> & row_indices,
1281  const std::vector<size_type> & col_indices,
1283  const bool elide_zero_values)
1284  {
1285  Assert(row_indices.size() == values.m(),
1286  ExcDimensionMismatch(row_indices.size(), values.m()));
1287  Assert(col_indices.size() == values.n(),
1288  ExcDimensionMismatch(col_indices.size(), values.n()));
1289 
1290  for (size_type i = 0; i < row_indices.size(); ++i)
1291  set(row_indices[i],
1292  col_indices.size(),
1293  col_indices.data(),
1294  &values(i, 0),
1295  elide_zero_values);
1296  }
1297 
1298 
1299 
1300  inline void
1301  MatrixBase::set(const size_type row,
1302  const std::vector<size_type> & col_indices,
1303  const std::vector<PetscScalar> &values,
1304  const bool elide_zero_values)
1305  {
1306  Assert(col_indices.size() == values.size(),
1307  ExcDimensionMismatch(col_indices.size(), values.size()));
1308 
1309  set(row,
1310  col_indices.size(),
1311  col_indices.data(),
1312  values.data(),
1313  elide_zero_values);
1314  }
1315 
1316 
1317 
1318  inline void
1319  MatrixBase::set(const size_type row,
1320  const size_type n_cols,
1321  const size_type * col_indices,
1322  const PetscScalar *values,
1323  const bool elide_zero_values)
1324  {
1325  prepare_action(VectorOperation::insert);
1326 
1327  const PetscInt petsc_i = row;
1328  PetscInt const *col_index_ptr;
1329 
1330  PetscScalar const *col_value_ptr;
1331  int n_columns;
1332 
1333  // If we don't elide zeros, the pointers are already available...
1334  if (elide_zero_values == false)
1335  {
1336  col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
1337  col_value_ptr = values;
1338  n_columns = n_cols;
1339  }
1340  else
1341  {
1342  // Otherwise, extract nonzero values in each row and get the
1343  // respective index.
1344  if (column_indices.size() < n_cols)
1345  {
1346  column_indices.resize(n_cols);
1347  column_values.resize(n_cols);
1348  }
1349 
1350  n_columns = 0;
1351  for (size_type j = 0; j < n_cols; ++j)
1352  {
1353  const PetscScalar value = values[j];
1354  AssertIsFinite(value);
1355  if (value != PetscScalar())
1356  {
1357  column_indices[n_columns] = col_indices[j];
1358  column_values[n_columns] = value;
1359  n_columns++;
1360  }
1361  }
1362  AssertIndexRange(n_columns, n_cols + 1);
1363 
1364  col_index_ptr = column_indices.data();
1365  col_value_ptr = column_values.data();
1366  }
1367 
1368  const PetscErrorCode ierr = MatSetValues(matrix,
1369  1,
1370  &petsc_i,
1371  n_columns,
1372  col_index_ptr,
1373  col_value_ptr,
1374  INSERT_VALUES);
1375  AssertThrow(ierr == 0, ExcPETScError(ierr));
1376  }
1377 
1378 
1379 
1380  inline void
1381  MatrixBase::add(const size_type i, const size_type j, const PetscScalar value)
1382  {
1383  AssertIsFinite(value);
1384 
1385  if (value == PetscScalar())
1386  {
1387  // we have to check after using Insert/Add in any case to be
1388  // consistent with the MPI communication model, but we can save
1389  // some work if the addend is zero. However, these actions are done
1390  // in case we pass on to the other function.
1391  prepare_action(VectorOperation::add);
1392 
1393  return;
1394  }
1395  else
1396  add(i, 1, &j, &value, false);
1397  }
1398 
1399 
1400 
1401  inline void
1402  MatrixBase::add(const std::vector<size_type> & indices,
1404  const bool elide_zero_values)
1405  {
1406  Assert(indices.size() == values.m(),
1407  ExcDimensionMismatch(indices.size(), values.m()));
1408  Assert(values.m() == values.n(), ExcNotQuadratic());
1409 
1410  for (size_type i = 0; i < indices.size(); ++i)
1411  add(indices[i],
1412  indices.size(),
1413  indices.data(),
1414  &values(i, 0),
1415  elide_zero_values);
1416  }
1417 
1418 
1419 
1420  inline void
1421  MatrixBase::add(const std::vector<size_type> & row_indices,
1422  const std::vector<size_type> & col_indices,
1424  const bool elide_zero_values)
1425  {
1426  Assert(row_indices.size() == values.m(),
1427  ExcDimensionMismatch(row_indices.size(), values.m()));
1428  Assert(col_indices.size() == values.n(),
1429  ExcDimensionMismatch(col_indices.size(), values.n()));
1430 
1431  for (size_type i = 0; i < row_indices.size(); ++i)
1432  add(row_indices[i],
1433  col_indices.size(),
1434  col_indices.data(),
1435  &values(i, 0),
1436  elide_zero_values);
1437  }
1438 
1439 
1440 
1441  inline void
1442  MatrixBase::add(const size_type row,
1443  const std::vector<size_type> & col_indices,
1444  const std::vector<PetscScalar> &values,
1445  const bool elide_zero_values)
1446  {
1447  Assert(col_indices.size() == values.size(),
1448  ExcDimensionMismatch(col_indices.size(), values.size()));
1449 
1450  add(row,
1451  col_indices.size(),
1452  col_indices.data(),
1453  values.data(),
1454  elide_zero_values);
1455  }
1456 
1457 
1458 
1459  inline void
1460  MatrixBase::add(const size_type row,
1461  const size_type n_cols,
1462  const size_type * col_indices,
1463  const PetscScalar *values,
1464  const bool elide_zero_values,
1465  const bool /*col_indices_are_sorted*/)
1466  {
1467  (void)elide_zero_values;
1468 
1469  prepare_action(VectorOperation::add);
1470 
1471  const PetscInt petsc_i = row;
1472  PetscInt const *col_index_ptr;
1473 
1474  PetscScalar const *col_value_ptr;
1475  int n_columns;
1476 
1477  // If we don't elide zeros, the pointers are already available...
1478  if (elide_zero_values == false)
1479  {
1480  col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
1481  col_value_ptr = values;
1482  n_columns = n_cols;
1483  }
1484  else
1485  {
1486  // Otherwise, extract nonzero values in each row and get the
1487  // respective index.
1488  if (column_indices.size() < n_cols)
1489  {
1490  column_indices.resize(n_cols);
1491  column_values.resize(n_cols);
1492  }
1493 
1494  n_columns = 0;
1495  for (size_type j = 0; j < n_cols; ++j)
1496  {
1497  const PetscScalar value = values[j];
1498  AssertIsFinite(value);
1499  if (value != PetscScalar())
1500  {
1501  column_indices[n_columns] = col_indices[j];
1502  column_values[n_columns] = value;
1503  n_columns++;
1504  }
1505  }
1506  AssertIndexRange(n_columns, n_cols + 1);
1507 
1508  col_index_ptr = column_indices.data();
1509  col_value_ptr = column_values.data();
1510  }
1511 
1512  const PetscErrorCode ierr = MatSetValues(
1513  matrix, 1, &petsc_i, n_columns, col_index_ptr, col_value_ptr, ADD_VALUES);
1514  AssertThrow(ierr == 0, ExcPETScError(ierr));
1515  }
1516 
1517 
1518 
1519  inline PetscScalar
1520  MatrixBase::operator()(const size_type i, const size_type j) const
1521  {
1522  return el(i, j);
1523  }
1524 
1525 
1526 
1527  inline MatrixBase::const_iterator
1528  MatrixBase::begin() const
1529  {
1530  Assert(
1531  (in_local_range(0) && in_local_range(m() - 1)),
1532  ExcMessage(
1533  "begin() and end() can only be called on a processor owning the entire matrix. If this is a distributed matrix, use begin(row) and end(row) instead."));
1534 
1535  // find the first non-empty row in order to make sure that the returned
1536  // iterator points to something useful
1537  size_type first_nonempty_row = 0;
1538  while ((first_nonempty_row < m()) && (row_length(first_nonempty_row) == 0))
1539  ++first_nonempty_row;
1540 
1541  return const_iterator(this, first_nonempty_row, 0);
1542  }
1543 
1544 
1545  inline MatrixBase::const_iterator
1546  MatrixBase::end() const
1547  {
1548  Assert(
1549  (in_local_range(0) && in_local_range(m() - 1)),
1550  ExcMessage(
1551  "begin() and end() can only be called on a processor owning the entire matrix. If this is a distributed matrix, use begin(row) and end(row) instead."));
1552 
1553  return const_iterator(this, m(), 0);
1554  }
1555 
1556 
1557  inline MatrixBase::const_iterator
1558  MatrixBase::begin(const size_type r) const
1559  {
1560  Assert(in_local_range(r),
1561  ExcIndexRange(r, local_range().first, local_range().second));
1562 
1563  if (row_length(r) > 0)
1564  return const_iterator(this, r, 0);
1565  else
1566  return end(r);
1567  }
1568 
1569 
1570  inline MatrixBase::const_iterator
1571  MatrixBase::end(const size_type r) const
1572  {
1573  Assert(in_local_range(r),
1574  ExcIndexRange(r, local_range().first, local_range().second));
1575 
1576  // place the iterator on the first entry past this line, or at the
1577  // end of the matrix
1578  //
1579  // in the parallel case, we need to put it on the first entry of
1580  // the first row after the locally owned range. this of course
1581  // doesn't exist, but we can nevertheless create such an
1582  // iterator. we need to check whether 'i' is past the locally
1583  // owned range of rows first, before we ask for the length of the
1584  // row since the latter query leads to an exception in case we ask
1585  // for a row that is not locally owned
1586  for (size_type i = r + 1; i < m(); ++i)
1587  if (i == local_range().second || (row_length(i) > 0))
1588  return const_iterator(this, i, 0);
1589 
1590  // if there is no such line, then take the
1591  // end iterator of the matrix
1592  // we don't allow calling end() directly for distributed matrices so we need
1593  // to copy the code without the assertion.
1594  return {this, m(), 0};
1595  }
1596 
1597 
1598 
1599  inline bool
1600  MatrixBase::in_local_range(const size_type index) const
1601  {
1602  PetscInt begin, end;
1603 
1604  const PetscErrorCode ierr =
1605  MatGetOwnershipRange(static_cast<const Mat &>(matrix), &begin, &end);
1606  AssertThrow(ierr == 0, ExcPETScError(ierr));
1607 
1608  return ((index >= static_cast<size_type>(begin)) &&
1609  (index < static_cast<size_type>(end)));
1610  }
1611 
1612 
1613 
1614  inline void
1615  MatrixBase::prepare_action(const VectorOperation::values new_action)
1616  {
1617  if (last_action == VectorOperation::unknown)
1618  last_action = new_action;
1619 
1620  Assert(last_action == new_action, ExcWrongMode(last_action, new_action));
1621  }
1622 
1623 
1624 
1625  inline void
1626  MatrixBase::assert_is_compressed()
1627  {
1628  // compress() sets the last action to none, which allows us to check if
1629  // there are pending add/insert operations:
1630  AssertThrow(last_action == VectorOperation::unknown,
1631  ExcMessage("Error: missing compress() call."));
1632  }
1633 
1634 
1635 
1636  inline void
1637  MatrixBase::prepare_add()
1638  {
1639  prepare_action(VectorOperation::add);
1640  }
1641 
1642 
1643 
1644  inline void
1645  MatrixBase::prepare_set()
1646  {
1647  prepare_action(VectorOperation::insert);
1648  }
1649 
1650  inline const MPI_Comm &
1651  MatrixBase::get_mpi_communicator() const
1652  {
1653  static MPI_Comm comm = PETSC_COMM_SELF;
1654  MPI_Comm pcomm = PetscObjectComm(reinterpret_cast<PetscObject>(matrix));
1655  if (pcomm != MPI_COMM_NULL)
1656  comm = pcomm;
1657  return comm;
1658  }
1659 
1660 # endif // DOXYGEN
1661 } // namespace PETScWrappers
1662 
1663 
1665 
1666 
1667 #endif // DEAL_II_WITH_PETSC
1668 
1669 #endif
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void add(const size_type i, const size_type j, const PetscScalar value)
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=true)
const_iterator end(const size_type r) const
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)
std::vector< PetscScalar > column_values
PetscScalar diag_element(const size_type i) const
MatrixBase(const MatrixBase &)=delete
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
const_iterator begin(const size_type r) const
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< PetscScalar > &values, const bool elide_zero_values=true)
PetscScalar el(const size_type i, const size_type j) const
PetscScalar operator()(const size_type i, const size_type j) const
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=false)
void prepare_action(const VectorOperation::values new_action)
void add(const std::vector< size_type > &indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=true)
bool in_local_range(const size_type index) const
PetscBool is_hermitian(const double tolerance=1.e-12)
PetscScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
std::vector< PetscInt > column_indices
void Tvmult(VectorBase &dst, const VectorBase &src) const
void clear_rows_columns(const std::vector< size_type > &row_and_column_indices, const PetscScalar new_diag_value=0)
MatrixBase & operator*=(const PetscScalar factor)
const MPI_Comm & get_mpi_communicator() const
PetscScalar residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const
void Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< PetscScalar > &values, const bool elide_zero_values=false)
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const PetscScalar *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
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
void compress(const VectorOperation::values operation)
PetscScalar matrix_norm_square(const VectorBase &v) const
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const PetscScalar *values, const bool elide_zero_values=false)
void set(const std::vector< size_type > &indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=false)
void set(const size_type i, const size_type j, const PetscScalar value)
std::uint64_t n_nonzero_elements() const
void clear_row(const size_type row, const PetscScalar new_diag_value=0)
std::shared_ptr< const std::vector< PetscScalar > > value_cache
std::shared_ptr< const std::vector< size_type > > colnum_cache
Accessor(const MatrixBase *matrix, const size_type row, const size_type index)
bool operator!=(const const_iterator &) const
const_iterator(const MatrixBase *matrix, const size_type row, const size_type index)
bool operator==(const const_iterator &) const
bool operator<(const const_iterator &) const
std::enable_if_t< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 2 > second
Definition: grid_out.cc:4616
Point< 2 > first
Definition: grid_out.cc:4615
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcInvalidIndexWithinRow(int arg1, int arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcBeyondEndOfMatrix()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcAccessToNonlocalRow(int arg1, int arg2, int arg3)
#define AssertIsFinite(number)
Definition: exceptions.h:1854
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:533
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:488
static ::ExceptionBase & ExcWrongMode(int arg1, int arg2)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:556
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
@ matrix
Contents is actually a matrix.
static const char A
static const char V
types::global_dof_index size_type
Definition: cuda_kernels.h:45
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)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int global_dof_index
Definition: types.h:82
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
const MPI_Comm & comm