Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
precondition.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 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 #ifndef dealii_precondition_h
17 #define dealii_precondition_h
18 
19 // This file contains simple preconditioners.
20 
21 #include <deal.II/base/config.h>
22 
23 #include <deal.II/base/memory_space.h>
24 #include <deal.II/base/parallel.h>
25 #include <deal.II/base/smartpointer.h>
26 #include <deal.II/base/template_constraints.h>
27 #include <deal.II/base/thread_management.h>
28 #include <deal.II/base/utilities.h>
29 
30 #include <deal.II/lac/diagonal_matrix.h>
31 #include <deal.II/lac/solver_cg.h>
32 #include <deal.II/lac/vector_memory.h>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 // forward declarations
37 
38 template <typename number>
39 class Vector;
40 template <typename number>
41 class SparseMatrix;
42 namespace LinearAlgebra
43 {
44  namespace distributed
45  {
46  template <typename, typename>
47  class Vector;
48  } // namespace distributed
49 } // namespace LinearAlgebra
50 
51 
52 
79 {
80 public:
85 
91  {
95  AdditionalData() = default;
96  };
97 
102 
107  template <typename MatrixType>
108  void
109  initialize(const MatrixType & matrix,
110  const AdditionalData &additional_data = AdditionalData());
111 
115  template <class VectorType>
116  void
117  vmult(VectorType &, const VectorType &) const;
118 
123  template <class VectorType>
124  void
125  Tvmult(VectorType &, const VectorType &) const;
126 
130  template <class VectorType>
131  void
132  vmult_add(VectorType &, const VectorType &) const;
133 
138  template <class VectorType>
139  void
140  Tvmult_add(VectorType &, const VectorType &) const;
141 
146  void
148  {}
149 
157  size_type
158  m() const;
159 
167  size_type
168  n() const;
169 
170 private:
175 
180 };
181 
182 
183 
198 {
199 public:
204 
209  {
210  public:
215  AdditionalData(const double relaxation = 1.);
216 
220  double relaxation;
221  };
222 
228 
232  void
233  initialize(const AdditionalData &parameters);
234 
240  template <typename MatrixType>
241  void
242  initialize(const MatrixType &matrix, const AdditionalData &parameters);
243 
247  template <class VectorType>
248  void
249  vmult(VectorType &, const VectorType &) const;
250 
255  template <class VectorType>
256  void
257  Tvmult(VectorType &, const VectorType &) const;
261  template <class VectorType>
262  void
263  vmult_add(VectorType &, const VectorType &) const;
264 
269  template <class VectorType>
270  void
271  Tvmult_add(VectorType &, const VectorType &) const;
272 
277  void
279  {}
280 
288  size_type
289  m() const;
290 
298  size_type
299  n() const;
300 
301 private:
305  double relaxation;
306 
311 
316 };
317 
318 
319 
361 template <typename MatrixType = SparseMatrix<double>,
362  class VectorType = Vector<double>>
364 {
365 public:
369  using function_ptr = void (MatrixType::*)(VectorType &,
370  const VectorType &) const;
371 
377  PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
378 
383  void
384  vmult(VectorType &dst, const VectorType &src) const;
385 
386 private:
390  const MatrixType &matrix;
391 
396 };
397 
398 
399 
408 template <typename MatrixType = SparseMatrix<double>>
410 {
411 public:
415  using size_type = typename MatrixType::size_type;
416 
421  {
422  public:
426  AdditionalData(const double relaxation = 1.);
427 
431  double relaxation;
432  };
433 
439  void
440  initialize(const MatrixType & A,
441  const AdditionalData &parameters = AdditionalData());
442 
446  void
447  clear();
448 
453  size_type
454  m() const;
455 
460  size_type
461  n() const;
462 
463 protected:
468 
472  double relaxation;
473 };
474 
475 
476 
505 template <typename MatrixType = SparseMatrix<double>>
506 class PreconditionJacobi : public PreconditionRelaxation<MatrixType>
507 {
508 public:
512  using AdditionalData =
514 
518  template <class VectorType>
519  void
520  vmult(VectorType &, const VectorType &) const;
521 
526  template <class VectorType>
527  void
528  Tvmult(VectorType &, const VectorType &) const;
529 
533  template <class VectorType>
534  void
535  step(VectorType &x, const VectorType &rhs) const;
536 
540  template <class VectorType>
541  void
542  Tstep(VectorType &x, const VectorType &rhs) const;
543 };
544 
545 
593 template <typename MatrixType = SparseMatrix<double>>
594 class PreconditionSOR : public PreconditionRelaxation<MatrixType>
595 {
596 public:
600  using AdditionalData =
602 
606  template <class VectorType>
607  void
608  vmult(VectorType &, const VectorType &) const;
609 
613  template <class VectorType>
614  void
615  Tvmult(VectorType &, const VectorType &) const;
616 
620  template <class VectorType>
621  void
622  step(VectorType &x, const VectorType &rhs) const;
623 
627  template <class VectorType>
628  void
629  Tstep(VectorType &x, const VectorType &rhs) const;
630 };
631 
632 
633 
662 template <typename MatrixType = SparseMatrix<double>>
663 class PreconditionSSOR : public PreconditionRelaxation<MatrixType>
664 {
665 public:
669  using AdditionalData =
671 
675  using size_type = typename MatrixType::size_type;
676 
681 
682 
688  void
689  initialize(const MatrixType & A,
690  const typename BaseClass::AdditionalData &parameters =
691  typename BaseClass::AdditionalData());
692 
696  template <class VectorType>
697  void
698  vmult(VectorType &, const VectorType &) const;
699 
704  template <class VectorType>
705  void
706  Tvmult(VectorType &, const VectorType &) const;
707 
708 
712  template <class VectorType>
713  void
714  step(VectorType &x, const VectorType &rhs) const;
715 
719  template <class VectorType>
720  void
721  Tstep(VectorType &x, const VectorType &rhs) const;
722 
723 private:
728  std::vector<std::size_t> pos_right_of_diagonal;
729 };
730 
731 
764 template <typename MatrixType = SparseMatrix<double>>
765 class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
766 {
767 public:
771  using size_type = typename MatrixType::size_type;
772 
777  {
778  public:
790  const std::vector<size_type> &permutation,
791  const std::vector<size_type> &inverse_permutation,
793  &parameters =
795 
799  const std::vector<size_type> &permutation;
803  const std::vector<size_type> &inverse_permutation;
808  };
809 
821  void
822  initialize(const MatrixType & A,
823  const std::vector<size_type> &permutation,
824  const std::vector<size_type> &inverse_permutation,
826  &parameters =
828 
839  void
840  initialize(const MatrixType &A, const AdditionalData &additional_data);
841 
845  template <class VectorType>
846  void
847  vmult(VectorType &, const VectorType &) const;
848 
852  template <class VectorType>
853  void
854  Tvmult(VectorType &, const VectorType &) const;
855 
856 private:
860  const std::vector<size_type> *permutation;
864  const std::vector<size_type> *inverse_permutation;
865 };
866 
867 
868 
1003 template <typename MatrixType = SparseMatrix<double>,
1004  typename VectorType = Vector<double>,
1005  typename PreconditionerType = DiagonalMatrix<VectorType>>
1007 {
1008 public:
1013 
1014  // avoid warning about use of deprecated variables
1015 
1021  {
1025  AdditionalData(const unsigned int degree = 1,
1026  const double smoothing_range = 0.,
1027  const bool nonzero_starting = false,
1028  const unsigned int eig_cg_n_iterations = 8,
1029  const double eig_cg_residual = 1e-2,
1030  const double max_eigenvalue = 1);
1031 
1043  unsigned int degree;
1044 
1057 
1070  bool nonzero_starting DEAL_II_DEPRECATED;
1071 
1078  unsigned int eig_cg_n_iterations;
1079 
1085 
1092 
1098  VectorType matrix_diagonal_inverse DEAL_II_DEPRECATED;
1099 
1103  std::shared_ptr<PreconditionerType> preconditioner;
1104  };
1105 
1106 
1108 
1120  void
1121  initialize(const MatrixType & matrix,
1122  const AdditionalData &additional_data = AdditionalData());
1123 
1128  void
1129  vmult(VectorType &dst, const VectorType &src) const;
1130 
1135  void
1136  Tvmult(VectorType &dst, const VectorType &src) const;
1137 
1141  void
1142  step(VectorType &dst, const VectorType &src) const;
1143 
1147  void
1148  Tstep(VectorType &dst, const VectorType &src) const;
1149 
1153  void
1154  clear();
1155 
1160  size_type
1161  m() const;
1162 
1167  size_type
1168  n() const;
1169 
1170 private:
1174  SmartPointer<
1175  const MatrixType,
1178 
1182  mutable VectorType solution_old;
1183 
1187  mutable VectorType temp_vector1;
1188 
1192  mutable VectorType temp_vector2;
1193 
1199 
1203  double theta;
1204 
1209  double delta;
1210 
1216 
1222 
1229  void
1230  estimate_eigenvalues(const VectorType &src) const;
1231 };
1232 
1233 
1234 
1236 /* ---------------------------------- Inline functions ------------------- */
1237 
1238 #ifndef DOXYGEN
1239 
1241  : n_rows(0)
1242  , n_columns(0)
1243 {}
1244 
1245 template <typename MatrixType>
1246 inline void
1247 PreconditionIdentity::initialize(const MatrixType &matrix,
1249 {
1250  n_rows = matrix.m();
1251  n_columns = matrix.n();
1252 }
1253 
1254 
1255 template <class VectorType>
1256 inline void
1257 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1258 {
1259  dst = src;
1260 }
1261 
1262 
1263 
1264 template <class VectorType>
1265 inline void
1266 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1267 {
1268  dst = src;
1269 }
1270 
1271 template <class VectorType>
1272 inline void
1273 PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
1274 {
1275  dst += src;
1276 }
1277 
1278 
1279 
1280 template <class VectorType>
1281 inline void
1282 PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
1283 {
1284  dst += src;
1285 }
1286 
1289 {
1290  Assert(n_rows != 0, ExcNotInitialized());
1291  return n_rows;
1292 }
1293 
1296 {
1298  return n_columns;
1299 }
1300 
1301 //---------------------------------------------------------------------------
1302 
1304  const double relaxation)
1305  : relaxation(relaxation)
1306 {}
1307 
1308 
1310  : relaxation(0)
1311  , n_rows(0)
1312  , n_columns(0)
1313 {
1314  AdditionalData add_data;
1315  relaxation = add_data.relaxation;
1316 }
1317 
1318 
1319 
1320 inline void
1322  const PreconditionRichardson::AdditionalData &parameters)
1323 {
1324  relaxation = parameters.relaxation;
1325 }
1326 
1327 
1328 
1329 template <typename MatrixType>
1330 inline void
1332  const MatrixType & matrix,
1333  const PreconditionRichardson::AdditionalData &parameters)
1334 {
1335  relaxation = parameters.relaxation;
1336  n_rows = matrix.m();
1337  n_columns = matrix.n();
1338 }
1339 
1340 
1341 
1342 template <class VectorType>
1343 inline void
1344 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
1345 {
1346  static_assert(
1347  std::is_same<size_type, typename VectorType::size_type>::value,
1348  "PreconditionRichardson and VectorType must have the same size_type.");
1349 
1350  dst.equ(relaxation, src);
1351 }
1352 
1353 
1354 
1355 template <class VectorType>
1356 inline void
1357 PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
1358 {
1359  static_assert(
1360  std::is_same<size_type, typename VectorType::size_type>::value,
1361  "PreconditionRichardson and VectorType must have the same size_type.");
1362 
1363  dst.equ(relaxation, src);
1364 }
1365 
1366 template <class VectorType>
1367 inline void
1368 PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
1369 {
1370  static_assert(
1371  std::is_same<size_type, typename VectorType::size_type>::value,
1372  "PreconditionRichardson and VectorType must have the same size_type.");
1373 
1374  dst.add(relaxation, src);
1375 }
1376 
1377 
1378 
1379 template <class VectorType>
1380 inline void
1381 PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
1382 {
1383  static_assert(
1384  std::is_same<size_type, typename VectorType::size_type>::value,
1385  "PreconditionRichardson and VectorType must have the same size_type.");
1386 
1387  dst.add(relaxation, src);
1388 }
1389 
1392 {
1393  Assert(n_rows != 0, ExcNotInitialized());
1394  return n_rows;
1395 }
1396 
1399 {
1401  return n_columns;
1402 }
1403 
1404 //---------------------------------------------------------------------------
1405 
1406 template <typename MatrixType>
1407 inline void
1409  const AdditionalData &parameters)
1410 {
1411  A = &rA;
1412  relaxation = parameters.relaxation;
1413 }
1414 
1415 
1416 template <typename MatrixType>
1417 inline void
1419 {
1420  A = nullptr;
1421 }
1422 
1423 template <typename MatrixType>
1426 {
1427  Assert(A != nullptr, ExcNotInitialized());
1428  return A->m();
1429 }
1430 
1431 template <typename MatrixType>
1434 {
1435  Assert(A != nullptr, ExcNotInitialized());
1436  return A->n();
1437 }
1438 
1439 //---------------------------------------------------------------------------
1440 
1441 template <typename MatrixType>
1442 template <class VectorType>
1443 inline void
1444 PreconditionJacobi<MatrixType>::vmult(VectorType & dst,
1445  const VectorType &src) const
1446 {
1447  static_assert(
1448  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1449  typename VectorType::size_type>::value,
1450  "PreconditionJacobi and VectorType must have the same size_type.");
1451 
1452  Assert(this->A != nullptr, ExcNotInitialized());
1453  this->A->precondition_Jacobi(dst, src, this->relaxation);
1454 }
1455 
1456 
1457 
1458 template <typename MatrixType>
1459 template <class VectorType>
1460 inline void
1462  const VectorType &src) const
1463 {
1464  static_assert(
1465  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1466  typename VectorType::size_type>::value,
1467  "PreconditionJacobi and VectorType must have the same size_type.");
1468 
1469  Assert(this->A != nullptr, ExcNotInitialized());
1470  this->A->precondition_Jacobi(dst, src, this->relaxation);
1471 }
1472 
1473 
1474 
1475 template <typename MatrixType>
1476 template <class VectorType>
1477 inline void
1478 PreconditionJacobi<MatrixType>::step(VectorType & dst,
1479  const VectorType &src) const
1480 {
1481  static_assert(
1482  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1483  typename VectorType::size_type>::value,
1484  "PreconditionJacobi and VectorType must have the same size_type.");
1485 
1486  Assert(this->A != nullptr, ExcNotInitialized());
1487  this->A->Jacobi_step(dst, src, this->relaxation);
1488 }
1489 
1490 
1491 
1492 template <typename MatrixType>
1493 template <class VectorType>
1494 inline void
1495 PreconditionJacobi<MatrixType>::Tstep(VectorType & dst,
1496  const VectorType &src) const
1497 {
1498  static_assert(
1499  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1500  typename VectorType::size_type>::value,
1501  "PreconditionJacobi and VectorType must have the same size_type.");
1502 
1503  step(dst, src);
1504 }
1505 
1506 
1507 
1508 //---------------------------------------------------------------------------
1509 
1510 template <typename MatrixType>
1511 template <class VectorType>
1512 inline void
1513 PreconditionSOR<MatrixType>::vmult(VectorType &dst, const VectorType &src) const
1514 {
1515  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1516  typename VectorType::size_type>::value,
1517  "PreconditionSOR and VectorType must have the same size_type.");
1518 
1519  Assert(this->A != nullptr, ExcNotInitialized());
1520  this->A->precondition_SOR(dst, src, this->relaxation);
1521 }
1522 
1523 
1524 
1525 template <typename MatrixType>
1526 template <class VectorType>
1527 inline void
1528 PreconditionSOR<MatrixType>::Tvmult(VectorType & dst,
1529  const VectorType &src) const
1530 {
1531  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1532  typename VectorType::size_type>::value,
1533  "PreconditionSOR and VectorType must have the same size_type.");
1534 
1535  Assert(this->A != nullptr, ExcNotInitialized());
1536  this->A->precondition_TSOR(dst, src, this->relaxation);
1537 }
1538 
1539 
1540 
1541 template <typename MatrixType>
1542 template <class VectorType>
1543 inline void
1544 PreconditionSOR<MatrixType>::step(VectorType &dst, const VectorType &src) const
1545 {
1546  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1547  typename VectorType::size_type>::value,
1548  "PreconditionSOR and VectorType must have the same size_type.");
1549 
1550  Assert(this->A != nullptr, ExcNotInitialized());
1551  this->A->SOR_step(dst, src, this->relaxation);
1552 }
1553 
1554 
1555 
1556 template <typename MatrixType>
1557 template <class VectorType>
1558 inline void
1559 PreconditionSOR<MatrixType>::Tstep(VectorType &dst, const VectorType &src) const
1560 {
1561  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1562  typename VectorType::size_type>::value,
1563  "PreconditionSOR and VectorType must have the same size_type.");
1564 
1565  Assert(this->A != nullptr, ExcNotInitialized());
1566  this->A->TSOR_step(dst, src, this->relaxation);
1567 }
1568 
1569 
1570 
1571 //---------------------------------------------------------------------------
1572 
1573 template <typename MatrixType>
1574 inline void
1576  const MatrixType & rA,
1577  const typename BaseClass::AdditionalData &parameters)
1578 {
1579  this->PreconditionRelaxation<MatrixType>::initialize(rA, parameters);
1580 
1581  // in case we have a SparseMatrix class, we can extract information about
1582  // the diagonal.
1584  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
1585  &*this->A);
1586 
1587  // calculate the positions first after the diagonal.
1588  if (mat != nullptr)
1589  {
1590  const size_type n = this->A->n();
1591  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1592  for (size_type row = 0; row < n; ++row)
1593  {
1594  // find the first element in this line which is on the right of the
1595  // diagonal. we need to precondition with the elements on the left
1596  // only. note: the first entry in each line denotes the diagonal
1597  // element, which we need not check.
1599  it = mat->begin(row) + 1;
1600  for (; it < mat->end(row); ++it)
1601  if (it->column() > row)
1602  break;
1603  pos_right_of_diagonal[row] = it - mat->begin();
1604  }
1605  }
1606 }
1607 
1608 
1609 template <typename MatrixType>
1610 template <class VectorType>
1611 inline void
1612 PreconditionSSOR<MatrixType>::vmult(VectorType & dst,
1613  const VectorType &src) const
1614 {
1615  static_assert(
1616  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1617  typename VectorType::size_type>::value,
1618  "PreconditionSSOR and VectorType must have the same size_type.");
1619 
1620  Assert(this->A != nullptr, ExcNotInitialized());
1621  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1622 }
1623 
1624 
1625 
1626 template <typename MatrixType>
1627 template <class VectorType>
1628 inline void
1629 PreconditionSSOR<MatrixType>::Tvmult(VectorType & dst,
1630  const VectorType &src) const
1631 {
1632  static_assert(
1633  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1634  typename VectorType::size_type>::value,
1635  "PreconditionSSOR and VectorType must have the same size_type.");
1636 
1637  Assert(this->A != nullptr, ExcNotInitialized());
1638  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1639 }
1640 
1641 
1642 
1643 template <typename MatrixType>
1644 template <class VectorType>
1645 inline void
1646 PreconditionSSOR<MatrixType>::step(VectorType &dst, const VectorType &src) const
1647 {
1648  static_assert(
1649  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1650  typename VectorType::size_type>::value,
1651  "PreconditionSSOR and VectorType must have the same size_type.");
1652 
1653  Assert(this->A != nullptr, ExcNotInitialized());
1654  this->A->SSOR_step(dst, src, this->relaxation);
1655 }
1656 
1657 
1658 
1659 template <typename MatrixType>
1660 template <class VectorType>
1661 inline void
1662 PreconditionSSOR<MatrixType>::Tstep(VectorType & dst,
1663  const VectorType &src) const
1664 {
1665  static_assert(
1666  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1667  typename VectorType::size_type>::value,
1668  "PreconditionSSOR and VectorType must have the same size_type.");
1669 
1670  step(dst, src);
1671 }
1672 
1673 
1674 
1675 //---------------------------------------------------------------------------
1676 
1677 template <typename MatrixType>
1678 inline void
1680  const MatrixType & rA,
1681  const std::vector<size_type> & p,
1682  const std::vector<size_type> & ip,
1683  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1684 {
1685  permutation = &p;
1686  inverse_permutation = &ip;
1688 }
1689 
1690 
1691 template <typename MatrixType>
1692 inline void
1693 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
1694  const AdditionalData &additional_data)
1695 {
1696  initialize(A,
1697  additional_data.permutation,
1698  additional_data.inverse_permutation,
1699  additional_data.parameters);
1700 }
1701 
1702 
1703 template <typename MatrixType>
1704 template <typename VectorType>
1705 inline void
1706 PreconditionPSOR<MatrixType>::vmult(VectorType & dst,
1707  const VectorType &src) const
1708 {
1709  static_assert(
1710  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1711  typename VectorType::size_type>::value,
1712  "PreconditionPSOR and VectorType must have the same size_type.");
1713 
1714  Assert(this->A != nullptr, ExcNotInitialized());
1715  dst = src;
1716  this->A->PSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1717 }
1718 
1719 
1720 
1721 template <typename MatrixType>
1722 template <class VectorType>
1723 inline void
1724 PreconditionPSOR<MatrixType>::Tvmult(VectorType & dst,
1725  const VectorType &src) const
1726 {
1727  static_assert(
1728  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1729  typename VectorType::size_type>::value,
1730  "PreconditionPSOR and VectorType must have the same size_type.");
1731 
1732  Assert(this->A != nullptr, ExcNotInitialized());
1733  dst = src;
1734  this->A->TPSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1735 }
1736 
1737 template <typename MatrixType>
1739  const std::vector<size_type> &permutation,
1740  const std::vector<size_type> &inverse_permutation,
1741  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1742  : permutation(permutation)
1743  , inverse_permutation(inverse_permutation)
1744  , parameters(parameters)
1745 {}
1746 
1747 
1748 //---------------------------------------------------------------------------
1749 
1750 
1751 template <typename MatrixType, class VectorType>
1753  const MatrixType & M,
1754  const function_ptr method)
1755  : matrix(M)
1756  , precondition(method)
1757 {}
1758 
1759 
1760 
1761 template <typename MatrixType, class VectorType>
1762 void
1764  VectorType & dst,
1765  const VectorType &src) const
1766 {
1767  (matrix.*precondition)(dst, src);
1768 }
1769 
1770 //---------------------------------------------------------------------------
1771 
1772 template <typename MatrixType>
1774  const double relaxation)
1775  : relaxation(relaxation)
1776 {}
1777 
1778 
1779 
1780 //---------------------------------------------------------------------------
1781 
1782 namespace internal
1783 {
1784  namespace PreconditionChebyshevImplementation
1785  {
1786  // for deal.II vectors, perform updates for Chebyshev preconditioner all
1787  // at once to reduce memory transfer. Here, we select between general
1788  // vectors and deal.II vectors where we expand the loop over the (local)
1789  // size of the vector
1790 
1791  // generic part for non-deal.II vectors
1792  template <typename VectorType, typename PreconditionerType>
1793  inline void
1794  vector_updates(const VectorType & rhs,
1795  const PreconditionerType &preconditioner,
1796  const unsigned int iteration_index,
1797  const double factor1,
1798  const double factor2,
1799  VectorType & solution_old,
1800  VectorType & temp_vector1,
1801  VectorType & temp_vector2,
1802  VectorType & solution)
1803  {
1804  if (iteration_index == 0)
1805  {
1806  solution.equ(factor2, rhs);
1807  preconditioner.vmult(solution_old, solution);
1808  }
1809  else if (iteration_index == 1)
1810  {
1811  // compute t = P^{-1} * (b-A*x^{n})
1812  temp_vector1.sadd(-1.0, 1.0, rhs);
1813  preconditioner.vmult(solution_old, temp_vector1);
1814 
1815  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
1816  solution_old.sadd(factor2, 1 + factor1, solution);
1817  }
1818  else
1819  {
1820  // compute t = P^{-1} * (b-A*x^{n})
1821  temp_vector1.sadd(-1.0, 1.0, rhs);
1822  preconditioner.vmult(temp_vector2, temp_vector1);
1823 
1824  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
1825  solution_old.sadd(-factor1, factor2, temp_vector2);
1826  solution_old.add(1 + factor1, solution);
1827  }
1828 
1829  solution.swap(solution_old);
1830  }
1831 
1832  // worker routine for deal.II vectors. Because of vectorization, we need
1833  // to put the loop into an extra structure because the virtual function of
1834  // VectorUpdatesRange prevents the compiler from applying vectorization.
1835  template <typename Number>
1836  struct VectorUpdater
1837  {
1838  VectorUpdater(const Number * rhs,
1839  const Number * matrix_diagonal_inverse,
1840  const unsigned int iteration_index,
1841  const Number factor1,
1842  const Number factor2,
1843  Number * solution_old,
1844  Number * tmp_vector,
1845  Number * solution)
1846  : rhs(rhs)
1847  , matrix_diagonal_inverse(matrix_diagonal_inverse)
1848  , iteration_index(iteration_index)
1849  , factor1(factor1)
1850  , factor2(factor2)
1851  , solution_old(solution_old)
1852  , tmp_vector(tmp_vector)
1853  , solution(solution)
1854  {}
1855 
1856  void
1857  apply_to_subrange(const std::size_t begin, const std::size_t end) const
1858  {
1859  // To circumvent a bug in gcc
1860  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
1861  // copies of the variables factor1 and factor2 and do not check based on
1862  // factor1.
1863  const Number factor1 = this->factor1;
1864  const Number factor1_plus_1 = 1. + this->factor1;
1865  const Number factor2 = this->factor2;
1866  if (iteration_index == 0)
1867  {
1868  DEAL_II_OPENMP_SIMD_PRAGMA
1869  for (std::size_t i = begin; i < end; ++i)
1870  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
1871  }
1872  else if (iteration_index == 1)
1873  {
1874  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
1875  DEAL_II_OPENMP_SIMD_PRAGMA
1876  for (std::size_t i = begin; i < end; ++i)
1877  // for efficiency reason, write back to temp_vector that is
1878  // already read (avoid read-for-ownership)
1879  tmp_vector[i] =
1880  factor1_plus_1 * solution[i] +
1881  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1882  }
1883  else
1884  {
1885  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
1886  // + f_2 * P^{-1} * (b-A*x^{n})
1887  DEAL_II_OPENMP_SIMD_PRAGMA
1888  for (std::size_t i = begin; i < end; ++i)
1889  solution_old[i] =
1890  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
1891  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1892  }
1893  }
1894 
1895  const Number * rhs;
1896  const Number * matrix_diagonal_inverse;
1897  const unsigned int iteration_index;
1898  const Number factor1;
1899  const Number factor2;
1900  mutable Number * solution_old;
1901  mutable Number * tmp_vector;
1902  mutable Number * solution;
1903  };
1904 
1905  template <typename Number>
1906  struct VectorUpdatesRange : public parallel::ParallelForInteger
1907  {
1908  VectorUpdatesRange(const VectorUpdater<Number> &updater,
1909  const std::size_t size)
1910  : updater(updater)
1911  {
1912  if (size < internal::VectorImplementation::minimum_parallel_grain_size)
1913  apply_to_subrange(0, size);
1914  else
1915  apply_parallel(
1916  0,
1917  size,
1918  internal::VectorImplementation::minimum_parallel_grain_size);
1919  }
1920 
1921  ~VectorUpdatesRange() override = default;
1922 
1923  virtual void
1924  apply_to_subrange(const std::size_t begin,
1925  const std::size_t end) const override
1926  {
1927  updater.apply_to_subrange(begin, end);
1928  }
1929 
1930  const VectorUpdater<Number> &updater;
1931  };
1932 
1933  // selection for diagonal matrix around deal.II vector
1934  template <typename Number>
1935  inline void
1936  vector_updates(const ::Vector<Number> & rhs,
1937  const DiagonalMatrix<::Vector<Number>> &jacobi,
1938  const unsigned int iteration_index,
1939  const double factor1,
1940  const double factor2,
1941  ::Vector<Number> &solution_old,
1942  ::Vector<Number> &temp_vector1,
1943  ::Vector<Number> &,
1944  ::Vector<Number> &solution)
1945  {
1946  VectorUpdater<Number> upd(rhs.begin(),
1947  jacobi.get_vector().begin(),
1948  iteration_index,
1949  factor1,
1950  factor2,
1951  solution_old.begin(),
1952  temp_vector1.begin(),
1953  solution.begin());
1954  VectorUpdatesRange<Number>(upd, rhs.size());
1955 
1956  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
1957  if (iteration_index == 0)
1958  {
1959  // nothing to do here because we can immediately write into the
1960  // solution vector without remembering any of the other vectors
1961  }
1962  else if (iteration_index == 1)
1963  {
1964  solution.swap(temp_vector1);
1965  solution_old.swap(temp_vector1);
1966  }
1967  else
1968  solution.swap(solution_old);
1969  }
1970 
1971  // selection for diagonal matrix around parallel deal.II vector
1972  template <typename Number, typename MemorySpace>
1973  inline void
1974  vector_updates(
1976  const DiagonalMatrix<
1978  const unsigned int iteration_index,
1979  const double factor1,
1980  const double factor2,
1985  {
1986  VectorUpdater<Number> upd(rhs.begin(),
1987  jacobi.get_vector().begin(),
1988  iteration_index,
1989  factor1,
1990  factor2,
1991  solution_old.begin(),
1992  temp_vector1.begin(),
1993  solution.begin());
1994  VectorUpdatesRange<Number>(upd, rhs.local_size());
1995 
1996  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
1997  if (iteration_index == 0)
1998  {
1999  // nothing to do here because we can immediately write into the
2000  // solution vector without remembering any of the other vectors
2001  }
2002  else if (iteration_index == 1)
2003  {
2004  solution.swap(temp_vector1);
2005  solution_old.swap(temp_vector1);
2006  }
2007  else
2008  solution.swap(solution_old);
2009  }
2010 
2011  template <typename MatrixType,
2012  typename VectorType,
2013  typename PreconditionerType>
2014  inline void
2015  initialize_preconditioner(
2016  const MatrixType & matrix,
2017  std::shared_ptr<PreconditionerType> &preconditioner,
2018  VectorType &)
2019  {
2020  (void)matrix;
2021  (void)preconditioner;
2022  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2023  }
2024 
2025  template <typename MatrixType, typename VectorType>
2026  inline void
2027  initialize_preconditioner(
2028  const MatrixType & matrix,
2029  std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner,
2030  VectorType & diagonal_inverse)
2031  {
2032  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2033  {
2034  if (preconditioner.get() == nullptr)
2035  preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2036 
2037  Assert(
2038  preconditioner->m() == 0,
2039  ExcMessage(
2040  "Preconditioner appears to be initialized but not sized correctly"));
2041 
2042  // Check if we can initialize from vector that then gets set to zero
2043  // as the matrix will own the memory
2044  preconditioner->reinit(diagonal_inverse);
2045  {
2046  VectorType empty_vector;
2047  diagonal_inverse.reinit(empty_vector);
2048  }
2049 
2050  // This part only works in serial
2051  if (preconditioner->m() != matrix.m())
2052  {
2053  preconditioner->get_vector().reinit(matrix.m());
2054  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2055  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2056  }
2057  }
2058  }
2059 
2060  template <typename VectorType>
2061  void
2062  set_initial_guess(VectorType &vector)
2063  {
2064  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2065  if (vector.locally_owned_elements().is_element(0))
2066  vector(0) = 0.;
2067  }
2068 
2069  template <typename Number>
2070  void
2071  set_initial_guess(::Vector<Number> &vector)
2072  {
2073  // Choose a high-frequency mode consisting of numbers between 0 and 1
2074  // that is cheap to compute (cheaper than random numbers) but avoids
2075  // obviously re-occurring numbers in multi-component systems by choosing
2076  // a period of 11
2077  for (unsigned int i = 0; i < vector.size(); ++i)
2078  vector(i) = i % 11;
2079 
2080  const Number mean_value = vector.mean_value();
2081  vector.add(-mean_value);
2082  }
2083 
2084  template <typename Number, typename MemorySpace>
2085  void
2086  set_initial_guess(
2088  {
2089  // Choose a high-frequency mode consisting of numbers between 0 and 1
2090  // that is cheap to compute (cheaper than random numbers) but avoids
2091  // obviously re-occurring numbers in multi-component systems by choosing
2092  // a period of 11.
2093  // Make initial guess robust with respect to number of processors
2094  // by operating on the global index.
2095  types::global_dof_index first_local_range = 0;
2096  if (!vector.locally_owned_elements().is_empty())
2097  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2098  for (unsigned int i = 0; i < vector.local_size(); ++i)
2099  vector.local_element(i) = (i + first_local_range) % 11;
2100 
2101  const Number mean_value = vector.mean_value();
2102  vector.add(-mean_value);
2103  }
2104 
2105  struct EigenvalueTracker
2106  {
2107  public:
2108  void
2109  slot(const std::vector<double> &eigenvalues)
2110  {
2111  values = eigenvalues;
2112  }
2113 
2114  std::vector<double> values;
2115  };
2116  } // namespace PreconditionChebyshevImplementation
2117 } // namespace internal
2118 
2119 
2120 
2121 // avoid warning about deprecated variable nonzero_starting
2122 
2123 template <typename MatrixType, class VectorType, typename PreconditionerType>
2125  AdditionalData::AdditionalData(const unsigned int degree,
2126  const double smoothing_range,
2127  const bool nonzero_starting,
2128  const unsigned int eig_cg_n_iterations,
2129  const double eig_cg_residual,
2130  const double max_eigenvalue)
2131  : degree(degree)
2132  , smoothing_range(smoothing_range)
2133  , nonzero_starting(nonzero_starting)
2134  , eig_cg_n_iterations(eig_cg_n_iterations)
2135  , eig_cg_residual(eig_cg_residual)
2136  , max_eigenvalue(max_eigenvalue)
2137 {}
2138 
2139 
2140 
2141 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2144  : theta(1.)
2145  , delta(1.)
2146  , eigenvalues_are_initialized(false)
2147 {
2148  static_assert(
2149  std::is_same<size_type, typename VectorType::size_type>::value,
2150  "PreconditionChebyshev and VectorType must have the same size_type.");
2151 }
2152 
2153 
2154 // avoid warning about deprecated variable
2155 // AdditionalData::matrix_diagonal_inverse
2156 
2157 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2158 inline void
2160  const MatrixType & matrix,
2161  const AdditionalData &additional_data)
2162 {
2163  matrix_ptr = &matrix;
2164  data = additional_data;
2165  Assert(data.degree > 0,
2166  ExcMessage("The degree of the Chebyshev method must be positive."));
2167  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2168  matrix, data.preconditioner, data.matrix_diagonal_inverse);
2169  eigenvalues_are_initialized = false;
2170 }
2171 
2172 
2173 
2174 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2175 inline void
2177 {
2178  eigenvalues_are_initialized = false;
2179  theta = delta = 1.0;
2180  matrix_ptr = nullptr;
2181  {
2182  VectorType empty_vector;
2183  data.matrix_diagonal_inverse.reinit(empty_vector);
2184  solution_old.reinit(empty_vector);
2185  temp_vector1.reinit(empty_vector);
2186  temp_vector2.reinit(empty_vector);
2187  }
2188  data.preconditioner.reset();
2189 }
2190 
2191 
2192 
2193 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2194 inline void
2196  estimate_eigenvalues(const VectorType &src) const
2197 {
2198  Assert(eigenvalues_are_initialized == false, ExcInternalError());
2199  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2200 
2201  solution_old.reinit(src);
2202  temp_vector1.reinit(src, true);
2203 
2204  // calculate largest eigenvalue using a hand-tuned CG iteration on the
2205  // matrix weighted by its diagonal. we start with a vector that consists of
2206  // ones only, weighted by the length.
2207  double max_eigenvalue, min_eigenvalue;
2208  if (data.eig_cg_n_iterations > 0)
2209  {
2210  Assert(data.eig_cg_n_iterations > 2,
2211  ExcMessage(
2212  "Need to set at least two iterations to find eigenvalues."));
2213 
2214  // set a very strict tolerance to force at least two iterations
2215  ReductionControl control(
2216  data.eig_cg_n_iterations,
2217  std::sqrt(
2218  std::numeric_limits<typename VectorType::value_type>::epsilon()),
2219  1e-10,
2220  false,
2221  false);
2222 
2223  internal::PreconditionChebyshevImplementation::EigenvalueTracker
2224  eigenvalue_tracker;
2225  SolverCG<VectorType> solver(control);
2226  solver.connect_eigenvalues_slot(std::bind(
2227  &internal::PreconditionChebyshevImplementation::EigenvalueTracker::slot,
2228  &eigenvalue_tracker,
2229  std::placeholders::_1));
2230 
2231  // set an initial guess which is close to the constant vector but where
2232  // one entry is different to trigger high frequencies
2233  internal::PreconditionChebyshevImplementation::set_initial_guess(
2234  temp_vector1);
2235 
2236  try
2237  {
2238  solver.solve(*matrix_ptr,
2239  solution_old,
2240  temp_vector1,
2241  *data.preconditioner);
2242  }
2244  {}
2245 
2246  // read the eigenvalues from the attached eigenvalue tracker
2247  if (eigenvalue_tracker.values.empty())
2248  min_eigenvalue = max_eigenvalue = 1;
2249  else
2250  {
2251  min_eigenvalue = eigenvalue_tracker.values.front();
2252 
2253  // include a safety factor since the CG method will in general not
2254  // be converged
2255  max_eigenvalue = 1.2 * eigenvalue_tracker.values.back();
2256  }
2257  }
2258  else
2259  {
2260  max_eigenvalue = data.max_eigenvalue;
2261  min_eigenvalue = data.max_eigenvalue / data.smoothing_range;
2262  }
2263 
2264  const double alpha = (data.smoothing_range > 1. ?
2265  max_eigenvalue / data.smoothing_range :
2266  std::min(0.9 * max_eigenvalue, min_eigenvalue));
2267 
2268  // in case the user set the degree to invalid unsigned int, we have to
2269  // determine the number of necessary iterations from the Chebyshev error
2270  // estimate, given the target tolerance specified by smoothing_range. This
2271  // estimate is based on the error formula given in section 5.1 of
2272  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
2273  if (data.degree == numbers::invalid_unsigned_int)
2274  {
2275  const double actual_range = max_eigenvalue / alpha;
2276  const double sigma = (1. - std::sqrt(1. / actual_range)) /
2277  (1. + std::sqrt(1. / actual_range));
2278  const double eps = data.smoothing_range;
2279  const_cast<
2281  this)
2282  ->data.degree =
2283  1 + static_cast<unsigned int>(
2284  std::log(1. / eps + std::sqrt(1. / eps / eps - 1)) /
2285  std::log(1. / sigma));
2286  }
2287 
2288  const_cast<
2290  ->delta = (max_eigenvalue - alpha) * 0.5;
2291  const_cast<
2293  ->theta = (max_eigenvalue + alpha) * 0.5;
2294 
2295  // We do not need the second temporary vector in case we have a
2296  // DiagonalMatrix as preconditioner and use deal.II's own vectors
2297  using NumberType = typename VectorType::value_type;
2298  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
2299  false ||
2300  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
2301  ((std::is_same<VectorType,
2303  Vector<NumberType, MemorySpace::Host>>::value ==
2304  false) ||
2305  (std::is_same<VectorType,
2306  LinearAlgebra::distributed::
2307  Vector<NumberType, MemorySpace::CUDA>>::value ==
2308  false))))
2309  temp_vector2.reinit(src, true);
2310  else
2311  {
2312  VectorType empty_vector;
2313  temp_vector2.reinit(empty_vector);
2314  }
2315 
2316  const_cast<
2318  ->eigenvalues_are_initialized = true;
2319 }
2320 
2321 
2322 
2323 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2324 inline void
2326  VectorType & solution,
2327  const VectorType &rhs) const
2328 {
2329  std::lock_guard<std::mutex> lock(mutex);
2330  if (eigenvalues_are_initialized == false)
2331  estimate_eigenvalues(rhs);
2332 
2333  internal::PreconditionChebyshevImplementation::vector_updates(
2334  rhs,
2335  *data.preconditioner,
2336  0,
2337  0.,
2338  1. / theta,
2339  solution_old,
2340  temp_vector1,
2341  temp_vector2,
2342  solution);
2343 
2344  // if delta is zero, we do not need to iterate because the updates will be
2345  // zero
2346  if (data.degree < 2 || std::abs(delta) < 1e-40)
2347  return;
2348 
2349  double rhok = delta / theta, sigma = theta / delta;
2350  for (unsigned int k = 0; k < data.degree - 1; ++k)
2351  {
2352  matrix_ptr->vmult(temp_vector1, solution);
2353  const double rhokp = 1. / (2. * sigma - rhok);
2354  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2355  rhok = rhokp;
2356  internal::PreconditionChebyshevImplementation::vector_updates(
2357  rhs,
2358  *data.preconditioner,
2359  k + 1,
2360  factor1,
2361  factor2,
2362  solution_old,
2363  temp_vector1,
2364  temp_vector2,
2365  solution);
2366  }
2367 }
2368 
2369 
2370 
2371 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2372 inline void
2374  VectorType & solution,
2375  const VectorType &rhs) const
2376 {
2377  std::lock_guard<std::mutex> lock(mutex);
2378  if (eigenvalues_are_initialized == false)
2379  estimate_eigenvalues(rhs);
2380 
2381  internal::PreconditionChebyshevImplementation::vector_updates(
2382  rhs,
2383  *data.preconditioner,
2384  0,
2385  0.,
2386  1. / theta,
2387  solution_old,
2388  temp_vector1,
2389  temp_vector2,
2390  solution);
2391 
2392  if (data.degree < 2 || std::abs(delta) < 1e-40)
2393  return;
2394 
2395  double rhok = delta / theta, sigma = theta / delta;
2396  for (unsigned int k = 0; k < data.degree - 1; ++k)
2397  {
2398  matrix_ptr->Tvmult(temp_vector1, solution);
2399  const double rhokp = 1. / (2. * sigma - rhok);
2400  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2401  rhok = rhokp;
2402  internal::PreconditionChebyshevImplementation::vector_updates(
2403  rhs,
2404  *data.preconditioner,
2405  k + 1,
2406  factor1,
2407  factor2,
2408  solution_old,
2409  temp_vector1,
2410  temp_vector2,
2411  solution);
2412  }
2413 }
2414 
2415 
2416 
2417 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2418 inline void
2420  VectorType & solution,
2421  const VectorType &rhs) const
2422 {
2423  std::lock_guard<std::mutex> lock(mutex);
2424  if (eigenvalues_are_initialized == false)
2425  estimate_eigenvalues(rhs);
2426 
2427  matrix_ptr->vmult(temp_vector1, solution);
2428  internal::PreconditionChebyshevImplementation::vector_updates(
2429  rhs,
2430  *data.preconditioner,
2431  1,
2432  0.,
2433  1. / theta,
2434  solution_old,
2435  temp_vector1,
2436  temp_vector2,
2437  solution);
2438 
2439  if (data.degree < 2 || std::abs(delta) < 1e-40)
2440  return;
2441 
2442  double rhok = delta / theta, sigma = theta / delta;
2443  for (unsigned int k = 0; k < data.degree - 1; ++k)
2444  {
2445  matrix_ptr->vmult(temp_vector1, solution);
2446  const double rhokp = 1. / (2. * sigma - rhok);
2447  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2448  rhok = rhokp;
2449  internal::PreconditionChebyshevImplementation::vector_updates(
2450  rhs,
2451  *data.preconditioner,
2452  k + 2,
2453  factor1,
2454  factor2,
2455  solution_old,
2456  temp_vector1,
2457  temp_vector2,
2458  solution);
2459  }
2460 }
2461 
2462 
2463 
2464 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2465 inline void
2467  VectorType & solution,
2468  const VectorType &rhs) const
2469 {
2470  std::lock_guard<std::mutex> lock(mutex);
2471  if (eigenvalues_are_initialized == false)
2472  estimate_eigenvalues(rhs);
2473 
2474  matrix_ptr->Tvmult(temp_vector1, solution);
2475  internal::PreconditionChebyshevImplementation::vector_updates(
2476  rhs,
2477  *data.preconditioner,
2478  1,
2479  0.,
2480  1. / theta,
2481  solution_old,
2482  temp_vector1,
2483  temp_vector2,
2484  solution);
2485 
2486  if (data.degree < 2 || std::abs(delta) < 1e-40)
2487  return;
2488 
2489  double rhok = delta / theta, sigma = theta / delta;
2490  for (unsigned int k = 0; k < data.degree - 1; ++k)
2491  {
2492  matrix_ptr->Tvmult(temp_vector1, solution);
2493  const double rhokp = 1. / (2. * sigma - rhok);
2494  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2495  rhok = rhokp;
2496  internal::PreconditionChebyshevImplementation::vector_updates(
2497  rhs,
2498  *data.preconditioner,
2499  k + 2,
2500  factor1,
2501  factor2,
2502  solution_old,
2503  temp_vector1,
2504  temp_vector2,
2505  solution);
2506  }
2507 }
2508 
2509 
2510 
2511 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2512 inline typename PreconditionChebyshev<MatrixType,
2513  VectorType,
2514  PreconditionerType>::size_type
2516 {
2517  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2518  return matrix_ptr->m();
2519 }
2520 
2521 
2522 
2523 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2524 inline typename PreconditionChebyshev<MatrixType,
2525  VectorType,
2526  PreconditionerType>::size_type
2528 {
2529  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2530  return matrix_ptr->n();
2531 }
2532 
2533 #endif // DOXYGEN
2534 
2535 DEAL_II_NAMESPACE_CLOSE
2536 
2537 #endif
void Tvmult(VectorType &, const VectorType &) const
static const unsigned int invalid_unsigned_int
Definition: types.h:173
size_type m() const
types::global_dof_index size_type
void initialize(const MatrixType &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
const_iterator end() const
Contents is actually a matrix.
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1780
void Tstep(VectorType &x, const VectorType &rhs) const
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
void vmult(VectorType &dst, const VectorType &src) const
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionRelaxation< MatrixType >::AdditionalData parameters
Definition: precondition.h:807
void vmult(VectorType &, const VectorType &) const
size_type m() const
const MatrixType & matrix
Definition: precondition.h:390
AdditionalData data
void reinit(const size_type size, const bool omit_zeroing_entries=false)
AdditionalData(const double relaxation=1.)
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition: precondition.h:370
static ::ExceptionBase & ExcNotInitialized()
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
std::shared_ptr< PreconditionerType > preconditioner
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
AdditionalData(const double relaxation=1.)
size_type m() const
const std::vector< size_type > * permutation
Definition: precondition.h:860
LinearAlgebra::distributed::Vector< Number > Vector
const function_ptr precondition
Definition: precondition.h:395
void vmult(VectorType &, const VectorType &) const
void initialize(const AdditionalData &parameters)
Number local_element(const size_type local_index) const
void Tvmult(VectorType &, const VectorType &) const
static ::ExceptionBase & ExcMessage(std::string arg1)
typename MatrixType::size_type size_type
Definition: precondition.h:415
void Tvmult(VectorType &dst, const VectorType &src) const
size_type n() const
void vmult(VectorType &, const VectorType &) const
const std::vector< size_type > & inverse_permutation
Definition: precondition.h:803
#define Assert(cond, exc)
Definition: exceptions.h:1407
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:670
void step(VectorType &x, const VectorType &rhs) const
typename MatrixType::size_type size_type
Definition: precondition.h:771
AdditionalData(const unsigned int degree=1, const double smoothing_range=0., const bool nonzero_starting=false, const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1)
void vmult_add(VectorType &, const VectorType &) const
void Tvmult(VectorType &, const VectorType &) const
size_type n() const
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
Threads::Mutex mutex
void Tvmult(VectorType &, const VectorType &) const
void estimate_eigenvalues(const VectorType &src) const
void Tstep(VectorType &x, const VectorType &rhs) const
void vmult(VectorType &, const VectorType &) const
void Tstep(VectorType &x, const VectorType &rhs) const
void vmult(VectorType &dst, const VectorType &src) const
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename PreconditionRelaxation< MatrixType >::AdditionalData &parameters=typename PreconditionRelaxation< MatrixType >::AdditionalData())
void Tvmult_add(VectorType &, const VectorType &) const
void Tvmult(VectorType &, const VectorType &) const
types::global_dof_index size_type
Definition: precondition.h:203
void vmult(VectorType &, const VectorType &) const
unsigned int global_dof_index
Definition: types.h:89
virtual ::IndexSet locally_owned_elements() const override
size_type n() const
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
Definition: precondition.h:467
void Tvmult(VectorType &, const VectorType &) const
void step(VectorType &x, const VectorType &rhs) const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
size_type n() const
void swap(Vector< Number, MemorySpace > &v)
typename MatrixType::size_type size_type
Definition: precondition.h:675
const std::vector< size_type > * inverse_permutation
Definition: precondition.h:864
void Tstep(VectorType &dst, const VectorType &src) const
void step(VectorType &dst, const VectorType &src) const
void step(VectorType &x, const VectorType &rhs) const
virtual void add(const Number a) override
types::global_dof_index size_type
Definition: precondition.h:84
size_type m() const
const_iterator begin() const
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:513
void vmult(VectorType &, const VectorType &) const
const std::vector< size_type > & permutation
Definition: precondition.h:799
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:601
virtual Number mean_value() const override
void vmult_add(VectorType &, const VectorType &) const
std::vector< std::size_t > pos_right_of_diagonal
Definition: precondition.h:728
bool is_empty() const
Definition: index_set.h:1724
SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
static ::ExceptionBase & ExcInternalError()
void initialize(const MatrixType &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename PreconditionRelaxation< MatrixType >::AdditionalData &parameters=typename PreconditionRelaxation< MatrixType >::AdditionalData())