Reference documentation for deal.II version GIT c14369f203 2023-10-01 07: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.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2023 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 
687  size_type
688  local_domain_size() const;
689 
696  std::pair<size_type, size_type>
697  local_domain() const;
698 
702  MPI_Comm
704 
710  std::uint64_t
711  n_nonzero_elements() const;
712 
716  size_type
717  row_length(const size_type row) const;
718 
726  PetscReal
727  l1_norm() const;
728 
736  PetscReal
737  linfty_norm() const;
738 
743  PetscReal
744  frobenius_norm() const;
745 
746 
766  PetscScalar
767  matrix_norm_square(const VectorBase &v) const;
768 
769 
783  PetscScalar
784  matrix_scalar_product(const VectorBase &u, const VectorBase &v) const;
785 
790  PetscScalar
791  trace() const;
792 
796  MatrixBase &
797  operator*=(const PetscScalar factor);
798 
802  MatrixBase &
803  operator/=(const PetscScalar factor);
804 
805 
810  MatrixBase &
811  add(const PetscScalar factor, const MatrixBase &other);
812 
824  void
825  vmult(VectorBase &dst, const VectorBase &src) const;
826 
839  void
840  Tvmult(VectorBase &dst, const VectorBase &src) const;
841 
853  void
854  vmult_add(VectorBase &dst, const VectorBase &src) const;
855 
868  void
869  Tvmult_add(VectorBase &dst, const VectorBase &src) const;
870 
883  PetscScalar
884  residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const;
885 
892  begin() const;
893 
900  end() const;
901 
911  begin(const size_type r) const;
912 
922  end(const size_type r) const;
923 
931  operator Mat() const;
932 
938  Mat &
939  petsc_matrix();
940 
944  void
945  transpose();
946 
951  PetscBool
952  is_symmetric(const double tolerance = 1.e-12);
953 
959  PetscBool
960  is_hermitian(const double tolerance = 1.e-12);
961 
968  void
969  write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
970 
978  void
979  print(std::ostream &out, const bool alternative_output = false) const;
980 
984  std::size_t
985  memory_consumption() const;
986 
991  "You are attempting an operation on two vectors that "
992  "are the same object, but the operation requires that the "
993  "two objects are in fact different.");
994 
999  int,
1000  int,
1001  << "You tried to do a "
1002  << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
1003  << " operation but the matrix is currently in "
1004  << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
1005  << " mode. You first have to call 'compress()'.");
1006 
1007  protected:
1012  Mat matrix;
1013 
1018 
1024  void
1026 
1032  void
1034 
1045  void
1052  void
1054 
1070  void
1071  mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
1072 
1089  void
1090  Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
1091 
1092  private:
1107  mutable std::vector<PetscInt> column_indices;
1108 
1117  mutable std::vector<PetscScalar> column_values;
1118 
1119 
1120  // To allow calling protected prepare_add() and prepare_set().
1121  template <class>
1122  friend class ::BlockMatrixBase;
1123  };
1124 
1125 
1126 
1127 # ifndef DOXYGEN
1128  // ---------------------- inline and template functions ---------------------
1129 
1130 
1131  namespace MatrixIterators
1132  {
1134  const size_type row,
1135  const size_type index)
1136  : matrix(const_cast<MatrixBase *>(matrix))
1137  , a_row(row)
1138  , a_index(index)
1139  {
1140  visit_present_row();
1141  }
1142 
1143 
1144 
1147  {
1148  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1149  return a_row;
1150  }
1151 
1152 
1155  {
1156  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1157  return (*colnum_cache)[a_index];
1158  }
1159 
1160 
1163  {
1164  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1165  return a_index;
1166  }
1167 
1168 
1169  inline PetscScalar
1171  {
1172  Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1173  return (*value_cache)[a_index];
1174  }
1175 
1176 
1177  inline const_iterator::const_iterator(const MatrixBase *matrix,
1178  const size_type row,
1179  const size_type index)
1180  : accessor(matrix, row, index)
1181  {}
1182 
1183 
1184 
1185  inline const_iterator &
1187  {
1188  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
1189 
1190  ++accessor.a_index;
1191 
1192  // if at end of line: do one step, then cycle until we find a
1193  // row with a nonzero number of entries
1194  if (accessor.a_index >= accessor.colnum_cache->size())
1195  {
1196  accessor.a_index = 0;
1197  ++accessor.a_row;
1198 
1199  while ((accessor.a_row < accessor.matrix->m()) &&
1200  (accessor.a_row < accessor.matrix->local_range().second) &&
1201  (accessor.matrix->row_length(accessor.a_row) == 0))
1202  ++accessor.a_row;
1203 
1204  accessor.visit_present_row();
1205  }
1206  return *this;
1207  }
1208 
1209 
1210  inline const_iterator
1212  {
1213  const const_iterator old_state = *this;
1214  ++(*this);
1215  return old_state;
1216  }
1217 
1218 
1219  inline const const_iterator::Accessor &
1221  {
1222  return accessor;
1223  }
1224 
1225 
1226  inline const const_iterator::Accessor *
1227  const_iterator::operator->() const
1228  {
1229  return &accessor;
1230  }
1231 
1232 
1233  inline bool
1234  const_iterator::operator==(const const_iterator &other) const
1235  {
1236  return (accessor.a_row == other.accessor.a_row &&
1237  accessor.a_index == other.accessor.a_index);
1238  }
1239 
1240 
1241  inline bool
1242  const_iterator::operator!=(const const_iterator &other) const
1243  {
1244  return !(*this == other);
1245  }
1246 
1247 
1248  inline bool
1249  const_iterator::operator<(const const_iterator &other) const
1250  {
1251  return (accessor.row() < other.accessor.row() ||
1252  (accessor.row() == other.accessor.row() &&
1253  accessor.index() < other.accessor.index()));
1254  }
1255 
1256  } // namespace MatrixIterators
1257 
1258 
1259 
1260  // Inline the set() and add()
1261  // functions, since they will be
1262  // called frequently, and the
1263  // compiler can optimize away
1264  // some unnecessary loops when
1265  // the sizes are given at
1266  // compile time.
1267  inline void
1268  MatrixBase::set(const size_type i, const size_type j, const PetscScalar value)
1269  {
1270  AssertIsFinite(value);
1271 
1272  set(i, 1, &j, &value, false);
1273  }
1274 
1275 
1276 
1277  inline void
1278  MatrixBase::set(const std::vector<size_type> &indices,
1280  const bool elide_zero_values)
1281  {
1282  Assert(indices.size() == values.m(),
1283  ExcDimensionMismatch(indices.size(), values.m()));
1284  Assert(values.m() == values.n(), ExcNotQuadratic());
1285 
1286  for (size_type i = 0; i < indices.size(); ++i)
1287  set(indices[i],
1288  indices.size(),
1289  indices.data(),
1290  &values(i, 0),
1291  elide_zero_values);
1292  }
1293 
1294 
1295 
1296  inline void
1297  MatrixBase::set(const std::vector<size_type> &row_indices,
1298  const std::vector<size_type> &col_indices,
1300  const bool elide_zero_values)
1301  {
1302  Assert(row_indices.size() == values.m(),
1303  ExcDimensionMismatch(row_indices.size(), values.m()));
1304  Assert(col_indices.size() == values.n(),
1305  ExcDimensionMismatch(col_indices.size(), values.n()));
1306 
1307  for (size_type i = 0; i < row_indices.size(); ++i)
1308  set(row_indices[i],
1309  col_indices.size(),
1310  col_indices.data(),
1311  &values(i, 0),
1312  elide_zero_values);
1313  }
1314 
1315 
1316 
1317  inline void
1318  MatrixBase::set(const size_type row,
1319  const std::vector<size_type> &col_indices,
1320  const std::vector<PetscScalar> &values,
1321  const bool elide_zero_values)
1322  {
1323  Assert(col_indices.size() == values.size(),
1324  ExcDimensionMismatch(col_indices.size(), values.size()));
1325 
1326  set(row,
1327  col_indices.size(),
1328  col_indices.data(),
1329  values.data(),
1330  elide_zero_values);
1331  }
1332 
1333 
1334 
1335  inline void
1336  MatrixBase::set(const size_type row,
1337  const size_type n_cols,
1338  const size_type *col_indices,
1339  const PetscScalar *values,
1340  const bool elide_zero_values)
1341  {
1342  prepare_action(VectorOperation::insert);
1343 
1344  const PetscInt petsc_i = row;
1345  const PetscInt *col_index_ptr;
1346 
1347  const PetscScalar *col_value_ptr;
1348  int n_columns;
1349 
1350  // If we don't elide zeros, the pointers are already available...
1351  if (elide_zero_values == false)
1352  {
1353  col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
1354  col_value_ptr = values;
1355  n_columns = n_cols;
1356  }
1357  else
1358  {
1359  // Otherwise, extract nonzero values in each row and get the
1360  // respective index.
1361  if (column_indices.size() < n_cols)
1362  {
1363  column_indices.resize(n_cols);
1364  column_values.resize(n_cols);
1365  }
1366 
1367  n_columns = 0;
1368  for (size_type j = 0; j < n_cols; ++j)
1369  {
1370  const PetscScalar value = values[j];
1371  AssertIsFinite(value);
1372  if (value != PetscScalar())
1373  {
1374  column_indices[n_columns] = col_indices[j];
1375  column_values[n_columns] = value;
1376  n_columns++;
1377  }
1378  }
1379  AssertIndexRange(n_columns, n_cols + 1);
1380 
1381  col_index_ptr = column_indices.data();
1382  col_value_ptr = column_values.data();
1383  }
1384 
1385  const PetscErrorCode ierr = MatSetValues(matrix,
1386  1,
1387  &petsc_i,
1388  n_columns,
1389  col_index_ptr,
1390  col_value_ptr,
1391  INSERT_VALUES);
1392  AssertThrow(ierr == 0, ExcPETScError(ierr));
1393  }
1394 
1395 
1396 
1397  inline void
1398  MatrixBase::add(const size_type i, const size_type j, const PetscScalar value)
1399  {
1400  AssertIsFinite(value);
1401 
1402  if (value == PetscScalar())
1403  {
1404  // we have to check after using Insert/Add in any case to be
1405  // consistent with the MPI communication model, but we can save
1406  // some work if the addend is zero. However, these actions are done
1407  // in case we pass on to the other function.
1408  prepare_action(VectorOperation::add);
1409 
1410  return;
1411  }
1412  else
1413  add(i, 1, &j, &value, false);
1414  }
1415 
1416 
1417 
1418  inline void
1419  MatrixBase::add(const std::vector<size_type> &indices,
1421  const bool elide_zero_values)
1422  {
1423  Assert(indices.size() == values.m(),
1424  ExcDimensionMismatch(indices.size(), values.m()));
1425  Assert(values.m() == values.n(), ExcNotQuadratic());
1426 
1427  for (size_type i = 0; i < indices.size(); ++i)
1428  add(indices[i],
1429  indices.size(),
1430  indices.data(),
1431  &values(i, 0),
1432  elide_zero_values);
1433  }
1434 
1435 
1436 
1437  inline void
1438  MatrixBase::add(const std::vector<size_type> &row_indices,
1439  const std::vector<size_type> &col_indices,
1441  const bool elide_zero_values)
1442  {
1443  Assert(row_indices.size() == values.m(),
1444  ExcDimensionMismatch(row_indices.size(), values.m()));
1445  Assert(col_indices.size() == values.n(),
1446  ExcDimensionMismatch(col_indices.size(), values.n()));
1447 
1448  for (size_type i = 0; i < row_indices.size(); ++i)
1449  add(row_indices[i],
1450  col_indices.size(),
1451  col_indices.data(),
1452  &values(i, 0),
1453  elide_zero_values);
1454  }
1455 
1456 
1457 
1458  inline void
1459  MatrixBase::add(const size_type row,
1460  const std::vector<size_type> &col_indices,
1461  const std::vector<PetscScalar> &values,
1462  const bool elide_zero_values)
1463  {
1464  Assert(col_indices.size() == values.size(),
1465  ExcDimensionMismatch(col_indices.size(), values.size()));
1466 
1467  add(row,
1468  col_indices.size(),
1469  col_indices.data(),
1470  values.data(),
1471  elide_zero_values);
1472  }
1473 
1474 
1475 
1476  inline void
1477  MatrixBase::add(const size_type row,
1478  const size_type n_cols,
1479  const size_type *col_indices,
1480  const PetscScalar *values,
1481  const bool elide_zero_values,
1482  const bool /*col_indices_are_sorted*/)
1483  {
1484  (void)elide_zero_values;
1485 
1486  prepare_action(VectorOperation::add);
1487 
1488  const PetscInt petsc_i = row;
1489  const PetscInt *col_index_ptr;
1490 
1491  const PetscScalar *col_value_ptr;
1492  int n_columns;
1493 
1494  // If we don't elide zeros, the pointers are already available...
1495  if (elide_zero_values == false)
1496  {
1497  col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
1498  col_value_ptr = values;
1499  n_columns = n_cols;
1500  }
1501  else
1502  {
1503  // Otherwise, extract nonzero values in each row and get the
1504  // respective index.
1505  if (column_indices.size() < n_cols)
1506  {
1507  column_indices.resize(n_cols);
1508  column_values.resize(n_cols);
1509  }
1510 
1511  n_columns = 0;
1512  for (size_type j = 0; j < n_cols; ++j)
1513  {
1514  const PetscScalar value = values[j];
1515  AssertIsFinite(value);
1516  if (value != PetscScalar())
1517  {
1518  column_indices[n_columns] = col_indices[j];
1519  column_values[n_columns] = value;
1520  n_columns++;
1521  }
1522  }
1523  AssertIndexRange(n_columns, n_cols + 1);
1524 
1525  col_index_ptr = column_indices.data();
1526  col_value_ptr = column_values.data();
1527  }
1528 
1529  const PetscErrorCode ierr = MatSetValues(
1530  matrix, 1, &petsc_i, n_columns, col_index_ptr, col_value_ptr, ADD_VALUES);
1531  AssertThrow(ierr == 0, ExcPETScError(ierr));
1532  }
1533 
1534 
1535 
1536  inline PetscScalar
1537  MatrixBase::operator()(const size_type i, const size_type j) const
1538  {
1539  return el(i, j);
1540  }
1541 
1542 
1543 
1544  inline MatrixBase::const_iterator
1545  MatrixBase::begin() const
1546  {
1547  Assert(
1548  (in_local_range(0) && in_local_range(m() - 1)),
1549  ExcMessage(
1550  "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."));
1551 
1552  // find the first non-empty row in order to make sure that the returned
1553  // iterator points to something useful
1554  size_type first_nonempty_row = 0;
1555  while ((first_nonempty_row < m()) && (row_length(first_nonempty_row) == 0))
1556  ++first_nonempty_row;
1557 
1558  return const_iterator(this, first_nonempty_row, 0);
1559  }
1560 
1561 
1562  inline MatrixBase::const_iterator
1563  MatrixBase::end() const
1564  {
1565  Assert(
1566  (in_local_range(0) && in_local_range(m() - 1)),
1567  ExcMessage(
1568  "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."));
1569 
1570  return const_iterator(this, m(), 0);
1571  }
1572 
1573 
1574  inline MatrixBase::const_iterator
1575  MatrixBase::begin(const size_type r) const
1576  {
1577  Assert(in_local_range(r),
1578  ExcIndexRange(r, local_range().first, local_range().second));
1579 
1580  if (row_length(r) > 0)
1581  return const_iterator(this, r, 0);
1582  else
1583  return end(r);
1584  }
1585 
1586 
1587  inline MatrixBase::const_iterator
1588  MatrixBase::end(const size_type r) const
1589  {
1590  Assert(in_local_range(r),
1591  ExcIndexRange(r, local_range().first, local_range().second));
1592 
1593  // place the iterator on the first entry past this line, or at the
1594  // end of the matrix
1595  //
1596  // in the parallel case, we need to put it on the first entry of
1597  // the first row after the locally owned range. this of course
1598  // doesn't exist, but we can nevertheless create such an
1599  // iterator. we need to check whether 'i' is past the locally
1600  // owned range of rows first, before we ask for the length of the
1601  // row since the latter query leads to an exception in case we ask
1602  // for a row that is not locally owned
1603  for (size_type i = r + 1; i < m(); ++i)
1604  if (i == local_range().second || (row_length(i) > 0))
1605  return const_iterator(this, i, 0);
1606 
1607  // if there is no such line, then take the
1608  // end iterator of the matrix
1609  // we don't allow calling end() directly for distributed matrices so we need
1610  // to copy the code without the assertion.
1611  return {this, m(), 0};
1612  }
1613 
1614 
1615 
1616  inline bool
1617  MatrixBase::in_local_range(const size_type index) const
1618  {
1619  PetscInt begin, end;
1620 
1621  const PetscErrorCode ierr =
1622  MatGetOwnershipRange(static_cast<const Mat &>(matrix), &begin, &end);
1623  AssertThrow(ierr == 0, ExcPETScError(ierr));
1624 
1625  return ((index >= static_cast<size_type>(begin)) &&
1626  (index < static_cast<size_type>(end)));
1627  }
1628 
1629 
1630 
1631  inline void
1632  MatrixBase::prepare_action(const VectorOperation::values new_action)
1633  {
1634  if (last_action == VectorOperation::unknown)
1635  last_action = new_action;
1636 
1637  Assert(last_action == new_action, ExcWrongMode(last_action, new_action));
1638  }
1639 
1640 
1641 
1642  inline void
1643  MatrixBase::assert_is_compressed()
1644  {
1645  // compress() sets the last action to none, which allows us to check if
1646  // there are pending add/insert operations:
1647  AssertThrow(last_action == VectorOperation::unknown,
1648  ExcMessage("Error: missing compress() call."));
1649  }
1650 
1651 
1652 
1653  inline void
1654  MatrixBase::prepare_add()
1655  {
1656  prepare_action(VectorOperation::add);
1657  }
1658 
1659 
1660 
1661  inline void
1662  MatrixBase::prepare_set()
1663  {
1664  prepare_action(VectorOperation::insert);
1665  }
1666 
1667  inline MPI_Comm
1668  MatrixBase::get_mpi_communicator() const
1669  {
1670  return PetscObjectComm(reinterpret_cast<PetscObject>(matrix));
1671  }
1672 
1673 # endif // DOXYGEN
1674 } // namespace PETScWrappers
1675 
1676 
1678 
1679 
1680 #endif // DEAL_II_WITH_PETSC
1681 
1682 #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
MPI_Comm get_mpi_communicator() const
PetscScalar diag_element(const size_type i) const
MatrixBase(const MatrixBase &)=delete
size_type local_domain_size() 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
std::pair< size_type, size_type > local_domain() const
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)
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_v< T > &&std::is_floating_point_v< U >, 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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 2 > second
Definition: grid_out.cc:4615
Point< 2 > first
Definition: grid_out.cc:4614
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition: exceptions.h:467
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:1616
static ::ExceptionBase & ExcAccessToNonlocalRow(int arg1, int arg2, int arg3)
#define AssertIsFinite(number)
Definition: exceptions.h:1884
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:535
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:490
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:558
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
@ 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)