Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 - 2020 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  } // namespace distributed
51 } // namespace LinearAlgebra
52 #endif
53 
54 
81 {
82 public:
87 
93  {
97  AdditionalData() = default;
98  };
99 
104 
109  template <typename MatrixType>
110  void
111  initialize(const MatrixType & matrix,
112  const AdditionalData &additional_data = AdditionalData());
113 
117  template <class VectorType>
118  void
119  vmult(VectorType &, const VectorType &) const;
120 
125  template <class VectorType>
126  void
127  Tvmult(VectorType &, const VectorType &) const;
128 
132  template <class VectorType>
133  void
134  vmult_add(VectorType &, const VectorType &) const;
135 
140  template <class VectorType>
141  void
142  Tvmult_add(VectorType &, const VectorType &) const;
143 
148  void
150  {}
151 
159  size_type
160  m() const;
161 
169  size_type
170  n() const;
171 
172 private:
177 
182 };
183 
184 
185 
200 {
201 public:
206 
211  {
212  public:
217  AdditionalData(const double relaxation = 1.);
218 
222  double relaxation;
223  };
224 
230 
234  void
235  initialize(const AdditionalData &parameters);
236 
242  template <typename MatrixType>
243  void
244  initialize(const MatrixType &matrix, const AdditionalData &parameters);
245 
249  template <class VectorType>
250  void
251  vmult(VectorType &, const VectorType &) const;
252 
257  template <class VectorType>
258  void
259  Tvmult(VectorType &, const VectorType &) const;
263  template <class VectorType>
264  void
265  vmult_add(VectorType &, const VectorType &) const;
266 
271  template <class VectorType>
272  void
273  Tvmult_add(VectorType &, const VectorType &) const;
274 
279  void
281  {}
282 
290  size_type
291  m() const;
292 
300  size_type
301  n() const;
302 
303 private:
307  double relaxation;
308 
313 
318 };
319 
320 
321 
363 template <typename MatrixType = SparseMatrix<double>,
364  class VectorType = Vector<double>>
366 {
367 public:
371  using function_ptr = void (MatrixType::*)(VectorType &,
372  const VectorType &) const;
373 
379  PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
380 
385  void
386  vmult(VectorType &dst, const VectorType &src) const;
387 
388 private:
392  const MatrixType &matrix;
393 
398 };
399 
400 
401 
410 template <typename MatrixType = SparseMatrix<double>>
412 {
413 public:
418 
423  {
424  public:
428  AdditionalData(const double relaxation = 1.);
429 
433  double relaxation;
434  };
435 
441  void
442  initialize(const MatrixType & A,
443  const AdditionalData &parameters = AdditionalData());
444 
448  void
449  clear();
450 
455  size_type
456  m() const;
457 
462  size_type
463  n() const;
464 
465 protected:
470 
474  double relaxation;
475 };
476 
477 
478 
507 template <typename MatrixType = SparseMatrix<double>>
508 class PreconditionJacobi : public PreconditionRelaxation<MatrixType>
509 {
510 public:
514  using AdditionalData =
516 
520  template <class VectorType>
521  void
522  vmult(VectorType &, const VectorType &) const;
523 
528  template <class VectorType>
529  void
530  Tvmult(VectorType &, const VectorType &) const;
531 
535  template <class VectorType>
536  void
537  step(VectorType &x, const VectorType &rhs) const;
538 
542  template <class VectorType>
543  void
544  Tstep(VectorType &x, const VectorType &rhs) const;
545 };
546 
547 
595 template <typename MatrixType = SparseMatrix<double>>
596 class PreconditionSOR : public PreconditionRelaxation<MatrixType>
597 {
598 public:
602  using AdditionalData =
604 
608  template <class VectorType>
609  void
610  vmult(VectorType &, const VectorType &) const;
611 
615  template <class VectorType>
616  void
617  Tvmult(VectorType &, const VectorType &) const;
618 
622  template <class VectorType>
623  void
624  step(VectorType &x, const VectorType &rhs) const;
625 
629  template <class VectorType>
630  void
631  Tstep(VectorType &x, const VectorType &rhs) const;
632 };
633 
634 
635 
664 template <typename MatrixType = SparseMatrix<double>>
665 class PreconditionSSOR : public PreconditionRelaxation<MatrixType>
666 {
667 public:
671  using AdditionalData =
673 
678 
683 
684 
690  void
691  initialize(const MatrixType & A,
692  const typename BaseClass::AdditionalData &parameters =
693  typename BaseClass::AdditionalData());
694 
698  template <class VectorType>
699  void
700  vmult(VectorType &, const VectorType &) const;
701 
706  template <class VectorType>
707  void
708  Tvmult(VectorType &, const VectorType &) const;
709 
710 
714  template <class VectorType>
715  void
716  step(VectorType &x, const VectorType &rhs) const;
717 
721  template <class VectorType>
722  void
723  Tstep(VectorType &x, const VectorType &rhs) const;
724 
725 private:
730  std::vector<std::size_t> pos_right_of_diagonal;
731 };
732 
733 
766 template <typename MatrixType = SparseMatrix<double>>
767 class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
768 {
769 public:
774 
779  {
780  public:
792  const std::vector<size_type> &permutation,
793  const std::vector<size_type> &inverse_permutation,
795  &parameters =
797 
801  const std::vector<size_type> &permutation;
805  const std::vector<size_type> &inverse_permutation;
810  };
811 
823  void
824  initialize(const MatrixType & A,
825  const std::vector<size_type> &permutation,
826  const std::vector<size_type> &inverse_permutation,
828  &parameters =
830 
841  void
842  initialize(const MatrixType &A, const AdditionalData &additional_data);
843 
847  template <class VectorType>
848  void
849  vmult(VectorType &, const VectorType &) const;
850 
854  template <class VectorType>
855  void
856  Tvmult(VectorType &, const VectorType &) const;
857 
858 private:
862  const std::vector<size_type> *permutation;
866  const std::vector<size_type> *inverse_permutation;
867 };
868 
869 
870 
1009 template <typename MatrixType = SparseMatrix<double>,
1010  typename VectorType = Vector<double>,
1011  typename PreconditionerType = DiagonalMatrix<VectorType>>
1013 {
1014 public:
1019 
1025  {
1029  AdditionalData(const unsigned int degree = 1,
1030  const double smoothing_range = 0.,
1031  const unsigned int eig_cg_n_iterations = 8,
1032  const double eig_cg_residual = 1e-2,
1033  const double max_eigenvalue = 1);
1034 
1038  AdditionalData &
1039  operator=(const AdditionalData &other_data);
1040 
1052  unsigned int degree;
1053 
1066 
1073  unsigned int eig_cg_n_iterations;
1074 
1080 
1087 
1093 
1097  std::shared_ptr<PreconditionerType> preconditioner;
1098  };
1099 
1100 
1105 
1117  void
1118  initialize(const MatrixType & matrix,
1119  const AdditionalData &additional_data = AdditionalData());
1120 
1125  void
1126  vmult(VectorType &dst, const VectorType &src) const;
1127 
1132  void
1133  Tvmult(VectorType &dst, const VectorType &src) const;
1134 
1138  void
1139  step(VectorType &dst, const VectorType &src) const;
1140 
1144  void
1145  Tstep(VectorType &dst, const VectorType &src) const;
1146 
1150  void
1151  clear();
1152 
1157  size_type
1158  m() const;
1159 
1164  size_type
1165  n() const;
1166 
1172  {
1184  unsigned int cg_iterations;
1189  unsigned int degree;
1194  : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
1195  , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
1196  , cg_iterations{0}
1197  , degree{0}
1198  {}
1199  };
1200 
1214  estimate_eigenvalues(const VectorType &src) const;
1215 
1216 private:
1220  SmartPointer<
1221  const MatrixType,
1224 
1229 
1234 
1239 
1245 
1249  double theta;
1250 
1255  double delta;
1256 
1262 
1268 };
1269 
1270 
1271 
1273 /* ---------------------------------- Inline functions ------------------- */
1274 
1275 #ifndef DOXYGEN
1276 
1278  : n_rows(0)
1279  , n_columns(0)
1280 {}
1281 
1282 template <typename MatrixType>
1283 inline void
1284 PreconditionIdentity::initialize(const MatrixType &matrix,
1286 {
1287  n_rows = matrix.m();
1288  n_columns = matrix.n();
1289 }
1290 
1291 
1292 template <class VectorType>
1293 inline void
1294 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1295 {
1296  dst = src;
1297 }
1298 
1299 
1300 
1301 template <class VectorType>
1302 inline void
1303 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1304 {
1305  dst = src;
1306 }
1307 
1308 template <class VectorType>
1309 inline void
1311 {
1312  dst += src;
1313 }
1314 
1315 
1316 
1317 template <class VectorType>
1318 inline void
1320 {
1321  dst += src;
1322 }
1323 
1326 {
1327  Assert(n_rows != 0, ExcNotInitialized());
1328  return n_rows;
1329 }
1330 
1333 {
1334  Assert(n_columns != 0, ExcNotInitialized());
1335  return n_columns;
1336 }
1337 
1338 //---------------------------------------------------------------------------
1339 
1341  const double relaxation)
1342  : relaxation(relaxation)
1343 {}
1344 
1345 
1347  : relaxation(0)
1348  , n_rows(0)
1349  , n_columns(0)
1350 {
1351  AdditionalData add_data;
1352  relaxation = add_data.relaxation;
1353 }
1354 
1355 
1356 
1357 inline void
1359  const PreconditionRichardson::AdditionalData &parameters)
1360 {
1361  relaxation = parameters.relaxation;
1362 }
1363 
1364 
1365 
1366 template <typename MatrixType>
1367 inline void
1369  const MatrixType & matrix,
1370  const PreconditionRichardson::AdditionalData &parameters)
1371 {
1372  relaxation = parameters.relaxation;
1373  n_rows = matrix.m();
1374  n_columns = matrix.n();
1375 }
1376 
1377 
1378 
1379 template <class VectorType>
1380 inline void
1381 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
1382 {
1383  static_assert(
1385  "PreconditionRichardson and VectorType must have the same size_type.");
1386 
1387  dst.equ(relaxation, src);
1388 }
1389 
1390 
1391 
1392 template <class VectorType>
1393 inline void
1395 {
1396  static_assert(
1398  "PreconditionRichardson and VectorType must have the same size_type.");
1399 
1400  dst.equ(relaxation, src);
1401 }
1402 
1403 template <class VectorType>
1404 inline void
1406 {
1407  static_assert(
1409  "PreconditionRichardson and VectorType must have the same size_type.");
1410 
1411  dst.add(relaxation, src);
1412 }
1413 
1414 
1415 
1416 template <class VectorType>
1417 inline void
1419 {
1420  static_assert(
1422  "PreconditionRichardson and VectorType must have the same size_type.");
1423 
1424  dst.add(relaxation, src);
1425 }
1426 
1429 {
1430  Assert(n_rows != 0, ExcNotInitialized());
1431  return n_rows;
1432 }
1433 
1436 {
1437  Assert(n_columns != 0, ExcNotInitialized());
1438  return n_columns;
1439 }
1440 
1441 //---------------------------------------------------------------------------
1442 
1443 template <typename MatrixType>
1444 inline void
1446  const AdditionalData &parameters)
1447 {
1448  A = &rA;
1449  relaxation = parameters.relaxation;
1450 }
1451 
1452 
1453 template <typename MatrixType>
1454 inline void
1456 {
1457  A = nullptr;
1458 }
1459 
1460 template <typename MatrixType>
1463 {
1464  Assert(A != nullptr, ExcNotInitialized());
1465  return A->m();
1466 }
1467 
1468 template <typename MatrixType>
1471 {
1472  Assert(A != nullptr, ExcNotInitialized());
1473  return A->n();
1474 }
1475 
1476 //---------------------------------------------------------------------------
1477 
1478 template <typename MatrixType>
1479 template <class VectorType>
1480 inline void
1482  const VectorType &src) const
1483 {
1484  static_assert(
1485  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1486  typename VectorType::size_type>::value,
1487  "PreconditionJacobi and VectorType must have the same size_type.");
1488 
1489  Assert(this->A != nullptr, ExcNotInitialized());
1490  this->A->precondition_Jacobi(dst, src, this->relaxation);
1491 }
1492 
1493 
1494 
1495 template <typename MatrixType>
1496 template <class VectorType>
1497 inline void
1499  const VectorType &src) const
1500 {
1501  static_assert(
1502  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1503  typename VectorType::size_type>::value,
1504  "PreconditionJacobi and VectorType must have the same size_type.");
1505 
1506  Assert(this->A != nullptr, ExcNotInitialized());
1507  this->A->precondition_Jacobi(dst, src, this->relaxation);
1508 }
1509 
1510 
1511 
1512 template <typename MatrixType>
1513 template <class VectorType>
1514 inline void
1516  const VectorType &src) const
1517 {
1518  static_assert(
1519  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1520  typename VectorType::size_type>::value,
1521  "PreconditionJacobi and VectorType must have the same size_type.");
1522 
1523  Assert(this->A != nullptr, ExcNotInitialized());
1524  this->A->Jacobi_step(dst, src, this->relaxation);
1525 }
1526 
1527 
1528 
1529 template <typename MatrixType>
1530 template <class VectorType>
1531 inline void
1533  const VectorType &src) const
1534 {
1535  static_assert(
1536  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1537  typename VectorType::size_type>::value,
1538  "PreconditionJacobi and VectorType must have the same size_type.");
1539 
1540  step(dst, src);
1541 }
1542 
1543 
1544 
1545 //---------------------------------------------------------------------------
1546 
1547 template <typename MatrixType>
1548 template <class VectorType>
1549 inline void
1551 {
1552  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1553  typename VectorType::size_type>::value,
1554  "PreconditionSOR and VectorType must have the same size_type.");
1555 
1556  Assert(this->A != nullptr, ExcNotInitialized());
1557  this->A->precondition_SOR(dst, src, this->relaxation);
1558 }
1559 
1560 
1561 
1562 template <typename MatrixType>
1563 template <class VectorType>
1564 inline void
1566  const VectorType &src) const
1567 {
1568  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1569  typename VectorType::size_type>::value,
1570  "PreconditionSOR and VectorType must have the same size_type.");
1571 
1572  Assert(this->A != nullptr, ExcNotInitialized());
1573  this->A->precondition_TSOR(dst, src, this->relaxation);
1574 }
1575 
1576 
1577 
1578 template <typename MatrixType>
1579 template <class VectorType>
1580 inline void
1582 {
1583  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1584  typename VectorType::size_type>::value,
1585  "PreconditionSOR and VectorType must have the same size_type.");
1586 
1587  Assert(this->A != nullptr, ExcNotInitialized());
1588  this->A->SOR_step(dst, src, this->relaxation);
1589 }
1590 
1591 
1592 
1593 template <typename MatrixType>
1594 template <class VectorType>
1595 inline void
1597 {
1598  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1599  typename VectorType::size_type>::value,
1600  "PreconditionSOR and VectorType must have the same size_type.");
1601 
1602  Assert(this->A != nullptr, ExcNotInitialized());
1603  this->A->TSOR_step(dst, src, this->relaxation);
1604 }
1605 
1606 
1607 
1608 //---------------------------------------------------------------------------
1609 
1610 template <typename MatrixType>
1611 inline void
1613  const MatrixType & rA,
1614  const typename BaseClass::AdditionalData &parameters)
1615 {
1616  this->PreconditionRelaxation<MatrixType>::initialize(rA, parameters);
1617 
1618  // in case we have a SparseMatrix class, we can extract information about
1619  // the diagonal.
1621  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
1622  &*this->A);
1623 
1624  // calculate the positions first after the diagonal.
1625  if (mat != nullptr)
1626  {
1627  const size_type n = this->A->n();
1628  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1629  for (size_type row = 0; row < n; ++row)
1630  {
1631  // find the first element in this line which is on the right of the
1632  // diagonal. we need to precondition with the elements on the left
1633  // only. note: the first entry in each line denotes the diagonal
1634  // element, which we need not check.
1636  it = mat->begin(row) + 1;
1637  for (; it < mat->end(row); ++it)
1638  if (it->column() > row)
1639  break;
1640  pos_right_of_diagonal[row] = it - mat->begin();
1641  }
1642  }
1643 }
1644 
1645 
1646 template <typename MatrixType>
1647 template <class VectorType>
1648 inline void
1650  const VectorType &src) const
1651 {
1652  static_assert(
1653  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1654  typename VectorType::size_type>::value,
1655  "PreconditionSSOR and VectorType must have the same size_type.");
1656 
1657  Assert(this->A != nullptr, ExcNotInitialized());
1658  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1659 }
1660 
1661 
1662 
1663 template <typename MatrixType>
1664 template <class VectorType>
1665 inline void
1667  const VectorType &src) const
1668 {
1669  static_assert(
1670  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1671  typename VectorType::size_type>::value,
1672  "PreconditionSSOR and VectorType must have the same size_type.");
1673 
1674  Assert(this->A != nullptr, ExcNotInitialized());
1675  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1676 }
1677 
1678 
1679 
1680 template <typename MatrixType>
1681 template <class VectorType>
1682 inline void
1684 {
1685  static_assert(
1686  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1687  typename VectorType::size_type>::value,
1688  "PreconditionSSOR and VectorType must have the same size_type.");
1689 
1690  Assert(this->A != nullptr, ExcNotInitialized());
1691  this->A->SSOR_step(dst, src, this->relaxation);
1692 }
1693 
1694 
1695 
1696 template <typename MatrixType>
1697 template <class VectorType>
1698 inline void
1700  const VectorType &src) const
1701 {
1702  static_assert(
1703  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1704  typename VectorType::size_type>::value,
1705  "PreconditionSSOR and VectorType must have the same size_type.");
1706 
1707  step(dst, src);
1708 }
1709 
1710 
1711 
1712 //---------------------------------------------------------------------------
1713 
1714 template <typename MatrixType>
1715 inline void
1717  const MatrixType & rA,
1718  const std::vector<size_type> & p,
1719  const std::vector<size_type> & ip,
1720  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1721 {
1722  permutation = &p;
1723  inverse_permutation = &ip;
1725 }
1726 
1727 
1728 template <typename MatrixType>
1729 inline void
1730 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
1731  const AdditionalData &additional_data)
1732 {
1733  initialize(A,
1734  additional_data.permutation,
1735  additional_data.inverse_permutation,
1736  additional_data.parameters);
1737 }
1738 
1739 
1740 template <typename MatrixType>
1741 template <typename VectorType>
1742 inline void
1744  const VectorType &src) const
1745 {
1746  static_assert(
1747  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1748  typename VectorType::size_type>::value,
1749  "PreconditionPSOR and VectorType must have the same size_type.");
1750 
1751  Assert(this->A != nullptr, ExcNotInitialized());
1752  dst = src;
1753  this->A->PSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1754 }
1755 
1756 
1757 
1758 template <typename MatrixType>
1759 template <class VectorType>
1760 inline void
1762  const VectorType &src) const
1763 {
1764  static_assert(
1765  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1766  typename VectorType::size_type>::value,
1767  "PreconditionPSOR and VectorType must have the same size_type.");
1768 
1769  Assert(this->A != nullptr, ExcNotInitialized());
1770  dst = src;
1771  this->A->TPSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1772 }
1773 
1774 template <typename MatrixType>
1776  const std::vector<size_type> &permutation,
1777  const std::vector<size_type> &inverse_permutation,
1778  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1779  : permutation(permutation)
1780  , inverse_permutation(inverse_permutation)
1781  , parameters(parameters)
1782 {}
1783 
1784 
1785 //---------------------------------------------------------------------------
1786 
1787 
1788 template <typename MatrixType, class VectorType>
1790  const MatrixType & M,
1791  const function_ptr method)
1792  : matrix(M)
1793  , precondition(method)
1794 {}
1795 
1796 
1797 
1798 template <typename MatrixType, class VectorType>
1799 void
1801  VectorType & dst,
1802  const VectorType &src) const
1803 {
1804  (matrix.*precondition)(dst, src);
1805 }
1806 
1807 //---------------------------------------------------------------------------
1808 
1809 template <typename MatrixType>
1811  const double relaxation)
1812  : relaxation(relaxation)
1813 {}
1814 
1815 
1816 
1817 //---------------------------------------------------------------------------
1818 
1819 namespace internal
1820 {
1821  namespace PreconditionChebyshevImplementation
1822  {
1823  // for deal.II vectors, perform updates for Chebyshev preconditioner all
1824  // at once to reduce memory transfer. Here, we select between general
1825  // vectors and deal.II vectors where we expand the loop over the (local)
1826  // size of the vector
1827 
1828  // generic part for non-deal.II vectors
1829  template <typename VectorType, typename PreconditionerType>
1830  inline void
1831  vector_updates(const VectorType & rhs,
1832  const PreconditionerType &preconditioner,
1833  const unsigned int iteration_index,
1834  const double factor1,
1835  const double factor2,
1839  VectorType & solution)
1840  {
1841  if (iteration_index == 0)
1842  {
1843  solution.equ(factor2, rhs);
1844  preconditioner.vmult(solution_old, solution);
1845  }
1846  else if (iteration_index == 1)
1847  {
1848  // compute t = P^{-1} * (b-A*x^{n})
1849  temp_vector1.sadd(-1.0, 1.0, rhs);
1850  preconditioner.vmult(solution_old, temp_vector1);
1851 
1852  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
1853  solution_old.sadd(factor2, 1 + factor1, solution);
1854  }
1855  else
1856  {
1857  // compute t = P^{-1} * (b-A*x^{n})
1858  temp_vector1.sadd(-1.0, 1.0, rhs);
1859  preconditioner.vmult(temp_vector2, temp_vector1);
1860 
1861  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
1862  solution_old.sadd(-factor1, factor2, temp_vector2);
1863  solution_old.add(1 + factor1, solution);
1864  }
1865 
1866  solution.swap(solution_old);
1867  }
1868 
1869  // worker routine for deal.II vectors. Because of vectorization, we need
1870  // to put the loop into an extra structure because the virtual function of
1871  // VectorUpdatesRange prevents the compiler from applying vectorization.
1872  template <typename Number>
1873  struct VectorUpdater
1874  {
1875  VectorUpdater(const Number * rhs,
1876  const Number * matrix_diagonal_inverse,
1877  const unsigned int iteration_index,
1878  const Number factor1,
1879  const Number factor2,
1880  Number * solution_old,
1881  Number * tmp_vector,
1882  Number * solution)
1883  : rhs(rhs)
1884  , matrix_diagonal_inverse(matrix_diagonal_inverse)
1885  , iteration_index(iteration_index)
1886  , factor1(factor1)
1887  , factor2(factor2)
1888  , solution_old(solution_old)
1889  , tmp_vector(tmp_vector)
1890  , solution(solution)
1891  {}
1892 
1893  void
1894  apply_to_subrange(const std::size_t begin, const std::size_t end) const
1895  {
1896  // To circumvent a bug in gcc
1897  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
1898  // copies of the variables factor1 and factor2 and do not check based on
1899  // factor1.
1900  const Number factor1 = this->factor1;
1901  const Number factor1_plus_1 = 1. + this->factor1;
1902  const Number factor2 = this->factor2;
1903  if (iteration_index == 0)
1904  {
1906  for (std::size_t i = begin; i < end; ++i)
1907  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
1908  }
1909  else if (iteration_index == 1)
1910  {
1911  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
1913  for (std::size_t i = begin; i < end; ++i)
1914  // for efficiency reason, write back to temp_vector that is
1915  // already read (avoid read-for-ownership)
1916  tmp_vector[i] =
1917  factor1_plus_1 * solution[i] +
1918  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1919  }
1920  else
1921  {
1922  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
1923  // + f_2 * P^{-1} * (b-A*x^{n})
1925  for (std::size_t i = begin; i < end; ++i)
1926  solution_old[i] =
1927  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
1928  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1929  }
1930  }
1931 
1932  const Number * rhs;
1933  const Number * matrix_diagonal_inverse;
1934  const unsigned int iteration_index;
1935  const Number factor1;
1936  const Number factor2;
1937  mutable Number * solution_old;
1938  mutable Number * tmp_vector;
1939  mutable Number * solution;
1940  };
1941 
1942  template <typename Number>
1943  struct VectorUpdatesRange : public ::parallel::ParallelForInteger
1944  {
1945  VectorUpdatesRange(const VectorUpdater<Number> &updater,
1946  const std::size_t size)
1947  : updater(updater)
1948  {
1950  VectorUpdatesRange::apply_to_subrange(0, size);
1951  else
1952  apply_parallel(
1953  0,
1954  size,
1956  }
1957 
1958  ~VectorUpdatesRange() override = default;
1959 
1960  virtual void
1961  apply_to_subrange(const std::size_t begin,
1962  const std::size_t end) const override
1963  {
1964  updater.apply_to_subrange(begin, end);
1965  }
1966 
1967  const VectorUpdater<Number> &updater;
1968  };
1969 
1970  // selection for diagonal matrix around deal.II vector
1971  template <typename Number>
1972  inline void
1973  vector_updates(const ::Vector<Number> & rhs,
1974  const DiagonalMatrix<::Vector<Number>> &jacobi,
1975  const unsigned int iteration_index,
1976  const double factor1,
1977  const double factor2,
1978  ::Vector<Number> &solution_old,
1979  ::Vector<Number> &temp_vector1,
1980  ::Vector<Number> &,
1981  ::Vector<Number> &solution)
1982  {
1983  VectorUpdater<Number> upd(rhs.begin(),
1984  jacobi.get_vector().begin(),
1985  iteration_index,
1986  factor1,
1987  factor2,
1988  solution_old.begin(),
1989  temp_vector1.begin(),
1990  solution.begin());
1991  VectorUpdatesRange<Number>(upd, rhs.size());
1992 
1993  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
1994  if (iteration_index == 0)
1995  {
1996  // nothing to do here because we can immediately write into the
1997  // solution vector without remembering any of the other vectors
1998  }
1999  else if (iteration_index == 1)
2000  {
2001  solution.swap(temp_vector1);
2002  solution_old.swap(temp_vector1);
2003  }
2004  else
2005  solution.swap(solution_old);
2006  }
2007 
2008  // selection for diagonal matrix around parallel deal.II vector
2009  template <typename Number>
2010  inline void
2011  vector_updates(
2013  const DiagonalMatrix<
2015  const unsigned int iteration_index,
2016  const double factor1,
2017  const double factor2,
2019  &solution_old,
2021  &temp_vector1,
2024  {
2025  VectorUpdater<Number> upd(rhs.begin(),
2026  jacobi.get_vector().begin(),
2027  iteration_index,
2028  factor1,
2029  factor2,
2030  solution_old.begin(),
2031  temp_vector1.begin(),
2032  solution.begin());
2033  VectorUpdatesRange<Number>(upd, rhs.local_size());
2034 
2035  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2036  if (iteration_index == 0)
2037  {
2038  // nothing to do here because we can immediately write into the
2039  // solution vector without remembering any of the other vectors
2040  }
2041  else if (iteration_index == 1)
2042  {
2043  solution.swap(temp_vector1);
2044  solution_old.swap(temp_vector1);
2045  }
2046  else
2047  solution.swap(solution_old);
2048  }
2049 
2050  template <typename MatrixType, typename PreconditionerType>
2051  inline void
2052  initialize_preconditioner(
2053  const MatrixType & matrix,
2054  std::shared_ptr<PreconditionerType> &preconditioner)
2055  {
2056  (void)matrix;
2057  (void)preconditioner;
2058  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2059  }
2060 
2061  template <typename MatrixType, typename VectorType>
2062  inline void
2063  initialize_preconditioner(
2064  const MatrixType & matrix,
2065  std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
2066  {
2067  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2068  {
2069  if (preconditioner.get() == nullptr)
2070  preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2071 
2072  Assert(
2073  preconditioner->m() == 0,
2074  ExcMessage(
2075  "Preconditioner appears to be initialized but not sized correctly"));
2076 
2077  // This part only works in serial
2078  if (preconditioner->m() != matrix.m())
2079  {
2080  preconditioner->get_vector().reinit(matrix.m());
2081  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2082  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2083  }
2084  }
2085  }
2086 
2087  template <typename VectorType>
2088  void
2089  set_initial_guess(VectorType &vector)
2090  {
2091  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2092  if (vector.locally_owned_elements().is_element(0))
2093  vector(0) = 0.;
2094  }
2095 
2096  template <typename Number>
2097  void
2098  set_initial_guess(::Vector<Number> &vector)
2099  {
2100  // Choose a high-frequency mode consisting of numbers between 0 and 1
2101  // that is cheap to compute (cheaper than random numbers) but avoids
2102  // obviously re-occurring numbers in multi-component systems by choosing
2103  // a period of 11
2104  for (unsigned int i = 0; i < vector.size(); ++i)
2105  vector(i) = i % 11;
2106 
2107  const Number mean_value = vector.mean_value();
2108  vector.add(-mean_value);
2109  }
2110 
2111  template <typename Number>
2112  void
2113  set_initial_guess(
2115  &vector)
2116  {
2117  // Choose a high-frequency mode consisting of numbers between 0 and 1
2118  // that is cheap to compute (cheaper than random numbers) but avoids
2119  // obviously re-occurring numbers in multi-component systems by choosing
2120  // a period of 11.
2121  // Make initial guess robust with respect to number of processors
2122  // by operating on the global index.
2123  types::global_dof_index first_local_range = 0;
2124  if (!vector.locally_owned_elements().is_empty())
2125  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2126  for (unsigned int i = 0; i < vector.local_size(); ++i)
2127  vector.local_element(i) = (i + first_local_range) % 11;
2128 
2129  const Number mean_value = vector.mean_value();
2130  vector.add(-mean_value);
2131  }
2132 
2133 
2134 # ifdef DEAL_II_COMPILER_CUDA_AWARE
2135  template <typename Number>
2136  __global__ void
2137  set_initial_guess_kernel(const types::global_dof_index offset,
2138  const unsigned int local_size,
2139  Number * values)
2140 
2141  {
2142  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
2143  if (index < local_size)
2144  values[index] = (index + offset) % 11;
2145  }
2146 
2147  template <typename Number>
2148  void
2149  set_initial_guess(
2151  &vector)
2152  {
2153  // Choose a high-frequency mode consisting of numbers between 0 and 1
2154  // that is cheap to compute (cheaper than random numbers) but avoids
2155  // obviously re-occurring numbers in multi-component systems by choosing
2156  // a period of 11.
2157  // Make initial guess robust with respect to number of processors
2158  // by operating on the global index.
2159  types::global_dof_index first_local_range = 0;
2160  if (!vector.locally_owned_elements().is_empty())
2161  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2162 
2163  const auto n_local_elements = vector.local_size();
2164  const int n_blocks =
2165  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
2166  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
2167  first_local_range, n_local_elements, vector.get_values());
2168  AssertCudaKernel();
2169 
2170  const Number mean_value = vector.mean_value();
2171  vector.add(-mean_value);
2172  }
2173 # endif // DEAL_II_COMPILER_CUDA_AWARE
2174 
2175  struct EigenvalueTracker
2176  {
2177  public:
2178  void
2179  slot(const std::vector<double> &eigenvalues)
2180  {
2181  values = eigenvalues;
2182  }
2183 
2184  std::vector<double> values;
2185  };
2186  } // namespace PreconditionChebyshevImplementation
2187 } // namespace internal
2188 
2189 
2190 
2191 template <typename MatrixType, class VectorType, typename PreconditionerType>
2193  AdditionalData::AdditionalData(const unsigned int degree,
2194  const double smoothing_range,
2195  const unsigned int eig_cg_n_iterations,
2196  const double eig_cg_residual,
2197  const double max_eigenvalue)
2198  : degree(degree)
2199  , smoothing_range(smoothing_range)
2200  , eig_cg_n_iterations(eig_cg_n_iterations)
2201  , eig_cg_residual(eig_cg_residual)
2202  , max_eigenvalue(max_eigenvalue)
2203 {}
2204 
2205 
2206 
2207 template <typename MatrixType, class VectorType, typename PreconditionerType>
2208 inline typename PreconditionChebyshev<MatrixType,
2209  VectorType,
2210  PreconditionerType>::AdditionalData &
2212  AdditionalData::operator=(const AdditionalData &other_data)
2213 {
2214  degree = other_data.degree;
2215  smoothing_range = other_data.smoothing_range;
2216  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
2217  eig_cg_residual = other_data.eig_cg_residual;
2218  max_eigenvalue = other_data.max_eigenvalue;
2219  preconditioner = other_data.preconditioner;
2220  constraints.copy_from(other_data.constraints);
2221 
2222  return *this;
2223 }
2224 
2225 
2226 
2227 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2230  : theta(1.)
2231  , delta(1.)
2233 {
2234  static_assert(
2236  "PreconditionChebyshev and VectorType must have the same size_type.");
2237 }
2238 
2239 
2240 
2241 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2242 inline void
2244  const MatrixType & matrix,
2245  const AdditionalData &additional_data)
2246 {
2247  matrix_ptr = &matrix;
2248  data = additional_data;
2249  Assert(data.degree > 0,
2250  ExcMessage("The degree of the Chebyshev method must be positive."));
2251  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2252  matrix, data.preconditioner);
2254 }
2255 
2256 
2257 
2258 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2259 inline void
2261 {
2263  theta = delta = 1.0;
2264  matrix_ptr = nullptr;
2265  {
2266  VectorType empty_vector;
2267  solution_old.reinit(empty_vector);
2268  temp_vector1.reinit(empty_vector);
2269  temp_vector2.reinit(empty_vector);
2270  }
2271  data.preconditioner.reset();
2272 }
2273 
2274 
2275 
2276 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2277 inline typename PreconditionChebyshev<MatrixType,
2278  VectorType,
2279  PreconditionerType>::EigenvalueInformation
2281  estimate_eigenvalues(const VectorType &src) const
2282 {
2284  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2285 
2287  EigenvalueInformation info{};
2288 
2289  solution_old.reinit(src);
2290  temp_vector1.reinit(src, true);
2291 
2292  // calculate largest eigenvalue using a hand-tuned CG iteration on the
2293  // matrix weighted by its diagonal. we start with a vector that consists of
2294  // ones only, weighted by the length.
2295  if (data.eig_cg_n_iterations > 0)
2296  {
2298  ExcMessage(
2299  "Need to set at least two iterations to find eigenvalues."));
2300 
2301  // set a very strict tolerance to force at least two iterations
2302  ReductionControl control(
2304  std::sqrt(
2306  1e-10,
2307  false,
2308  false);
2309 
2310  internal::PreconditionChebyshevImplementation::EigenvalueTracker
2311  eigenvalue_tracker;
2312  SolverCG<VectorType> solver(control);
2313  solver.connect_eigenvalues_slot(
2314  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2315  eigenvalue_tracker.slot(eigenvalues);
2316  });
2317 
2318  // set an initial guess which is close to the constant vector but where
2319  // one entry is different to trigger high frequencies
2320  internal::PreconditionChebyshevImplementation::set_initial_guess(
2321  temp_vector1);
2323 
2324  try
2325  {
2326  solver.solve(*matrix_ptr,
2327  solution_old,
2328  temp_vector1,
2329  *data.preconditioner);
2330  }
2332  {}
2333 
2334  // read the eigenvalues from the attached eigenvalue tracker
2335  if (eigenvalue_tracker.values.empty())
2336  info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2337  else
2338  {
2339  info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2340 
2341  // include a safety factor since the CG method will in general not
2342  // be converged
2343  info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
2344  }
2345 
2346  info.cg_iterations = control.last_step();
2347  }
2348  else
2349  {
2350  info.max_eigenvalue_estimate = data.max_eigenvalue;
2351  info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
2352  }
2353 
2354  const double alpha = (data.smoothing_range > 1. ?
2355  info.max_eigenvalue_estimate / data.smoothing_range :
2356  std::min(0.9 * info.max_eigenvalue_estimate,
2357  info.min_eigenvalue_estimate));
2358 
2359  // in case the user set the degree to invalid unsigned int, we have to
2360  // determine the number of necessary iterations from the Chebyshev error
2361  // estimate, given the target tolerance specified by smoothing_range. This
2362  // estimate is based on the error formula given in section 5.1 of
2363  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
2365  {
2366  const double actual_range = info.max_eigenvalue_estimate / alpha;
2367  const double sigma = (1. - std::sqrt(1. / actual_range)) /
2368  (1. + std::sqrt(1. / actual_range));
2369  const double eps = data.smoothing_range;
2370  const_cast<
2372  this)
2373  ->data.degree =
2374  1 + static_cast<unsigned int>(
2375  std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
2376  std::log(1. / sigma));
2377 
2378  info.degree = data.degree;
2379  }
2380 
2381  const_cast<
2383  ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
2384  const_cast<
2386  ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
2387 
2388  // We do not need the second temporary vector in case we have a
2389  // DiagonalMatrix as preconditioner and use deal.II's own vectors
2390  using NumberType = typename VectorType::value_type;
2391  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
2392  false ||
2393  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
2394  ((std::is_same<VectorType,
2396  Vector<NumberType, MemorySpace::Host>>::value ==
2397  false) ||
2398  (std::is_same<VectorType,
2399  LinearAlgebra::distributed::
2400  Vector<NumberType, MemorySpace::CUDA>>::value ==
2401  false))))
2402  temp_vector2.reinit(src, true);
2403  else
2404  {
2405  VectorType empty_vector;
2406  temp_vector2.reinit(empty_vector);
2407  }
2408 
2409  const_cast<
2411  ->eigenvalues_are_initialized = true;
2412 
2413  return info;
2414 }
2415 
2416 
2417 
2418 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2419 inline void
2421  VectorType & solution,
2422  const VectorType &rhs) const
2423 {
2424  std::lock_guard<std::mutex> lock(mutex);
2425  if (eigenvalues_are_initialized == false)
2426  estimate_eigenvalues(rhs);
2427 
2428  internal::PreconditionChebyshevImplementation::vector_updates(
2429  rhs,
2431  0,
2432  0.,
2433  1. / theta,
2434  solution_old,
2435  temp_vector1,
2436  temp_vector2,
2437  solution);
2438 
2439  // if delta is zero, we do not need to iterate because the updates will be
2440  // zero
2441  if (data.degree < 2 || std::abs(delta) < 1e-40)
2442  return;
2443 
2444  double rhok = delta / theta, sigma = theta / delta;
2445  for (unsigned int k = 0; k < data.degree - 1; ++k)
2446  {
2447  matrix_ptr->vmult(temp_vector1, solution);
2448  const double rhokp = 1. / (2. * sigma - rhok);
2449  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2450  rhok = rhokp;
2451  internal::PreconditionChebyshevImplementation::vector_updates(
2452  rhs,
2454  k + 1,
2455  factor1,
2456  factor2,
2457  solution_old,
2458  temp_vector1,
2459  temp_vector2,
2460  solution);
2461  }
2462 }
2463 
2464 
2465 
2466 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2467 inline void
2469  VectorType & solution,
2470  const VectorType &rhs) const
2471 {
2472  std::lock_guard<std::mutex> lock(mutex);
2473  if (eigenvalues_are_initialized == false)
2474  estimate_eigenvalues(rhs);
2475 
2476  internal::PreconditionChebyshevImplementation::vector_updates(
2477  rhs,
2479  0,
2480  0.,
2481  1. / theta,
2482  solution_old,
2483  temp_vector1,
2484  temp_vector2,
2485  solution);
2486 
2487  if (data.degree < 2 || std::abs(delta) < 1e-40)
2488  return;
2489 
2490  double rhok = delta / theta, sigma = theta / delta;
2491  for (unsigned int k = 0; k < data.degree - 1; ++k)
2492  {
2493  matrix_ptr->Tvmult(temp_vector1, solution);
2494  const double rhokp = 1. / (2. * sigma - rhok);
2495  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2496  rhok = rhokp;
2497  internal::PreconditionChebyshevImplementation::vector_updates(
2498  rhs,
2500  k + 1,
2501  factor1,
2502  factor2,
2503  solution_old,
2504  temp_vector1,
2505  temp_vector2,
2506  solution);
2507  }
2508 }
2509 
2510 
2511 
2512 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2513 inline void
2515  VectorType & solution,
2516  const VectorType &rhs) const
2517 {
2518  std::lock_guard<std::mutex> lock(mutex);
2519  if (eigenvalues_are_initialized == false)
2520  estimate_eigenvalues(rhs);
2521 
2522  matrix_ptr->vmult(temp_vector1, solution);
2523  internal::PreconditionChebyshevImplementation::vector_updates(
2524  rhs,
2526  1,
2527  0.,
2528  1. / theta,
2529  solution_old,
2530  temp_vector1,
2531  temp_vector2,
2532  solution);
2533 
2534  if (data.degree < 2 || std::abs(delta) < 1e-40)
2535  return;
2536 
2537  double rhok = delta / theta, sigma = theta / delta;
2538  for (unsigned int k = 0; k < data.degree - 1; ++k)
2539  {
2540  matrix_ptr->vmult(temp_vector1, solution);
2541  const double rhokp = 1. / (2. * sigma - rhok);
2542  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2543  rhok = rhokp;
2544  internal::PreconditionChebyshevImplementation::vector_updates(
2545  rhs,
2547  k + 2,
2548  factor1,
2549  factor2,
2550  solution_old,
2551  temp_vector1,
2552  temp_vector2,
2553  solution);
2554  }
2555 }
2556 
2557 
2558 
2559 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2560 inline void
2562  VectorType & solution,
2563  const VectorType &rhs) const
2564 {
2565  std::lock_guard<std::mutex> lock(mutex);
2566  if (eigenvalues_are_initialized == false)
2567  estimate_eigenvalues(rhs);
2568 
2569  matrix_ptr->Tvmult(temp_vector1, solution);
2570  internal::PreconditionChebyshevImplementation::vector_updates(
2571  rhs,
2573  1,
2574  0.,
2575  1. / theta,
2576  solution_old,
2577  temp_vector1,
2578  temp_vector2,
2579  solution);
2580 
2581  if (data.degree < 2 || std::abs(delta) < 1e-40)
2582  return;
2583 
2584  double rhok = delta / theta, sigma = theta / delta;
2585  for (unsigned int k = 0; k < data.degree - 1; ++k)
2586  {
2587  matrix_ptr->Tvmult(temp_vector1, solution);
2588  const double rhokp = 1. / (2. * sigma - rhok);
2589  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2590  rhok = rhokp;
2591  internal::PreconditionChebyshevImplementation::vector_updates(
2592  rhs,
2594  k + 2,
2595  factor1,
2596  factor2,
2597  solution_old,
2598  temp_vector1,
2599  temp_vector2,
2600  solution);
2601  }
2602 }
2603 
2604 
2605 
2606 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2607 inline typename PreconditionChebyshev<MatrixType,
2608  VectorType,
2609  PreconditionerType>::size_type
2611 {
2612  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2613  return matrix_ptr->m();
2614 }
2615 
2616 
2617 
2618 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2619 inline typename PreconditionChebyshev<MatrixType,
2620  VectorType,
2621  PreconditionerType>::size_type
2623 {
2624  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2625  return matrix_ptr->n();
2626 }
2627 
2628 #endif // DOXYGEN
2629 
2631 
2632 #endif
void Tvmult(VectorType &, const VectorType &) const
static const unsigned int invalid_unsigned_int
Definition: types.h:191
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:809
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:392
AdditionalData data
AdditionalData(const double relaxation=1.)
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition: precondition.h:372
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:1531
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:1791
size_type m() const
const std::vector< size_type > * permutation
Definition: precondition.h:862
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
const function_ptr precondition
Definition: precondition.h:397
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:417
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:805
#define Assert(cond, exc)
Definition: exceptions.h:1419
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:672
void step(VectorType &x, const VectorType &rhs) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
VectorType::value_type * end(VectorType &V)
typename MatrixType::size_type size_type
Definition: precondition.h:773
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:469
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:361
VectorType::value_type * begin(VectorType &V)
T min(const T &t, const MPI_Comm &mpi_communicator)
typename MatrixType::size_type size_type
Definition: precondition.h:677
const std::vector< size_type > * inverse_permutation
Definition: precondition.h:866
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:515
void vmult(VectorType &, const VectorType &) const
const std::vector< size_type > & permutation
Definition: precondition.h:801
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:603
Eigenvalue vector is filled.
virtual Number mean_value() const override
void vmult_add(VectorType &, const VectorType &) const
static const bool value
std::vector< std::size_t > pos_right_of_diagonal
Definition: precondition.h:730
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1881
bool is_empty() const
Definition: index_set.h:1825
T max(const T &t, const MPI_Comm &mpi_communicator)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:133
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())