Reference documentation for deal.II version 9.3.0
\(\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\}}\)
precondition.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2021 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/cuda_size.h>
25 #include <deal.II/base/parallel.h>
29 #include <deal.II/base/utilities.h>
30 
33 #include <deal.II/lac/solver_cg.h>
35 
37 
38 // forward declarations
39 #ifndef DOXYGEN
40 template <typename number>
41 class Vector;
42 template <typename number>
43 class SparseMatrix;
44 namespace LinearAlgebra
45 {
46  namespace distributed
47  {
48  template <typename, typename>
49  class Vector;
50  template <typename>
51  class BlockVector;
52  } // namespace distributed
53 } // namespace LinearAlgebra
54 #endif
55 
56 
80 {
81 public:
86 
92  {
96  AdditionalData() = default;
97  };
98 
103 
108  template <typename MatrixType>
109  void
110  initialize(const MatrixType & matrix,
111  const AdditionalData &additional_data = AdditionalData());
112 
116  template <class VectorType>
117  void
118  vmult(VectorType &, const VectorType &) const;
119 
124  template <class VectorType>
125  void
126  Tvmult(VectorType &, const VectorType &) const;
127 
131  template <class VectorType>
132  void
133  vmult_add(VectorType &, const VectorType &) const;
134 
139  template <class VectorType>
140  void
141  Tvmult_add(VectorType &, const VectorType &) const;
142 
147  void
149  {}
150 
158  size_type
159  m() const;
160 
168  size_type
169  n() const;
170 
171 private:
176 
181 };
182 
183 
184 
196 {
197 public:
202 
207  {
208  public:
213  AdditionalData(const double relaxation = 1.);
214 
218  double relaxation;
219  };
220 
226 
230  void
231  initialize(const AdditionalData &parameters);
232 
238  template <typename MatrixType>
239  void
240  initialize(const MatrixType &matrix, const AdditionalData &parameters);
241 
245  template <class VectorType>
246  void
247  vmult(VectorType &, const VectorType &) const;
248 
253  template <class VectorType>
254  void
255  Tvmult(VectorType &, const VectorType &) const;
259  template <class VectorType>
260  void
261  vmult_add(VectorType &, const VectorType &) const;
262 
267  template <class VectorType>
268  void
269  Tvmult_add(VectorType &, const VectorType &) const;
270 
275  void
277  {}
278 
286  size_type
287  m() const;
288 
296  size_type
297  n() const;
298 
299 private:
303  double relaxation;
304 
309 
314 };
315 
316 
317 
357 template <typename MatrixType = SparseMatrix<double>,
358  class VectorType = Vector<double>>
360 {
361 public:
365  using function_ptr = void (MatrixType::*)(VectorType &,
366  const VectorType &) const;
367 
373  PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
374 
379  void
380  vmult(VectorType &dst, const VectorType &src) const;
381 
382 private:
386  const MatrixType &matrix;
387 
392 };
393 
394 
395 
401 template <typename MatrixType = SparseMatrix<double>>
403 {
404 public:
409 
414  {
415  public:
419  AdditionalData(const double relaxation = 1.);
420 
424  double relaxation;
425  };
426 
432  void
433  initialize(const MatrixType & A,
434  const AdditionalData &parameters = AdditionalData());
435 
439  void
440  clear();
441 
446  size_type
447  m() const;
448 
453  size_type
454  n() const;
455 
456 protected:
461 
465  double relaxation;
466 };
467 
468 
469 
496 template <typename MatrixType = SparseMatrix<double>>
497 class PreconditionJacobi : public PreconditionRelaxation<MatrixType>
498 {
499 public:
503  using AdditionalData =
505 
509  template <class VectorType>
510  void
511  vmult(VectorType &, const VectorType &) const;
512 
517  template <class VectorType>
518  void
519  Tvmult(VectorType &, const VectorType &) const;
520 
524  template <class VectorType>
525  void
526  step(VectorType &x, const VectorType &rhs) const;
527 
531  template <class VectorType>
532  void
533  Tstep(VectorType &x, const VectorType &rhs) const;
534 };
535 
536 
582 template <typename MatrixType = SparseMatrix<double>>
583 class PreconditionSOR : public PreconditionRelaxation<MatrixType>
584 {
585 public:
589  using AdditionalData =
591 
595  template <class VectorType>
596  void
597  vmult(VectorType &, const VectorType &) const;
598 
602  template <class VectorType>
603  void
604  Tvmult(VectorType &, const VectorType &) const;
605 
609  template <class VectorType>
610  void
611  step(VectorType &x, const VectorType &rhs) const;
612 
616  template <class VectorType>
617  void
618  Tstep(VectorType &x, const VectorType &rhs) const;
619 };
620 
621 
622 
649 template <typename MatrixType = SparseMatrix<double>>
650 class PreconditionSSOR : public PreconditionRelaxation<MatrixType>
651 {
652 public:
656  using AdditionalData =
658 
663 
668 
669 
675  void
676  initialize(const MatrixType & A,
677  const typename BaseClass::AdditionalData &parameters =
678  typename BaseClass::AdditionalData());
679 
683  template <class VectorType>
684  void
685  vmult(VectorType &, const VectorType &) const;
686 
691  template <class VectorType>
692  void
693  Tvmult(VectorType &, const VectorType &) const;
694 
695 
699  template <class VectorType>
700  void
701  step(VectorType &x, const VectorType &rhs) const;
702 
706  template <class VectorType>
707  void
708  Tstep(VectorType &x, const VectorType &rhs) const;
709 
710 private:
715  std::vector<std::size_t> pos_right_of_diagonal;
716 };
717 
718 
748 template <typename MatrixType = SparseMatrix<double>>
749 class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
750 {
751 public:
756 
761  {
762  public:
774  const std::vector<size_type> &permutation,
775  const std::vector<size_type> &inverse_permutation,
777  &parameters =
779 
783  const std::vector<size_type> &permutation;
787  const std::vector<size_type> &inverse_permutation;
792  };
793 
805  void
806  initialize(const MatrixType & A,
807  const std::vector<size_type> &permutation,
808  const std::vector<size_type> &inverse_permutation,
810  &parameters =
812 
823  void
824  initialize(const MatrixType &A, const AdditionalData &additional_data);
825 
829  template <class VectorType>
830  void
831  vmult(VectorType &, const VectorType &) const;
832 
836  template <class VectorType>
837  void
838  Tvmult(VectorType &, const VectorType &) const;
839 
840 private:
844  const std::vector<size_type> *permutation;
848  const std::vector<size_type> *inverse_permutation;
849 };
850 
851 
852 
988 template <typename MatrixType = SparseMatrix<double>,
989  typename VectorType = Vector<double>,
990  typename PreconditionerType = DiagonalMatrix<VectorType>>
992 {
993 public:
998 
1004  {
1008  AdditionalData(const unsigned int degree = 1,
1009  const double smoothing_range = 0.,
1010  const unsigned int eig_cg_n_iterations = 8,
1011  const double eig_cg_residual = 1e-2,
1012  const double max_eigenvalue = 1);
1013 
1017  AdditionalData &
1018  operator=(const AdditionalData &other_data);
1019 
1031  unsigned int degree;
1032 
1045 
1052  unsigned int eig_cg_n_iterations;
1053 
1059 
1066 
1072 
1076  std::shared_ptr<PreconditionerType> preconditioner;
1077  };
1078 
1079 
1084 
1096  void
1097  initialize(const MatrixType & matrix,
1098  const AdditionalData &additional_data = AdditionalData());
1099 
1104  void
1105  vmult(VectorType &dst, const VectorType &src) const;
1106 
1111  void
1112  Tvmult(VectorType &dst, const VectorType &src) const;
1113 
1117  void
1118  step(VectorType &dst, const VectorType &src) const;
1119 
1123  void
1124  Tstep(VectorType &dst, const VectorType &src) const;
1125 
1129  void
1130  clear();
1131 
1136  size_type
1137  m() const;
1138 
1143  size_type
1144  n() const;
1145 
1151  {
1163  unsigned int cg_iterations;
1168  unsigned int degree;
1173  : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
1174  , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
1175  , cg_iterations{0}
1176  , degree{0}
1177  {}
1178  };
1179 
1193  estimate_eigenvalues(const VectorType &src) const;
1194 
1195 private:
1199  SmartPointer<
1200  const MatrixType,
1203 
1208 
1213 
1218 
1224 
1228  double theta;
1229 
1234  double delta;
1235 
1241 
1247 };
1248 
1249 
1250 
1252 /* ---------------------------------- Inline functions ------------------- */
1253 
1254 #ifndef DOXYGEN
1255 
1257  : n_rows(0)
1258  , n_columns(0)
1259 {}
1260 
1261 template <typename MatrixType>
1262 inline void
1263 PreconditionIdentity::initialize(const MatrixType &matrix,
1265 {
1266  n_rows = matrix.m();
1267  n_columns = matrix.n();
1268 }
1269 
1270 
1271 template <class VectorType>
1272 inline void
1273 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1274 {
1275  dst = src;
1276 }
1277 
1278 
1279 
1280 template <class VectorType>
1281 inline void
1282 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1283 {
1284  dst = src;
1285 }
1286 
1287 template <class VectorType>
1288 inline void
1290 {
1291  dst += src;
1292 }
1293 
1294 
1295 
1296 template <class VectorType>
1297 inline void
1299 {
1300  dst += src;
1301 }
1302 
1305 {
1306  Assert(n_rows != 0, ExcNotInitialized());
1307  return n_rows;
1308 }
1309 
1312 {
1313  Assert(n_columns != 0, ExcNotInitialized());
1314  return n_columns;
1315 }
1316 
1317 //---------------------------------------------------------------------------
1318 
1320  const double relaxation)
1321  : relaxation(relaxation)
1322 {}
1323 
1324 
1326  : relaxation(0)
1327  , n_rows(0)
1328  , n_columns(0)
1329 {
1330  AdditionalData add_data;
1331  relaxation = add_data.relaxation;
1332 }
1333 
1334 
1335 
1336 inline void
1338  const PreconditionRichardson::AdditionalData &parameters)
1339 {
1340  relaxation = parameters.relaxation;
1341 }
1342 
1343 
1344 
1345 template <typename MatrixType>
1346 inline void
1348  const MatrixType & matrix,
1349  const PreconditionRichardson::AdditionalData &parameters)
1350 {
1351  relaxation = parameters.relaxation;
1352  n_rows = matrix.m();
1353  n_columns = matrix.n();
1354 }
1355 
1356 
1357 
1358 template <class VectorType>
1359 inline void
1360 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
1361 {
1362  static_assert(
1363  std::is_same<size_type, typename VectorType::size_type>::value,
1364  "PreconditionRichardson and VectorType must have the same size_type.");
1365 
1366  dst.equ(relaxation, src);
1367 }
1368 
1369 
1370 
1371 template <class VectorType>
1372 inline void
1374 {
1375  static_assert(
1376  std::is_same<size_type, typename VectorType::size_type>::value,
1377  "PreconditionRichardson and VectorType must have the same size_type.");
1378 
1379  dst.equ(relaxation, src);
1380 }
1381 
1382 template <class VectorType>
1383 inline void
1385 {
1386  static_assert(
1387  std::is_same<size_type, typename VectorType::size_type>::value,
1388  "PreconditionRichardson and VectorType must have the same size_type.");
1389 
1390  dst.add(relaxation, src);
1391 }
1392 
1393 
1394 
1395 template <class VectorType>
1396 inline void
1398 {
1399  static_assert(
1400  std::is_same<size_type, typename VectorType::size_type>::value,
1401  "PreconditionRichardson and VectorType must have the same size_type.");
1402 
1403  dst.add(relaxation, src);
1404 }
1405 
1408 {
1409  Assert(n_rows != 0, ExcNotInitialized());
1410  return n_rows;
1411 }
1412 
1415 {
1416  Assert(n_columns != 0, ExcNotInitialized());
1417  return n_columns;
1418 }
1419 
1420 //---------------------------------------------------------------------------
1421 
1422 template <typename MatrixType>
1423 inline void
1425  const AdditionalData &parameters)
1426 {
1427  A = &rA;
1428  relaxation = parameters.relaxation;
1429 }
1430 
1431 
1432 template <typename MatrixType>
1433 inline void
1435 {
1436  A = nullptr;
1437 }
1438 
1439 template <typename MatrixType>
1442 {
1443  Assert(A != nullptr, ExcNotInitialized());
1444  return A->m();
1445 }
1446 
1447 template <typename MatrixType>
1450 {
1451  Assert(A != nullptr, ExcNotInitialized());
1452  return A->n();
1453 }
1454 
1455 //---------------------------------------------------------------------------
1456 
1457 template <typename MatrixType>
1458 template <class VectorType>
1459 inline void
1461  const VectorType &src) const
1462 {
1463  static_assert(
1464  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1465  typename VectorType::size_type>::value,
1466  "PreconditionJacobi and VectorType must have the same size_type.");
1467 
1468  Assert(this->A != nullptr, ExcNotInitialized());
1469  this->A->precondition_Jacobi(dst, src, this->relaxation);
1470 }
1471 
1472 
1473 
1474 template <typename MatrixType>
1475 template <class VectorType>
1476 inline void
1478  const VectorType &src) const
1479 {
1480  static_assert(
1481  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1482  typename VectorType::size_type>::value,
1483  "PreconditionJacobi and VectorType must have the same size_type.");
1484 
1485  Assert(this->A != nullptr, ExcNotInitialized());
1486  this->A->precondition_Jacobi(dst, src, this->relaxation);
1487 }
1488 
1489 
1490 
1491 template <typename MatrixType>
1492 template <class VectorType>
1493 inline void
1495  const VectorType &src) const
1496 {
1497  static_assert(
1498  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1499  typename VectorType::size_type>::value,
1500  "PreconditionJacobi and VectorType must have the same size_type.");
1501 
1502  Assert(this->A != nullptr, ExcNotInitialized());
1503  this->A->Jacobi_step(dst, src, this->relaxation);
1504 }
1505 
1506 
1507 
1508 template <typename MatrixType>
1509 template <class VectorType>
1510 inline void
1512  const VectorType &src) const
1513 {
1514  static_assert(
1515  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1516  typename VectorType::size_type>::value,
1517  "PreconditionJacobi and VectorType must have the same size_type.");
1518 
1519  step(dst, src);
1520 }
1521 
1522 
1523 
1524 //---------------------------------------------------------------------------
1525 
1526 template <typename MatrixType>
1527 template <class VectorType>
1528 inline void
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_SOR(dst, src, this->relaxation);
1537 }
1538 
1539 
1540 
1541 template <typename MatrixType>
1542 template <class VectorType>
1543 inline void
1545  const VectorType &src) const
1546 {
1547  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1548  typename VectorType::size_type>::value,
1549  "PreconditionSOR and VectorType must have the same size_type.");
1550 
1551  Assert(this->A != nullptr, ExcNotInitialized());
1552  this->A->precondition_TSOR(dst, src, this->relaxation);
1553 }
1554 
1555 
1556 
1557 template <typename MatrixType>
1558 template <class VectorType>
1559 inline void
1561 {
1562  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1563  typename VectorType::size_type>::value,
1564  "PreconditionSOR and VectorType must have the same size_type.");
1565 
1566  Assert(this->A != nullptr, ExcNotInitialized());
1567  this->A->SOR_step(dst, src, this->relaxation);
1568 }
1569 
1570 
1571 
1572 template <typename MatrixType>
1573 template <class VectorType>
1574 inline void
1576 {
1577  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1578  typename VectorType::size_type>::value,
1579  "PreconditionSOR and VectorType must have the same size_type.");
1580 
1581  Assert(this->A != nullptr, ExcNotInitialized());
1582  this->A->TSOR_step(dst, src, this->relaxation);
1583 }
1584 
1585 
1586 
1587 //---------------------------------------------------------------------------
1588 
1589 template <typename MatrixType>
1590 inline void
1592  const MatrixType & rA,
1593  const typename BaseClass::AdditionalData &parameters)
1594 {
1595  this->PreconditionRelaxation<MatrixType>::initialize(rA, parameters);
1596 
1597  // in case we have a SparseMatrix class, we can extract information about
1598  // the diagonal.
1600  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
1601  &*this->A);
1602 
1603  // calculate the positions first after the diagonal.
1604  if (mat != nullptr)
1605  {
1606  const size_type n = this->A->n();
1607  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1608  for (size_type row = 0; row < n; ++row)
1609  {
1610  // find the first element in this line which is on the right of the
1611  // diagonal. we need to precondition with the elements on the left
1612  // only. note: the first entry in each line denotes the diagonal
1613  // element, which we need not check.
1615  it = mat->begin(row) + 1;
1616  for (; it < mat->end(row); ++it)
1617  if (it->column() > row)
1618  break;
1619  pos_right_of_diagonal[row] = it - mat->begin();
1620  }
1621  }
1622 }
1623 
1624 
1625 template <typename MatrixType>
1626 template <class VectorType>
1627 inline void
1629  const VectorType &src) const
1630 {
1631  static_assert(
1632  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1633  typename VectorType::size_type>::value,
1634  "PreconditionSSOR and VectorType must have the same size_type.");
1635 
1636  Assert(this->A != nullptr, ExcNotInitialized());
1637  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1638 }
1639 
1640 
1641 
1642 template <typename MatrixType>
1643 template <class VectorType>
1644 inline void
1646  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->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1655 }
1656 
1657 
1658 
1659 template <typename MatrixType>
1660 template <class VectorType>
1661 inline void
1663 {
1664  static_assert(
1665  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1666  typename VectorType::size_type>::value,
1667  "PreconditionSSOR and VectorType must have the same size_type.");
1668 
1669  Assert(this->A != nullptr, ExcNotInitialized());
1670  this->A->SSOR_step(dst, src, this->relaxation);
1671 }
1672 
1673 
1674 
1675 template <typename MatrixType>
1676 template <class VectorType>
1677 inline void
1679  const VectorType &src) const
1680 {
1681  static_assert(
1682  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1683  typename VectorType::size_type>::value,
1684  "PreconditionSSOR and VectorType must have the same size_type.");
1685 
1686  step(dst, src);
1687 }
1688 
1689 
1690 
1691 //---------------------------------------------------------------------------
1692 
1693 template <typename MatrixType>
1694 inline void
1696  const MatrixType & rA,
1697  const std::vector<size_type> & p,
1698  const std::vector<size_type> & ip,
1699  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1700 {
1701  permutation = &p;
1702  inverse_permutation = &ip;
1704 }
1705 
1706 
1707 template <typename MatrixType>
1708 inline void
1709 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
1710  const AdditionalData &additional_data)
1711 {
1712  initialize(A,
1713  additional_data.permutation,
1714  additional_data.inverse_permutation,
1715  additional_data.parameters);
1716 }
1717 
1718 
1719 template <typename MatrixType>
1720 template <typename VectorType>
1721 inline void
1723  const VectorType &src) const
1724 {
1725  static_assert(
1726  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1727  typename VectorType::size_type>::value,
1728  "PreconditionPSOR and VectorType must have the same size_type.");
1729 
1730  Assert(this->A != nullptr, ExcNotInitialized());
1731  dst = src;
1732  this->A->PSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1733 }
1734 
1735 
1736 
1737 template <typename MatrixType>
1738 template <class VectorType>
1739 inline void
1741  const VectorType &src) const
1742 {
1743  static_assert(
1744  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1745  typename VectorType::size_type>::value,
1746  "PreconditionPSOR and VectorType must have the same size_type.");
1747 
1748  Assert(this->A != nullptr, ExcNotInitialized());
1749  dst = src;
1750  this->A->TPSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1751 }
1752 
1753 template <typename MatrixType>
1755  const std::vector<size_type> &permutation,
1756  const std::vector<size_type> &inverse_permutation,
1757  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1758  : permutation(permutation)
1759  , inverse_permutation(inverse_permutation)
1760  , parameters(parameters)
1761 {}
1762 
1763 
1764 //---------------------------------------------------------------------------
1765 
1766 
1767 template <typename MatrixType, class VectorType>
1769  const MatrixType & M,
1770  const function_ptr method)
1771  : matrix(M)
1772  , precondition(method)
1773 {}
1774 
1775 
1776 
1777 template <typename MatrixType, class VectorType>
1778 void
1780  VectorType & dst,
1781  const VectorType &src) const
1782 {
1783  (matrix.*precondition)(dst, src);
1784 }
1785 
1786 //---------------------------------------------------------------------------
1787 
1788 template <typename MatrixType>
1790  const double relaxation)
1791  : relaxation(relaxation)
1792 {}
1793 
1794 
1795 
1796 //---------------------------------------------------------------------------
1797 
1798 namespace internal
1799 {
1800  namespace PreconditionChebyshevImplementation
1801  {
1802  // for deal.II vectors, perform updates for Chebyshev preconditioner all
1803  // at once to reduce memory transfer. Here, we select between general
1804  // vectors and deal.II vectors where we expand the loop over the (local)
1805  // size of the vector
1806 
1807  // generic part for non-deal.II vectors
1808  template <typename VectorType, typename PreconditionerType>
1809  inline void
1810  vector_updates(const VectorType & rhs,
1811  const PreconditionerType &preconditioner,
1812  const unsigned int iteration_index,
1813  const double factor1,
1814  const double factor2,
1818  VectorType & solution)
1819  {
1820  if (iteration_index == 0)
1821  {
1822  solution.equ(factor2, rhs);
1823  preconditioner.vmult(solution_old, solution);
1824  }
1825  else if (iteration_index == 1)
1826  {
1827  // compute t = P^{-1} * (b-A*x^{n})
1828  temp_vector1.sadd(-1.0, 1.0, rhs);
1829  preconditioner.vmult(solution_old, temp_vector1);
1830 
1831  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
1832  solution_old.sadd(factor2, 1 + factor1, solution);
1833  }
1834  else
1835  {
1836  // compute t = P^{-1} * (b-A*x^{n})
1837  temp_vector1.sadd(-1.0, 1.0, rhs);
1838  preconditioner.vmult(temp_vector2, temp_vector1);
1839 
1840  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
1841  solution_old.sadd(-factor1, factor2, temp_vector2);
1842  solution_old.add(1 + factor1, solution);
1843  }
1844 
1845  solution.swap(solution_old);
1846  }
1847 
1848  // worker routine for deal.II vectors. Because of vectorization, we need
1849  // to put the loop into an extra structure because the virtual function of
1850  // VectorUpdatesRange prevents the compiler from applying vectorization.
1851  template <typename Number>
1852  struct VectorUpdater
1853  {
1854  VectorUpdater(const Number * rhs,
1855  const Number * matrix_diagonal_inverse,
1856  const unsigned int iteration_index,
1857  const Number factor1,
1858  const Number factor2,
1859  Number * solution_old,
1860  Number * tmp_vector,
1861  Number * solution)
1862  : rhs(rhs)
1863  , matrix_diagonal_inverse(matrix_diagonal_inverse)
1864  , iteration_index(iteration_index)
1865  , factor1(factor1)
1866  , factor2(factor2)
1867  , solution_old(solution_old)
1868  , tmp_vector(tmp_vector)
1869  , solution(solution)
1870  {}
1871 
1872  void
1873  apply_to_subrange(const std::size_t begin, const std::size_t end) const
1874  {
1875  // To circumvent a bug in gcc
1876  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
1877  // copies of the variables factor1 and factor2 and do not check based on
1878  // factor1.
1879  const Number factor1 = this->factor1;
1880  const Number factor1_plus_1 = 1. + this->factor1;
1881  const Number factor2 = this->factor2;
1882  if (iteration_index == 0)
1883  {
1885  for (std::size_t i = begin; i < end; ++i)
1886  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
1887  }
1888  else if (iteration_index == 1)
1889  {
1890  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
1892  for (std::size_t i = begin; i < end; ++i)
1893  // for efficiency reason, write back to temp_vector that is
1894  // already read (avoid read-for-ownership)
1895  tmp_vector[i] =
1896  factor1_plus_1 * solution[i] +
1897  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1898  }
1899  else
1900  {
1901  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
1902  // + f_2 * P^{-1} * (b-A*x^{n})
1904  for (std::size_t i = begin; i < end; ++i)
1905  solution_old[i] =
1906  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
1907  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1908  }
1909  }
1910 
1911  const Number * rhs;
1912  const Number * matrix_diagonal_inverse;
1913  const unsigned int iteration_index;
1914  const Number factor1;
1915  const Number factor2;
1916  mutable Number * solution_old;
1917  mutable Number * tmp_vector;
1918  mutable Number * solution;
1919  };
1920 
1921  template <typename Number>
1922  struct VectorUpdatesRange : public ::parallel::ParallelForInteger
1923  {
1924  VectorUpdatesRange(const VectorUpdater<Number> &updater,
1925  const std::size_t size)
1926  : updater(updater)
1927  {
1929  VectorUpdatesRange::apply_to_subrange(0, size);
1930  else
1931  apply_parallel(
1932  0,
1933  size,
1935  }
1936 
1937  ~VectorUpdatesRange() override = default;
1938 
1939  virtual void
1940  apply_to_subrange(const std::size_t begin,
1941  const std::size_t end) const override
1942  {
1943  updater.apply_to_subrange(begin, end);
1944  }
1945 
1946  const VectorUpdater<Number> &updater;
1947  };
1948 
1949  // selection for diagonal matrix around deal.II vector
1950  template <typename Number>
1951  inline void
1952  vector_updates(const ::Vector<Number> & rhs,
1953  const DiagonalMatrix<::Vector<Number>> &jacobi,
1954  const unsigned int iteration_index,
1955  const double factor1,
1956  const double factor2,
1957  ::Vector<Number> &solution_old,
1958  ::Vector<Number> &temp_vector1,
1959  ::Vector<Number> &,
1960  ::Vector<Number> &solution)
1961  {
1962  VectorUpdater<Number> upd(rhs.begin(),
1963  jacobi.get_vector().begin(),
1964  iteration_index,
1965  factor1,
1966  factor2,
1967  solution_old.begin(),
1968  temp_vector1.begin(),
1969  solution.begin());
1970  VectorUpdatesRange<Number>(upd, rhs.size());
1971 
1972  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
1973  if (iteration_index == 0)
1974  {
1975  // nothing to do here because we can immediately write into the
1976  // solution vector without remembering any of the other vectors
1977  }
1978  else if (iteration_index == 1)
1979  {
1980  solution.swap(temp_vector1);
1981  solution_old.swap(temp_vector1);
1982  }
1983  else
1984  solution.swap(solution_old);
1985  }
1986 
1987  // selection for diagonal matrix around parallel deal.II vector
1988  template <typename Number>
1989  inline void
1990  vector_updates(
1992  const DiagonalMatrix<
1994  const unsigned int iteration_index,
1995  const double factor1,
1996  const double factor2,
1998  &solution_old,
2000  &temp_vector1,
2003  {
2004  VectorUpdater<Number> upd(rhs.begin(),
2005  jacobi.get_vector().begin(),
2006  iteration_index,
2007  factor1,
2008  factor2,
2009  solution_old.begin(),
2010  temp_vector1.begin(),
2011  solution.begin());
2012  VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
2013 
2014  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2015  if (iteration_index == 0)
2016  {
2017  // nothing to do here because we can immediately write into the
2018  // solution vector without remembering any of the other vectors
2019  }
2020  else if (iteration_index == 1)
2021  {
2022  solution.swap(temp_vector1);
2023  solution_old.swap(temp_vector1);
2024  }
2025  else
2026  solution.swap(solution_old);
2027  }
2028 
2029  template <typename MatrixType, typename PreconditionerType>
2030  inline void
2031  initialize_preconditioner(
2032  const MatrixType & matrix,
2033  std::shared_ptr<PreconditionerType> &preconditioner)
2034  {
2035  (void)matrix;
2036  (void)preconditioner;
2037  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2038  }
2039 
2040  template <typename MatrixType, typename VectorType>
2041  inline void
2042  initialize_preconditioner(
2043  const MatrixType & matrix,
2044  std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
2045  {
2046  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2047  {
2048  if (preconditioner.get() == nullptr)
2049  preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2050 
2051  Assert(
2052  preconditioner->m() == 0,
2053  ExcMessage(
2054  "Preconditioner appears to be initialized but not sized correctly"));
2055 
2056  // This part only works in serial
2057  if (preconditioner->m() != matrix.m())
2058  {
2059  preconditioner->get_vector().reinit(matrix.m());
2060  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2061  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2062  }
2063  }
2064  }
2065 
2066  template <typename VectorType>
2067  void
2068  set_initial_guess(VectorType &vector)
2069  {
2070  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2071  if (vector.locally_owned_elements().is_element(0))
2072  vector(0) = 0.;
2073  }
2074 
2075  template <typename Number>
2076  void
2077  set_initial_guess(::Vector<Number> &vector)
2078  {
2079  // Choose a high-frequency mode consisting of numbers between 0 and 1
2080  // that is cheap to compute (cheaper than random numbers) but avoids
2081  // obviously re-occurring numbers in multi-component systems by choosing
2082  // a period of 11
2083  for (unsigned int i = 0; i < vector.size(); ++i)
2084  vector(i) = i % 11;
2085 
2086  const Number mean_value = vector.mean_value();
2087  vector.add(-mean_value);
2088  }
2089 
2090  template <typename Number>
2091  void
2092  set_initial_guess(
2094  &vector)
2095  {
2096  // Choose a high-frequency mode consisting of numbers between 0 and 1
2097  // that is cheap to compute (cheaper than random numbers) but avoids
2098  // obviously re-occurring numbers in multi-component systems by choosing
2099  // a period of 11.
2100  // Make initial guess robust with respect to number of processors
2101  // by operating on the global index.
2102  types::global_dof_index first_local_range = 0;
2103  if (!vector.locally_owned_elements().is_empty())
2104  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2105  for (unsigned int i = 0; i < vector.locally_owned_size(); ++i)
2106  vector.local_element(i) = (i + first_local_range) % 11;
2107 
2108  const Number mean_value = vector.mean_value();
2109  vector.add(-mean_value);
2110  }
2111 
2112  template <typename Number>
2113  void
2114  set_initial_guess(
2116  {
2117  for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2118  set_initial_guess(vector.block(block));
2119  }
2120 
2121 
2122 # ifdef DEAL_II_COMPILER_CUDA_AWARE
2123  template <typename Number>
2124  __global__ void
2125  set_initial_guess_kernel(const types::global_dof_index offset,
2126  const unsigned int locally_owned_size,
2127  Number * values)
2128 
2129  {
2130  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
2131  if (index < locally_owned_size)
2132  values[index] = (index + offset) % 11;
2133  }
2134 
2135  template <typename Number>
2136  void
2137  set_initial_guess(
2139  &vector)
2140  {
2141  // Choose a high-frequency mode consisting of numbers between 0 and 1
2142  // that is cheap to compute (cheaper than random numbers) but avoids
2143  // obviously re-occurring numbers in multi-component systems by choosing
2144  // a period of 11.
2145  // Make initial guess robust with respect to number of processors
2146  // by operating on the global index.
2147  types::global_dof_index first_local_range = 0;
2148  if (!vector.locally_owned_elements().is_empty())
2149  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2150 
2151  const auto n_local_elements = vector.locally_owned_size();
2152  const int n_blocks =
2153  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
2154  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
2155  first_local_range, n_local_elements, vector.get_values());
2156  AssertCudaKernel();
2157 
2158  const Number mean_value = vector.mean_value();
2159  vector.add(-mean_value);
2160  }
2161 # endif // DEAL_II_COMPILER_CUDA_AWARE
2162 
2163  struct EigenvalueTracker
2164  {
2165  public:
2166  void
2167  slot(const std::vector<double> &eigenvalues)
2168  {
2169  values = eigenvalues;
2170  }
2171 
2172  std::vector<double> values;
2173  };
2174  } // namespace PreconditionChebyshevImplementation
2175 } // namespace internal
2176 
2177 
2178 
2179 template <typename MatrixType, class VectorType, typename PreconditionerType>
2181  AdditionalData::AdditionalData(const unsigned int degree,
2182  const double smoothing_range,
2183  const unsigned int eig_cg_n_iterations,
2184  const double eig_cg_residual,
2185  const double max_eigenvalue)
2186  : degree(degree)
2187  , smoothing_range(smoothing_range)
2188  , eig_cg_n_iterations(eig_cg_n_iterations)
2189  , eig_cg_residual(eig_cg_residual)
2190  , max_eigenvalue(max_eigenvalue)
2191 {}
2192 
2193 
2194 
2195 template <typename MatrixType, class VectorType, typename PreconditionerType>
2196 inline typename PreconditionChebyshev<MatrixType,
2197  VectorType,
2198  PreconditionerType>::AdditionalData &
2200  AdditionalData::operator=(const AdditionalData &other_data)
2201 {
2202  degree = other_data.degree;
2203  smoothing_range = other_data.smoothing_range;
2204  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
2205  eig_cg_residual = other_data.eig_cg_residual;
2206  max_eigenvalue = other_data.max_eigenvalue;
2207  preconditioner = other_data.preconditioner;
2208  constraints.copy_from(other_data.constraints);
2209 
2210  return *this;
2211 }
2212 
2213 
2214 
2215 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2218  : theta(1.)
2219  , delta(1.)
2221 {
2222  static_assert(
2223  std::is_same<size_type, typename VectorType::size_type>::value,
2224  "PreconditionChebyshev and VectorType must have the same size_type.");
2225 }
2226 
2227 
2228 
2229 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2230 inline void
2232  const MatrixType & matrix,
2233  const AdditionalData &additional_data)
2234 {
2235  matrix_ptr = &matrix;
2236  data = additional_data;
2237  Assert(data.degree > 0,
2238  ExcMessage("The degree of the Chebyshev method must be positive."));
2239  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2240  matrix, data.preconditioner);
2242 }
2243 
2244 
2245 
2246 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2247 inline void
2249 {
2251  theta = delta = 1.0;
2252  matrix_ptr = nullptr;
2253  {
2254  VectorType empty_vector;
2255  solution_old.reinit(empty_vector);
2256  temp_vector1.reinit(empty_vector);
2257  temp_vector2.reinit(empty_vector);
2258  }
2259  data.preconditioner.reset();
2260 }
2261 
2262 
2263 
2264 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2265 inline typename PreconditionChebyshev<MatrixType,
2266  VectorType,
2267  PreconditionerType>::EigenvalueInformation
2269  estimate_eigenvalues(const VectorType &src) const
2270 {
2272  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2273 
2275  EigenvalueInformation info{};
2276 
2277  solution_old.reinit(src);
2278  temp_vector1.reinit(src, true);
2279 
2280  if (data.eig_cg_n_iterations > 0)
2281  {
2283  ExcMessage(
2284  "Need to set at least two iterations to find eigenvalues."));
2285 
2286  // set a very strict tolerance to force at least two iterations
2287  ReductionControl control(
2289  std::sqrt(
2291  1e-10,
2292  false,
2293  false);
2294 
2295  internal::PreconditionChebyshevImplementation::EigenvalueTracker
2296  eigenvalue_tracker;
2297  SolverCG<VectorType> solver(control);
2298  solver.connect_eigenvalues_slot(
2299  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2300  eigenvalue_tracker.slot(eigenvalues);
2301  });
2302 
2303  // set an initial guess that contains some high-frequency parts (to the
2304  // extent possible without knowing the discretization and the numbering)
2305  // to trigger high eigenvalues according to the external function
2306  internal::PreconditionChebyshevImplementation::set_initial_guess(
2307  temp_vector1);
2309 
2310  try
2311  {
2312  solver.solve(*matrix_ptr,
2313  solution_old,
2314  temp_vector1,
2315  *data.preconditioner);
2316  }
2318  {}
2319 
2320  // read the eigenvalues from the attached eigenvalue tracker
2321  if (eigenvalue_tracker.values.empty())
2322  info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2323  else
2324  {
2325  info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2326 
2327  // include a safety factor since the CG method will in general not
2328  // be converged
2329  info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
2330  }
2331 
2332  info.cg_iterations = control.last_step();
2333  }
2334  else
2335  {
2336  info.max_eigenvalue_estimate = data.max_eigenvalue;
2337  info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
2338  }
2339 
2340  const double alpha = (data.smoothing_range > 1. ?
2341  info.max_eigenvalue_estimate / data.smoothing_range :
2342  std::min(0.9 * info.max_eigenvalue_estimate,
2343  info.min_eigenvalue_estimate));
2344 
2345  // in case the user set the degree to invalid unsigned int, we have to
2346  // determine the number of necessary iterations from the Chebyshev error
2347  // estimate, given the target tolerance specified by smoothing_range. This
2348  // estimate is based on the error formula given in section 5.1 of
2349  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
2351  {
2352  const double actual_range = info.max_eigenvalue_estimate / alpha;
2353  const double sigma = (1. - std::sqrt(1. / actual_range)) /
2354  (1. + std::sqrt(1. / actual_range));
2355  const double eps = data.smoothing_range;
2356  const_cast<
2358  this)
2359  ->data.degree =
2360  1 + static_cast<unsigned int>(
2361  std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
2362  std::log(1. / sigma));
2363  }
2364 
2365  info.degree = data.degree;
2366 
2367  const_cast<
2369  ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
2370  const_cast<
2372  ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
2373 
2374  // We do not need the second temporary vector in case we have a
2375  // DiagonalMatrix as preconditioner and use deal.II's own vectors
2376  using NumberType = typename VectorType::value_type;
2377  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
2378  false ||
2379  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
2380  ((std::is_same<VectorType,
2382  Vector<NumberType, MemorySpace::Host>>::value ==
2383  false) ||
2384  (std::is_same<VectorType,
2385  LinearAlgebra::distributed::
2386  Vector<NumberType, MemorySpace::CUDA>>::value ==
2387  false))))
2388  temp_vector2.reinit(src, true);
2389  else
2390  {
2391  VectorType empty_vector;
2392  temp_vector2.reinit(empty_vector);
2393  }
2394 
2395  const_cast<
2397  ->eigenvalues_are_initialized = true;
2398 
2399  return info;
2400 }
2401 
2402 
2403 
2404 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2405 inline void
2407  VectorType & solution,
2408  const VectorType &rhs) const
2409 {
2410  std::lock_guard<std::mutex> lock(mutex);
2411  if (eigenvalues_are_initialized == false)
2412  estimate_eigenvalues(rhs);
2413 
2414  internal::PreconditionChebyshevImplementation::vector_updates(
2415  rhs,
2417  0,
2418  0.,
2419  1. / theta,
2420  solution_old,
2421  temp_vector1,
2422  temp_vector2,
2423  solution);
2424 
2425  // if delta is zero, we do not need to iterate because the updates will be
2426  // zero
2427  if (data.degree < 2 || std::abs(delta) < 1e-40)
2428  return;
2429 
2430  double rhok = delta / theta, sigma = theta / delta;
2431  for (unsigned int k = 0; k < data.degree - 1; ++k)
2432  {
2433  matrix_ptr->vmult(temp_vector1, solution);
2434  const double rhokp = 1. / (2. * sigma - rhok);
2435  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2436  rhok = rhokp;
2437  internal::PreconditionChebyshevImplementation::vector_updates(
2438  rhs,
2440  k + 1,
2441  factor1,
2442  factor2,
2443  solution_old,
2444  temp_vector1,
2445  temp_vector2,
2446  solution);
2447  }
2448 }
2449 
2450 
2451 
2452 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2453 inline void
2455  VectorType & solution,
2456  const VectorType &rhs) const
2457 {
2458  std::lock_guard<std::mutex> lock(mutex);
2459  if (eigenvalues_are_initialized == false)
2460  estimate_eigenvalues(rhs);
2461 
2462  internal::PreconditionChebyshevImplementation::vector_updates(
2463  rhs,
2465  0,
2466  0.,
2467  1. / theta,
2468  solution_old,
2469  temp_vector1,
2470  temp_vector2,
2471  solution);
2472 
2473  if (data.degree < 2 || std::abs(delta) < 1e-40)
2474  return;
2475 
2476  double rhok = delta / theta, sigma = theta / delta;
2477  for (unsigned int k = 0; k < data.degree - 1; ++k)
2478  {
2479  matrix_ptr->Tvmult(temp_vector1, solution);
2480  const double rhokp = 1. / (2. * sigma - rhok);
2481  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2482  rhok = rhokp;
2483  internal::PreconditionChebyshevImplementation::vector_updates(
2484  rhs,
2486  k + 1,
2487  factor1,
2488  factor2,
2489  solution_old,
2490  temp_vector1,
2491  temp_vector2,
2492  solution);
2493  }
2494 }
2495 
2496 
2497 
2498 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2499 inline void
2501  VectorType & solution,
2502  const VectorType &rhs) const
2503 {
2504  std::lock_guard<std::mutex> lock(mutex);
2505  if (eigenvalues_are_initialized == false)
2506  estimate_eigenvalues(rhs);
2507 
2508  matrix_ptr->vmult(temp_vector1, solution);
2509  internal::PreconditionChebyshevImplementation::vector_updates(
2510  rhs,
2512  1,
2513  0.,
2514  1. / theta,
2515  solution_old,
2516  temp_vector1,
2517  temp_vector2,
2518  solution);
2519 
2520  if (data.degree < 2 || std::abs(delta) < 1e-40)
2521  return;
2522 
2523  double rhok = delta / theta, sigma = theta / delta;
2524  for (unsigned int k = 0; k < data.degree - 1; ++k)
2525  {
2526  matrix_ptr->vmult(temp_vector1, solution);
2527  const double rhokp = 1. / (2. * sigma - rhok);
2528  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2529  rhok = rhokp;
2530  internal::PreconditionChebyshevImplementation::vector_updates(
2531  rhs,
2533  k + 2,
2534  factor1,
2535  factor2,
2536  solution_old,
2537  temp_vector1,
2538  temp_vector2,
2539  solution);
2540  }
2541 }
2542 
2543 
2544 
2545 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2546 inline void
2548  VectorType & solution,
2549  const VectorType &rhs) const
2550 {
2551  std::lock_guard<std::mutex> lock(mutex);
2552  if (eigenvalues_are_initialized == false)
2553  estimate_eigenvalues(rhs);
2554 
2555  matrix_ptr->Tvmult(temp_vector1, solution);
2556  internal::PreconditionChebyshevImplementation::vector_updates(
2557  rhs,
2559  1,
2560  0.,
2561  1. / theta,
2562  solution_old,
2563  temp_vector1,
2564  temp_vector2,
2565  solution);
2566 
2567  if (data.degree < 2 || std::abs(delta) < 1e-40)
2568  return;
2569 
2570  double rhok = delta / theta, sigma = theta / delta;
2571  for (unsigned int k = 0; k < data.degree - 1; ++k)
2572  {
2573  matrix_ptr->Tvmult(temp_vector1, solution);
2574  const double rhokp = 1. / (2. * sigma - rhok);
2575  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2576  rhok = rhokp;
2577  internal::PreconditionChebyshevImplementation::vector_updates(
2578  rhs,
2580  k + 2,
2581  factor1,
2582  factor2,
2583  solution_old,
2584  temp_vector1,
2585  temp_vector2,
2586  solution);
2587  }
2588 }
2589 
2590 
2591 
2592 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2593 inline typename PreconditionChebyshev<MatrixType,
2594  VectorType,
2595  PreconditionerType>::size_type
2597 {
2598  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2599  return matrix_ptr->m();
2600 }
2601 
2602 
2603 
2604 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2605 inline typename PreconditionChebyshev<MatrixType,
2606  VectorType,
2607  PreconditionerType>::size_type
2609 {
2610  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2611  return matrix_ptr->n();
2612 }
2613 
2614 #endif // DOXYGEN
2615 
2617 
2618 #endif
void Tvmult(VectorType &, const VectorType &) const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
void set_zero(VectorType &vec) const
size_type m() const
void initialize(const MatrixType &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
const_iterator end() const
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void Tstep(VectorType &x, const VectorType &rhs) const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void vmult(VectorType &dst, const VectorType &src) const
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionRelaxation< MatrixType >::AdditionalData parameters
Definition: precondition.h:791
void vmult(VectorType &, const VectorType &) const
size_type m() const
constexpr int block_size
Definition: cuda_size.h:29
const MatrixType & matrix
Definition: precondition.h:386
AdditionalData data
AdditionalData(const double relaxation=1.)
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition: precondition.h:366
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:1575
AdditionalData(const unsigned int degree=1, const double smoothing_range=0., const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1)
AdditionalData(const double relaxation=1.)
#define AssertCudaKernel()
Definition: exceptions.h:1827
size_type m() const
size_type locally_owned_size() const
const std::vector< size_type > * permutation
Definition: precondition.h:844
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
const function_ptr precondition
Definition: precondition.h:391
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:408
AdditionalData & operator=(const AdditionalData &other_data)
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:787
#define Assert(cond, exc)
Definition: exceptions.h:1465
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:657
void step(VectorType &x, const VectorType &rhs) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
VectorType::value_type * end(VectorType &V)
typename MatrixType::size_type size_type
Definition: precondition.h:755
void vmult_add(VectorType &, const VectorType &) const
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
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 Tstep(VectorType &x, const VectorType &rhs) const
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
void vmult(VectorType &, const VectorType &) const
void Tstep(VectorType &x, const VectorType &rhs) const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
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
static const char A
void vmult(VectorType &, const VectorType &) const
unsigned int global_dof_index
Definition: types.h:76
virtual ::IndexSet locally_owned_elements() const override
size_type n() const
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
Definition: precondition.h:460
void Tvmult(VectorType &, const VectorType &) const
void step(VectorType &x, const VectorType &rhs) const
AffineConstraints< double > constraints
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< double > &)> &slot, const bool every_iteration=false)
size_type n() const
void swap(Vector< Number, MemorySpace > &v)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
VectorType::value_type * begin(VectorType &V)
typename MatrixType::size_type size_type
Definition: precondition.h:662
const std::vector< size_type > * inverse_permutation
Definition: precondition.h:848
void Tstep(VectorType &dst, const VectorType &src) const
void step(VectorType &dst, const VectorType &src) const
void step(VectorType &x, const VectorType &rhs) const
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
virtual void add(const Number a) override
size_type m() const
const_iterator begin() const
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:504
void vmult(VectorType &, const VectorType &) const
const std::vector< size_type > & permutation
Definition: precondition.h:783
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:590
Eigenvalue vector is filled.
virtual Number mean_value() const override
void vmult_add(VectorType &, const VectorType &) const
BlockType & block(const unsigned int i)
std::vector< std::size_t > pos_right_of_diagonal
Definition: precondition.h:715
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1880
bool is_empty() const
Definition: index_set.h:1824
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:135
SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
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())