deal.II version GIT relicensing-2017-g7ae82fad8b 2024-10-22 19:20:00+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\}}\)
Loading...
Searching...
No Matches
petsc_matrix_base.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_petsc_matrix_base_h
16#define dealii_petsc_matrix_base_h
17
18
19#include <deal.II/base/config.h>
20
21#ifdef DEAL_II_WITH_PETSC
22
24
30
31# include <boost/container/small_vector.hpp>
32
33# include <petscmat.h>
34
35# include <cmath>
36# include <memory>
37# include <optional>
38# include <vector>
39
41
42// Forward declarations
43# ifndef DOXYGEN
44template <typename Matrix>
45class BlockMatrixBase;
46# endif
47
48
49namespace PETScWrappers
50{
51 // forward declarations
52 class MatrixBase;
53
54 namespace MatrixIterators
55 {
70 {
71# ifdef __CUDA_ARCH__
72 // NVCC, at least until 12.6, fails to compile the
73 // implementations of the nested Accessor class if
74 // it is declared as private. Work around this by
75 // making it public.
76 public:
77# else
78 private:
79# endif
84 {
85 public:
90
96 const size_type row,
97 const size_type index);
98
103 row() const;
104
109 index() const;
110
115 column() const;
116
120 PetscScalar
121 value() const;
122
131 int,
132 int,
133 int,
134 << "You tried to access row " << arg1
135 << " of a distributed matrix, but only rows " << arg2
136 << " through " << arg3
137 << " are stored locally and can be accessed.");
138
139 private:
144
149
154
167 std::shared_ptr<const std::vector<size_type>> colnum_cache;
168
172 std::shared_ptr<const std::vector<PetscScalar>> value_cache;
173
179 void
181
182 // Make enclosing class a friend.
183 friend class const_iterator;
184 };
185
186 public:
191
197 const size_type row,
198 const size_type index);
199
205
211
215 const Accessor &
216 operator*() const;
217
221 const Accessor *
222 operator->() const;
223
228 bool
233 bool
235
241 bool
242 operator<(const const_iterator &) const;
243
248 int,
249 int,
250 << "Attempt to access element " << arg2 << " of row "
251 << arg1 << " which doesn't have that many elements.");
252
253 private:
258 };
259
260 } // namespace MatrixIterators
261
262
295 {
296 public:
301
306
310 using value_type = PetscScalar;
311
315 MatrixBase();
316
335 explicit MatrixBase(const Mat &);
336
342 MatrixBase(const MatrixBase &) = delete;
343
349 MatrixBase &
350 operator=(const MatrixBase &) = delete;
351
355 virtual ~MatrixBase() override;
356
363 void
364 reinit(Mat A);
365
375 MatrixBase &
376 operator=(const value_type d);
377
382 void
383 clear();
384
394 void
395 set(const size_type i, const size_type j, const PetscScalar value);
396
416 void
417 set(const std::vector<size_type> &indices,
418 const FullMatrix<PetscScalar> &full_matrix,
419 const bool elide_zero_values = false);
420
426 void
427 set(const std::vector<size_type> &row_indices,
428 const std::vector<size_type> &col_indices,
429 const FullMatrix<PetscScalar> &full_matrix,
430 const bool elide_zero_values = false);
431
446 void
447 set(const size_type row,
448 const std::vector<size_type> &col_indices,
449 const std::vector<PetscScalar> &values,
450 const bool elide_zero_values = false);
451
466 void
467 set(const size_type row,
468 const size_type n_cols,
469 const size_type *col_indices,
470 const PetscScalar *values,
471 const bool elide_zero_values = false);
472
482 void
483 add(const size_type i, const size_type j, const PetscScalar value);
484
504 void
505 add(const std::vector<size_type> &indices,
506 const FullMatrix<PetscScalar> &full_matrix,
507 const bool elide_zero_values = true);
508
514 void
515 add(const std::vector<size_type> &row_indices,
516 const std::vector<size_type> &col_indices,
517 const FullMatrix<PetscScalar> &full_matrix,
518 const bool elide_zero_values = true);
519
534 void
535 add(const size_type row,
536 const std::vector<size_type> &col_indices,
537 const std::vector<PetscScalar> &values,
538 const bool elide_zero_values = true);
539
554 void
555 add(const size_type row,
556 const size_type n_cols,
557 const size_type *col_indices,
558 const PetscScalar *values,
559 const bool elide_zero_values = true,
560 const bool col_indices_are_sorted = false);
561
578 void
579 clear_row(const size_type row, const PetscScalar new_diag_value = 0);
580
589 void
591 const PetscScalar new_diag_value = 0);
592
596 void
597 clear_rows_columns(const std::vector<size_type> &row_and_column_indices,
598 const PetscScalar new_diag_value = 0);
599
611 void
612 compress(const VectorOperation::values operation);
613
625 PetscScalar
626 operator()(const size_type i, const size_type j) const;
627
635 PetscScalar
636 el(const size_type i, const size_type j) const;
637
647 PetscScalar
648 diag_element(const size_type i) const;
649
654 m() const;
655
660 n() const;
661
671 local_size() const;
672
681 std::pair<size_type, size_type>
682 local_range() const;
683
688 bool
689 in_local_range(const size_type index) const;
690
698 local_domain_size() const;
699
706 std::pair<size_type, size_type>
707 local_domain() const;
708
714
720 std::uint64_t
721 n_nonzero_elements() const;
722
727 row_length(const size_type row) const;
728
736 PetscReal
737 l1_norm() const;
738
746 PetscReal
747 linfty_norm() const;
748
753 PetscReal
754 frobenius_norm() const;
755
756
776 PetscScalar
777 matrix_norm_square(const VectorBase &v) const;
778
779
793 PetscScalar
794 matrix_scalar_product(const VectorBase &u, const VectorBase &v) const;
795
800 PetscScalar
801 trace() const;
802
806 MatrixBase &
807 operator*=(const PetscScalar factor);
808
812 MatrixBase &
813 operator/=(const PetscScalar factor);
814
815
820 MatrixBase &
821 add(const PetscScalar factor, const MatrixBase &other);
822
834 void
835 vmult(VectorBase &dst, const VectorBase &src) const;
836
849 void
850 Tvmult(VectorBase &dst, const VectorBase &src) const;
851
863 void
864 vmult_add(VectorBase &dst, const VectorBase &src) const;
865
878 void
879 Tvmult_add(VectorBase &dst, const VectorBase &src) const;
880
893 PetscScalar
894 residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const;
895
902 begin() const;
903
910 end() const;
911
921 begin(const size_type r) const;
922
932 end(const size_type r) const;
933
941 operator Mat() const;
942
948 Mat &
949 petsc_matrix();
950
954 void
955 transpose();
956
961 PetscBool
962 is_symmetric(const double tolerance = 1.e-12);
963
969 PetscBool
970 is_hermitian(const double tolerance = 1.e-12);
971
978 void
979 write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
980
988 void
989 print(std::ostream &out, const bool alternative_output = false) const;
990
994 std::size_t
995 memory_consumption() const;
996
1001 "You are attempting an operation on two vectors that "
1002 "are the same object, but the operation requires that the "
1003 "two objects are in fact different.");
1004
1009 int,
1010 int,
1011 << "You tried to do a "
1012 << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
1013 << " operation but the matrix is currently in "
1014 << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
1015 << " mode. You first have to call 'compress()'.");
1016
1017 protected:
1023
1028
1034 void
1036
1042 void
1044
1055 void
1062 void
1064
1080 void
1081 mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
1082
1099 void
1100 Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
1101
1102 private:
1103 // To allow calling protected prepare_add() and prepare_set().
1104 template <class>
1105 friend class ::BlockMatrixBase;
1106
1107
1111 void
1113 const size_type row,
1114 const size_type n_cols,
1115 const size_type *col_indices,
1116 const PetscScalar *values,
1117 const bool elide_zero_values);
1118 };
1119
1120
1121
1122# ifndef DOXYGEN
1123 // ---------------------- inline and template functions ---------------------
1124
1125
1126 namespace MatrixIterators
1127 {
1129 const size_type row,
1130 const size_type index)
1131 : matrix(const_cast<MatrixBase *>(matrix))
1132 , a_row(row)
1133 , a_index(index)
1134 {
1135 visit_present_row();
1136 }
1137
1138
1139
1140 inline const_iterator::Accessor::size_type
1141 const_iterator::Accessor::row() const
1142 {
1143 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1144 return a_row;
1145 }
1146
1147
1148 inline const_iterator::Accessor::size_type
1149 const_iterator::Accessor::column() const
1150 {
1151 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1152 return (*colnum_cache)[a_index];
1153 }
1154
1155
1156 inline const_iterator::Accessor::size_type
1157 const_iterator::Accessor::index() const
1158 {
1159 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1160 return a_index;
1161 }
1162
1163
1164 inline PetscScalar
1165 const_iterator::Accessor::value() const
1166 {
1167 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1168 return (*value_cache)[a_index];
1169 }
1170
1171
1172 inline const_iterator::const_iterator(const MatrixBase *matrix,
1173 const size_type row,
1174 const size_type index)
1175 : accessor(matrix, row, index)
1176 {}
1177
1178
1179
1180 inline const_iterator &
1181 const_iterator::operator++()
1182 {
1183 Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
1184
1185 ++accessor.a_index;
1186
1187 // if at end of line: do one step, then cycle until we find a
1188 // row with a nonzero number of entries
1189 if (accessor.a_index >= accessor.colnum_cache->size())
1190 {
1191 accessor.a_index = 0;
1192 ++accessor.a_row;
1193
1194 while ((accessor.a_row < accessor.matrix->m()) &&
1195 (accessor.a_row < accessor.matrix->local_range().second) &&
1196 (accessor.matrix->row_length(accessor.a_row) == 0))
1197 ++accessor.a_row;
1198
1199 accessor.visit_present_row();
1200 }
1201 return *this;
1202 }
1203
1204
1205 inline const_iterator
1206 const_iterator::operator++(int)
1207 {
1208 const const_iterator old_state = *this;
1209 ++(*this);
1210 return old_state;
1211 }
1212
1213
1214 inline const const_iterator::Accessor &
1215 const_iterator::operator*() const
1216 {
1217 return accessor;
1218 }
1219
1220
1221 inline const const_iterator::Accessor *
1222 const_iterator::operator->() const
1223 {
1224 return &accessor;
1225 }
1226
1227
1228 inline bool
1229 const_iterator::operator==(const const_iterator &other) const
1230 {
1231 return (accessor.a_row == other.accessor.a_row &&
1232 accessor.a_index == other.accessor.a_index);
1233 }
1234
1235
1236 inline bool
1237 const_iterator::operator!=(const const_iterator &other) const
1238 {
1239 return !(*this == other);
1240 }
1241
1242
1243 inline bool
1244 const_iterator::operator<(const const_iterator &other) const
1245 {
1246 return (accessor.row() < other.accessor.row() ||
1247 (accessor.row() == other.accessor.row() &&
1248 accessor.index() < other.accessor.index()));
1249 }
1250
1251 } // namespace MatrixIterators
1252
1253
1254
1255 // Inline the set() and add()
1256 // functions, since they will be
1257 // called frequently, and the
1258 // compiler can optimize away
1259 // some unnecessary loops when
1260 // the sizes are given at
1261 // compile time.
1262 inline void
1263 MatrixBase::set(const size_type i, const size_type j, const PetscScalar value)
1264 {
1265 AssertIsFinite(value);
1266
1267 set(i, 1, &j, &value, false);
1268 }
1269
1270
1271
1272 inline void
1273 MatrixBase::set(const std::vector<size_type> &indices,
1274 const FullMatrix<PetscScalar> &values,
1275 const bool elide_zero_values)
1276 {
1277 Assert(indices.size() == values.m(),
1278 ExcDimensionMismatch(indices.size(), values.m()));
1279 Assert(values.m() == values.n(), ExcNotQuadratic());
1280
1281 for (size_type i = 0; i < indices.size(); ++i)
1282 set(indices[i],
1283 indices.size(),
1284 indices.data(),
1285 &values(i, 0),
1286 elide_zero_values);
1287 }
1288
1289
1290
1291 inline void
1292 MatrixBase::set(const std::vector<size_type> &row_indices,
1293 const std::vector<size_type> &col_indices,
1294 const FullMatrix<PetscScalar> &values,
1295 const bool elide_zero_values)
1296 {
1297 Assert(row_indices.size() == values.m(),
1298 ExcDimensionMismatch(row_indices.size(), values.m()));
1299 Assert(col_indices.size() == values.n(),
1300 ExcDimensionMismatch(col_indices.size(), values.n()));
1301
1302 for (size_type i = 0; i < row_indices.size(); ++i)
1303 set(row_indices[i],
1304 col_indices.size(),
1305 col_indices.data(),
1306 &values(i, 0),
1307 elide_zero_values);
1308 }
1309
1310
1311
1312 inline void
1313 MatrixBase::set(const size_type row,
1314 const std::vector<size_type> &col_indices,
1315 const std::vector<PetscScalar> &values,
1316 const bool elide_zero_values)
1317 {
1318 Assert(col_indices.size() == values.size(),
1319 ExcDimensionMismatch(col_indices.size(), values.size()));
1320
1321 set(row,
1322 col_indices.size(),
1323 col_indices.data(),
1324 values.data(),
1325 elide_zero_values);
1326 }
1327
1328
1329
1330 inline void
1331 MatrixBase::set(const size_type row,
1332 const size_type n_cols,
1333 const size_type *col_indices,
1334 const PetscScalar *values,
1335 const bool elide_zero_values)
1336 {
1337 add_or_set(VectorOperation::insert,
1338 row,
1339 n_cols,
1340 col_indices,
1341 values,
1342 elide_zero_values);
1343 }
1344
1345
1346
1347 inline void
1348 MatrixBase::add(const size_type i, const size_type j, const PetscScalar value)
1349 {
1350 AssertIsFinite(value);
1351
1352 if (value == PetscScalar())
1353 {
1354 // we have to check after using Insert/Add in any case to be
1355 // consistent with the MPI communication model, but we can save
1356 // some work if the addend is zero. However, these actions are done
1357 // in case we pass on to the other function.
1358 prepare_action(VectorOperation::add);
1359
1360 return;
1361 }
1362 else
1363 add(i, 1, &j, &value, false);
1364 }
1365
1366
1367
1368 inline void
1369 MatrixBase::add(const std::vector<size_type> &indices,
1370 const FullMatrix<PetscScalar> &values,
1371 const bool elide_zero_values)
1372 {
1373 Assert(indices.size() == values.m(),
1374 ExcDimensionMismatch(indices.size(), values.m()));
1375 Assert(values.m() == values.n(), ExcNotQuadratic());
1376
1377 for (size_type i = 0; i < indices.size(); ++i)
1378 add(indices[i],
1379 indices.size(),
1380 indices.data(),
1381 &values(i, 0),
1382 elide_zero_values);
1383 }
1384
1385
1386
1387 inline void
1388 MatrixBase::add(const std::vector<size_type> &row_indices,
1389 const std::vector<size_type> &col_indices,
1390 const FullMatrix<PetscScalar> &values,
1391 const bool elide_zero_values)
1392 {
1393 Assert(row_indices.size() == values.m(),
1394 ExcDimensionMismatch(row_indices.size(), values.m()));
1395 Assert(col_indices.size() == values.n(),
1396 ExcDimensionMismatch(col_indices.size(), values.n()));
1397
1398 for (size_type i = 0; i < row_indices.size(); ++i)
1399 add(row_indices[i],
1400 col_indices.size(),
1401 col_indices.data(),
1402 &values(i, 0),
1403 elide_zero_values);
1404 }
1405
1406
1407
1408 inline void
1409 MatrixBase::add(const size_type row,
1410 const std::vector<size_type> &col_indices,
1411 const std::vector<PetscScalar> &values,
1412 const bool elide_zero_values)
1413 {
1414 Assert(col_indices.size() == values.size(),
1415 ExcDimensionMismatch(col_indices.size(), values.size()));
1416
1417 add(row,
1418 col_indices.size(),
1419 col_indices.data(),
1420 values.data(),
1421 elide_zero_values);
1422 }
1423
1424
1425
1426 inline void
1427 MatrixBase::add(const size_type row,
1428 const size_type n_cols,
1429 const size_type *col_indices,
1430 const PetscScalar *values,
1431 const bool elide_zero_values,
1432 const bool /*col_indices_are_sorted*/)
1433 {
1434 add_or_set(VectorOperation::add,
1435 row,
1436 n_cols,
1437 col_indices,
1438 values,
1439 elide_zero_values);
1440 }
1441
1442
1443
1444 inline void
1445 MatrixBase::add_or_set(const VectorOperation::values &operation,
1446 const size_type row,
1447 const size_type n_cols,
1448 const size_type *col_indices,
1449 const PetscScalar *values,
1450 const bool elide_zero_values)
1451 {
1452 prepare_action(operation);
1453
1454 const auto petsc_row = static_cast<PetscInt>(row);
1455 AssertIntegerConversion(petsc_row, row);
1456
1457 // Use 100 entries so that we can store a row of a matrix constructed with
1458 // FESystem<3>(FE_Q<3>(2), 3, FE_Q<3>(1), 1) (i.e., a common Stokes FE,
1459 // which has 89 DoFs) with room to spare without allocating memory in the
1460 // free store.
1461 //
1462 // Setting up small_vectors isn't free so only set up column_values if we
1463 // actually use it.
1464 boost::container::small_vector<PetscInt, 100> column_indices;
1465 std::optional<boost::container::small_vector<PetscScalar, 100>>
1466 column_values;
1467
1468 const PetscScalar *values_ptr = nullptr;
1469 if (elide_zero_values == false)
1470 {
1471 column_indices.resize(n_cols);
1472
1473 for (size_type j = 0; j < n_cols; ++j)
1474 {
1475 AssertIsFinite(values[j]);
1476 column_indices[j] = static_cast<PetscInt>(col_indices[j]);
1477 AssertIntegerConversion(column_indices[j], col_indices[j]);
1478 }
1479
1480 values_ptr = values;
1481 }
1482 else
1483 {
1484 // Otherwise, extract nonzero values in each row and get the
1485 // respective index.
1486 column_values.emplace();
1487 for (size_type j = 0; j < n_cols; ++j)
1488 {
1489 const PetscScalar value = values[j];
1490 AssertIsFinite(value);
1491 if (value != PetscScalar())
1492 {
1493 column_indices.push_back(static_cast<PetscInt>(col_indices[j]));
1494 AssertIntegerConversion(column_indices.back(), col_indices[j]);
1495 column_values->push_back(value);
1496 }
1497 }
1498 values_ptr = column_values->data();
1499 }
1500
1501 const auto petsc_n_columns = static_cast<PetscInt>(column_indices.size());
1502 AssertIntegerConversion(petsc_n_columns, column_indices.size());
1503
1504 Assert(operation == VectorOperation::insert ||
1505 operation == VectorOperation::add,
1507 const PetscErrorCode ierr =
1508 MatSetValues(matrix,
1509 1,
1510 &petsc_row,
1511 petsc_n_columns,
1512 column_indices.data(),
1513 values_ptr,
1514 operation == VectorOperation::insert ? INSERT_VALUES :
1515 ADD_VALUES);
1516 AssertThrow(ierr == 0, ExcPETScError(ierr));
1517 }
1518
1519
1520
1521 inline PetscScalar
1522 MatrixBase::operator()(const size_type i, const size_type j) const
1523 {
1524 return el(i, j);
1525 }
1526
1527
1528
1529 inline MatrixBase::const_iterator
1530 MatrixBase::begin() const
1531 {
1532 Assert(
1533 (in_local_range(0) && in_local_range(m() - 1)),
1534 ExcMessage(
1535 "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."));
1536
1537 // find the first non-empty row in order to make sure that the returned
1538 // iterator points to something useful
1539 size_type first_nonempty_row = 0;
1540 while ((first_nonempty_row < m()) && (row_length(first_nonempty_row) == 0))
1541 ++first_nonempty_row;
1542
1543 return const_iterator(this, first_nonempty_row, 0);
1544 }
1545
1546
1547 inline MatrixBase::const_iterator
1548 MatrixBase::end() const
1549 {
1550 Assert(
1551 (in_local_range(0) && in_local_range(m() - 1)),
1552 ExcMessage(
1553 "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."));
1554
1555 return const_iterator(this, m(), 0);
1556 }
1557
1558
1559 inline MatrixBase::const_iterator
1560 MatrixBase::begin(const size_type r) const
1561 {
1562 Assert(in_local_range(r),
1563 ExcIndexRange(r, local_range().first, local_range().second));
1564
1565 if (row_length(r) > 0)
1566 return const_iterator(this, r, 0);
1567 else
1568 return end(r);
1569 }
1570
1571
1572 inline MatrixBase::const_iterator
1573 MatrixBase::end(const size_type r) const
1574 {
1575 Assert(in_local_range(r),
1576 ExcIndexRange(r, local_range().first, local_range().second));
1577
1578 // place the iterator on the first entry past this line, or at the
1579 // end of the matrix
1580 //
1581 // in the parallel case, we need to put it on the first entry of
1582 // the first row after the locally owned range. this of course
1583 // doesn't exist, but we can nevertheless create such an
1584 // iterator. we need to check whether 'i' is past the locally
1585 // owned range of rows first, before we ask for the length of the
1586 // row since the latter query leads to an exception in case we ask
1587 // for a row that is not locally owned
1588 for (size_type i = r + 1; i < m(); ++i)
1589 if (i == local_range().second || (row_length(i) > 0))
1590 return const_iterator(this, i, 0);
1591
1592 // if there is no such line, then take the
1593 // end iterator of the matrix
1594 // we don't allow calling end() directly for distributed matrices so we need
1595 // to copy the code without the assertion.
1596 return {this, m(), 0};
1597 }
1598
1599
1600
1601 inline bool
1602 MatrixBase::in_local_range(const size_type index) const
1603 {
1604 PetscInt petsc_begin, petsc_end;
1605
1606 const PetscErrorCode ierr =
1607 MatGetOwnershipRange(static_cast<const Mat &>(matrix),
1608 &petsc_begin,
1609 &petsc_end);
1610 AssertThrow(ierr == 0, ExcPETScError(ierr));
1611
1612 const auto begin = static_cast<size_type>(petsc_begin);
1613 AssertIntegerConversion(begin, petsc_begin);
1614 const auto end = static_cast<size_type>(petsc_end);
1615 AssertIntegerConversion(end, petsc_end);
1616
1617 return ((index >= begin) && (index < end));
1618 }
1619
1620
1621
1622 inline void
1623 MatrixBase::prepare_action(const VectorOperation::values new_action)
1624 {
1625 if (last_action == VectorOperation::unknown)
1626 last_action = new_action;
1627
1628 Assert(last_action == new_action, ExcWrongMode(last_action, new_action));
1629 }
1630
1631
1632
1633 inline void
1634 MatrixBase::assert_is_compressed()
1635 {
1636 // compress() sets the last action to none, which allows us to check if
1637 // there are pending add/insert operations:
1639 ExcMessage("Error: missing compress() call."));
1640 }
1641
1642
1643
1644 inline void
1645 MatrixBase::prepare_add()
1646 {
1647 prepare_action(VectorOperation::add);
1648 }
1649
1650
1651
1652 inline void
1653 MatrixBase::prepare_set()
1654 {
1655 prepare_action(VectorOperation::insert);
1656 }
1657
1658 inline MPI_Comm
1659 MatrixBase::get_mpi_communicator() const
1660 {
1661 return PetscObjectComm(reinterpret_cast<PetscObject>(matrix));
1662 }
1663
1664# endif // DOXYGEN
1665} // namespace PETScWrappers
1666
1667
1669
1670
1671#endif // DEAL_II_WITH_PETSC
1672
1673#endif
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
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)
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_or_set(const VectorOperation::values &operation, const size_type row, const size_type n_cols, const size_type *col_indices, const PetscScalar *values, const bool elide_zero_values)
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
MatrixBase & operator=(const MatrixBase &)=delete
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
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 clear_rows(const ArrayView< const size_type > &rows, const PetscScalar new_diag_value=0)
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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#define DeclException0(Exception0)
Definition exceptions.h:466
static ::ExceptionBase & ExcAccessToNonlocalRow(int arg1, int arg2, int arg3)
#define Assert(cond, exc)
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcInvalidIndexWithinRow(int arg1, int arg2)
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:534
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:489
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcWrongMode(int arg1, int arg2)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition exceptions.h:557
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define AssertIntegerConversion(index1, index2)
@ matrix
Contents is actually a matrix.
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int global_dof_index
Definition types.h:81