Reference documentation for deal.II version Git 24fafe1087 2021-11-29 17:11:51 +0100
\(\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:
504 
508  using AdditionalData =
510 
514  template <class VectorType>
515  void
516  vmult(VectorType &, const VectorType &) const;
517 
522  template <class VectorType>
523  void
524  Tvmult(VectorType &, const VectorType &) const;
525 
529  template <class VectorType>
530  void
531  step(VectorType &x, const VectorType &rhs) const;
532 
536  template <class VectorType>
537  void
538  Tstep(VectorType &x, const VectorType &rhs) const;
539 };
540 
541 
587 template <typename MatrixType = SparseMatrix<double>>
588 class PreconditionSOR : public PreconditionRelaxation<MatrixType>
589 {
590 public:
595 
599  using AdditionalData =
601 
605  template <class VectorType>
606  void
607  vmult(VectorType &, const VectorType &) const;
608 
612  template <class VectorType>
613  void
614  Tvmult(VectorType &, const VectorType &) const;
615 
619  template <class VectorType>
620  void
621  step(VectorType &x, const VectorType &rhs) const;
622 
626  template <class VectorType>
627  void
628  Tstep(VectorType &x, const VectorType &rhs) const;
629 };
630 
631 
632 
659 template <typename MatrixType = SparseMatrix<double>>
660 class PreconditionSSOR : public PreconditionRelaxation<MatrixType>
661 {
662 public:
667 
671  using AdditionalData =
673 
678 
679 
685  void
686  initialize(const MatrixType & A,
687  const typename BaseClass::AdditionalData &parameters =
688  typename BaseClass::AdditionalData());
689 
693  template <class VectorType>
694  void
695  vmult(VectorType &, const VectorType &) const;
696 
701  template <class VectorType>
702  void
703  Tvmult(VectorType &, const VectorType &) const;
704 
705 
709  template <class VectorType>
710  void
711  step(VectorType &x, const VectorType &rhs) const;
712 
716  template <class VectorType>
717  void
718  Tstep(VectorType &x, const VectorType &rhs) const;
719 
720 private:
725  std::vector<std::size_t> pos_right_of_diagonal;
726 };
727 
728 
758 template <typename MatrixType = SparseMatrix<double>>
759 class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
760 {
761 public:
766 
771  {
772  public:
784  const std::vector<size_type> &permutation,
785  const std::vector<size_type> &inverse_permutation,
787  &parameters =
789 
793  const std::vector<size_type> &permutation;
797  const std::vector<size_type> &inverse_permutation;
802  };
803 
815  void
816  initialize(const MatrixType & A,
817  const std::vector<size_type> &permutation,
818  const std::vector<size_type> &inverse_permutation,
820  &parameters =
822 
833  void
834  initialize(const MatrixType &A, const AdditionalData &additional_data);
835 
839  template <class VectorType>
840  void
841  vmult(VectorType &, const VectorType &) const;
842 
846  template <class VectorType>
847  void
848  Tvmult(VectorType &, const VectorType &) const;
849 
850 private:
854  const std::vector<size_type> *permutation;
858  const std::vector<size_type> *inverse_permutation;
859 };
860 
861 
862 
1055 template <typename MatrixType = SparseMatrix<double>,
1056  typename VectorType = Vector<double>,
1057  typename PreconditionerType = DiagonalMatrix<VectorType>>
1059 {
1060 public:
1065 
1071  {
1075  AdditionalData(const unsigned int degree = 1,
1076  const double smoothing_range = 0.,
1077  const unsigned int eig_cg_n_iterations = 8,
1078  const double eig_cg_residual = 1e-2,
1079  const double max_eigenvalue = 1);
1080 
1084  AdditionalData &
1085  operator=(const AdditionalData &other_data);
1086 
1098  unsigned int degree;
1099 
1112 
1119  unsigned int eig_cg_n_iterations;
1120 
1126 
1133 
1139 
1143  std::shared_ptr<PreconditionerType> preconditioner;
1144  };
1145 
1146 
1151 
1163  void
1164  initialize(const MatrixType & matrix,
1165  const AdditionalData &additional_data = AdditionalData());
1166 
1171  void
1172  vmult(VectorType &dst, const VectorType &src) const;
1173 
1178  void
1179  Tvmult(VectorType &dst, const VectorType &src) const;
1180 
1184  void
1185  step(VectorType &dst, const VectorType &src) const;
1186 
1190  void
1191  Tstep(VectorType &dst, const VectorType &src) const;
1192 
1196  void
1197  clear();
1198 
1203  size_type
1204  m() const;
1205 
1210  size_type
1211  n() const;
1212 
1218  {
1230  unsigned int cg_iterations;
1235  unsigned int degree;
1240  : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
1241  , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
1242  , cg_iterations{0}
1243  , degree{0}
1244  {}
1245  };
1246 
1260  estimate_eigenvalues(const VectorType &src) const;
1261 
1262 private:
1266  SmartPointer<
1267  const MatrixType,
1270 
1275 
1280 
1285 
1291 
1295  double theta;
1296 
1301  double delta;
1302 
1308 
1314 };
1315 
1316 
1317 
1319 /* ---------------------------------- Inline functions ------------------- */
1320 
1321 #ifndef DOXYGEN
1322 
1324  : n_rows(0)
1325  , n_columns(0)
1326 {}
1327 
1328 template <typename MatrixType>
1329 inline void
1330 PreconditionIdentity::initialize(const MatrixType &matrix,
1332 {
1333  n_rows = matrix.m();
1334  n_columns = matrix.n();
1335 }
1336 
1337 
1338 template <class VectorType>
1339 inline void
1340 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1341 {
1342  dst = src;
1343 }
1344 
1345 
1346 
1347 template <class VectorType>
1348 inline void
1349 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1350 {
1351  dst = src;
1352 }
1353 
1354 template <class VectorType>
1355 inline void
1357 {
1358  dst += src;
1359 }
1360 
1361 
1362 
1363 template <class VectorType>
1364 inline void
1366 {
1367  dst += src;
1368 }
1369 
1372 {
1373  Assert(n_rows != 0, ExcNotInitialized());
1374  return n_rows;
1375 }
1376 
1379 {
1380  Assert(n_columns != 0, ExcNotInitialized());
1381  return n_columns;
1382 }
1383 
1384 //---------------------------------------------------------------------------
1385 
1387  const double relaxation)
1388  : relaxation(relaxation)
1389 {}
1390 
1391 
1393  : relaxation(0)
1394  , n_rows(0)
1395  , n_columns(0)
1396 {
1397  AdditionalData add_data;
1398  relaxation = add_data.relaxation;
1399 }
1400 
1401 
1402 
1403 inline void
1405  const PreconditionRichardson::AdditionalData &parameters)
1406 {
1407  relaxation = parameters.relaxation;
1408 }
1409 
1410 
1411 
1412 template <typename MatrixType>
1413 inline void
1415  const MatrixType & matrix,
1416  const PreconditionRichardson::AdditionalData &parameters)
1417 {
1418  relaxation = parameters.relaxation;
1419  n_rows = matrix.m();
1420  n_columns = matrix.n();
1421 }
1422 
1423 
1424 
1425 template <class VectorType>
1426 inline void
1427 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
1428 {
1429  static_assert(
1430  std::is_same<size_type, typename VectorType::size_type>::value,
1431  "PreconditionRichardson and VectorType must have the same size_type.");
1432 
1433  dst.equ(relaxation, src);
1434 }
1435 
1436 
1437 
1438 template <class VectorType>
1439 inline void
1441 {
1442  static_assert(
1443  std::is_same<size_type, typename VectorType::size_type>::value,
1444  "PreconditionRichardson and VectorType must have the same size_type.");
1445 
1446  dst.equ(relaxation, src);
1447 }
1448 
1449 template <class VectorType>
1450 inline void
1452 {
1453  static_assert(
1454  std::is_same<size_type, typename VectorType::size_type>::value,
1455  "PreconditionRichardson and VectorType must have the same size_type.");
1456 
1457  dst.add(relaxation, src);
1458 }
1459 
1460 
1461 
1462 template <class VectorType>
1463 inline void
1465 {
1466  static_assert(
1467  std::is_same<size_type, typename VectorType::size_type>::value,
1468  "PreconditionRichardson and VectorType must have the same size_type.");
1469 
1470  dst.add(relaxation, src);
1471 }
1472 
1475 {
1476  Assert(n_rows != 0, ExcNotInitialized());
1477  return n_rows;
1478 }
1479 
1482 {
1483  Assert(n_columns != 0, ExcNotInitialized());
1484  return n_columns;
1485 }
1486 
1487 //---------------------------------------------------------------------------
1488 
1489 template <typename MatrixType>
1490 inline void
1492  const AdditionalData &parameters)
1493 {
1494  A = &rA;
1495  relaxation = parameters.relaxation;
1496 }
1497 
1498 
1499 template <typename MatrixType>
1500 inline void
1502 {
1503  A = nullptr;
1504 }
1505 
1506 template <typename MatrixType>
1509 {
1510  Assert(A != nullptr, ExcNotInitialized());
1511  return A->m();
1512 }
1513 
1514 template <typename MatrixType>
1517 {
1518  Assert(A != nullptr, ExcNotInitialized());
1519  return A->n();
1520 }
1521 
1522 //---------------------------------------------------------------------------
1523 
1524 template <typename MatrixType>
1525 template <class VectorType>
1526 inline void
1528  const VectorType &src) const
1529 {
1530  static_assert(
1531  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1532  typename VectorType::size_type>::value,
1533  "PreconditionJacobi and VectorType must have the same size_type.");
1534 
1535  Assert(this->A != nullptr, ExcNotInitialized());
1536  this->A->precondition_Jacobi(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(
1548  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1549  typename VectorType::size_type>::value,
1550  "PreconditionJacobi and VectorType must have the same size_type.");
1551 
1552  Assert(this->A != nullptr, ExcNotInitialized());
1553  this->A->precondition_Jacobi(dst, src, this->relaxation);
1554 }
1555 
1556 
1557 
1558 template <typename MatrixType>
1559 template <class VectorType>
1560 inline void
1562  const VectorType &src) const
1563 {
1564  static_assert(
1565  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1566  typename VectorType::size_type>::value,
1567  "PreconditionJacobi and VectorType must have the same size_type.");
1568 
1569  Assert(this->A != nullptr, ExcNotInitialized());
1570  this->A->Jacobi_step(dst, src, this->relaxation);
1571 }
1572 
1573 
1574 
1575 template <typename MatrixType>
1576 template <class VectorType>
1577 inline void
1579  const VectorType &src) const
1580 {
1581  static_assert(
1582  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1583  typename VectorType::size_type>::value,
1584  "PreconditionJacobi and VectorType must have the same size_type.");
1585 
1586  step(dst, src);
1587 }
1588 
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->precondition_SOR(dst, src, this->relaxation);
1604 }
1605 
1606 
1607 
1608 template <typename MatrixType>
1609 template <class VectorType>
1610 inline void
1612  const VectorType &src) const
1613 {
1614  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1615  typename VectorType::size_type>::value,
1616  "PreconditionSOR and VectorType must have the same size_type.");
1617 
1618  Assert(this->A != nullptr, ExcNotInitialized());
1619  this->A->precondition_TSOR(dst, src, this->relaxation);
1620 }
1621 
1622 
1623 
1624 template <typename MatrixType>
1625 template <class VectorType>
1626 inline void
1628 {
1629  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1630  typename VectorType::size_type>::value,
1631  "PreconditionSOR and VectorType must have the same size_type.");
1632 
1633  Assert(this->A != nullptr, ExcNotInitialized());
1634  this->A->SOR_step(dst, src, this->relaxation);
1635 }
1636 
1637 
1638 
1639 template <typename MatrixType>
1640 template <class VectorType>
1641 inline void
1643 {
1644  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1645  typename VectorType::size_type>::value,
1646  "PreconditionSOR and VectorType must have the same size_type.");
1647 
1648  Assert(this->A != nullptr, ExcNotInitialized());
1649  this->A->TSOR_step(dst, src, this->relaxation);
1650 }
1651 
1652 
1653 
1654 //---------------------------------------------------------------------------
1655 
1656 template <typename MatrixType>
1657 inline void
1659  const MatrixType & rA,
1660  const typename BaseClass::AdditionalData &parameters)
1661 {
1662  this->PreconditionRelaxation<MatrixType>::initialize(rA, parameters);
1663 
1664  // in case we have a SparseMatrix class, we can extract information about
1665  // the diagonal.
1667  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
1668  &*this->A);
1669 
1670  // calculate the positions first after the diagonal.
1671  if (mat != nullptr)
1672  {
1673  const size_type n = this->A->n();
1674  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1675  for (size_type row = 0; row < n; ++row)
1676  {
1677  // find the first element in this line which is on the right of the
1678  // diagonal. we need to precondition with the elements on the left
1679  // only. note: the first entry in each line denotes the diagonal
1680  // element, which we need not check.
1682  it = mat->begin(row) + 1;
1683  for (; it < mat->end(row); ++it)
1684  if (it->column() > row)
1685  break;
1686  pos_right_of_diagonal[row] = it - mat->begin();
1687  }
1688  }
1689 }
1690 
1691 
1692 template <typename MatrixType>
1693 template <class VectorType>
1694 inline void
1696  const VectorType &src) const
1697 {
1698  static_assert(
1699  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1700  typename VectorType::size_type>::value,
1701  "PreconditionSSOR and VectorType must have the same size_type.");
1702 
1703  Assert(this->A != nullptr, ExcNotInitialized());
1704  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1705 }
1706 
1707 
1708 
1709 template <typename MatrixType>
1710 template <class VectorType>
1711 inline void
1713  const VectorType &src) const
1714 {
1715  static_assert(
1716  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1717  typename VectorType::size_type>::value,
1718  "PreconditionSSOR and VectorType must have the same size_type.");
1719 
1720  Assert(this->A != nullptr, ExcNotInitialized());
1721  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1722 }
1723 
1724 
1725 
1726 template <typename MatrixType>
1727 template <class VectorType>
1728 inline void
1730 {
1731  static_assert(
1732  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1733  typename VectorType::size_type>::value,
1734  "PreconditionSSOR and VectorType must have the same size_type.");
1735 
1736  Assert(this->A != nullptr, ExcNotInitialized());
1737  this->A->SSOR_step(dst, src, this->relaxation);
1738 }
1739 
1740 
1741 
1742 template <typename MatrixType>
1743 template <class VectorType>
1744 inline void
1746  const VectorType &src) const
1747 {
1748  static_assert(
1749  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1750  typename VectorType::size_type>::value,
1751  "PreconditionSSOR and VectorType must have the same size_type.");
1752 
1753  step(dst, src);
1754 }
1755 
1756 
1757 
1758 //---------------------------------------------------------------------------
1759 
1760 template <typename MatrixType>
1761 inline void
1763  const MatrixType & rA,
1764  const std::vector<size_type> & p,
1765  const std::vector<size_type> & ip,
1766  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1767 {
1768  permutation = &p;
1769  inverse_permutation = &ip;
1771 }
1772 
1773 
1774 template <typename MatrixType>
1775 inline void
1776 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
1777  const AdditionalData &additional_data)
1778 {
1779  initialize(A,
1780  additional_data.permutation,
1781  additional_data.inverse_permutation,
1782  additional_data.parameters);
1783 }
1784 
1785 
1786 template <typename MatrixType>
1787 template <typename VectorType>
1788 inline void
1790  const VectorType &src) const
1791 {
1792  static_assert(
1793  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1794  typename VectorType::size_type>::value,
1795  "PreconditionPSOR and VectorType must have the same size_type.");
1796 
1797  Assert(this->A != nullptr, ExcNotInitialized());
1798  dst = src;
1799  this->A->PSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1800 }
1801 
1802 
1803 
1804 template <typename MatrixType>
1805 template <class VectorType>
1806 inline void
1808  const VectorType &src) const
1809 {
1810  static_assert(
1811  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1812  typename VectorType::size_type>::value,
1813  "PreconditionPSOR and VectorType must have the same size_type.");
1814 
1815  Assert(this->A != nullptr, ExcNotInitialized());
1816  dst = src;
1817  this->A->TPSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1818 }
1819 
1820 template <typename MatrixType>
1822  const std::vector<size_type> &permutation,
1823  const std::vector<size_type> &inverse_permutation,
1824  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1825  : permutation(permutation)
1826  , inverse_permutation(inverse_permutation)
1827  , parameters(parameters)
1828 {}
1829 
1830 
1831 //---------------------------------------------------------------------------
1832 
1833 
1834 template <typename MatrixType, class VectorType>
1836  const MatrixType & M,
1837  const function_ptr method)
1838  : matrix(M)
1839  , precondition(method)
1840 {}
1841 
1842 
1843 
1844 template <typename MatrixType, class VectorType>
1845 void
1847  VectorType & dst,
1848  const VectorType &src) const
1849 {
1850  (matrix.*precondition)(dst, src);
1851 }
1852 
1853 //---------------------------------------------------------------------------
1854 
1855 template <typename MatrixType>
1857  const double relaxation)
1858  : relaxation(relaxation)
1859 {}
1860 
1861 
1862 
1863 //---------------------------------------------------------------------------
1864 
1865 namespace internal
1866 {
1867  namespace PreconditionChebyshevImplementation
1868  {
1869  // for deal.II vectors, perform updates for Chebyshev preconditioner all
1870  // at once to reduce memory transfer. Here, we select between general
1871  // vectors and deal.II vectors where we expand the loop over the (local)
1872  // size of the vector
1873 
1874  // generic part for non-deal.II vectors
1875  template <typename VectorType, typename PreconditionerType>
1876  inline void
1877  vector_updates(const VectorType & rhs,
1878  const PreconditionerType &preconditioner,
1879  const unsigned int iteration_index,
1880  const double factor1,
1881  const double factor2,
1885  VectorType & solution)
1886  {
1887  if (iteration_index == 0)
1888  {
1889  solution.equ(factor2, rhs);
1890  preconditioner.vmult(solution_old, solution);
1891  }
1892  else if (iteration_index == 1)
1893  {
1894  // compute t = P^{-1} * (b-A*x^{n})
1895  temp_vector1.sadd(-1.0, 1.0, rhs);
1896  preconditioner.vmult(solution_old, temp_vector1);
1897 
1898  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
1899  solution_old.sadd(factor2, 1 + factor1, solution);
1900  }
1901  else
1902  {
1903  // compute t = P^{-1} * (b-A*x^{n})
1904  temp_vector1.sadd(-1.0, 1.0, rhs);
1905  preconditioner.vmult(temp_vector2, temp_vector1);
1906 
1907  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
1908  solution_old.sadd(-factor1, factor2, temp_vector2);
1909  solution_old.add(1 + factor1, solution);
1910  }
1911 
1912  solution.swap(solution_old);
1913  }
1914 
1915  // worker routine for deal.II vectors. Because of vectorization, we need
1916  // to put the loop into an extra structure because the virtual function of
1917  // VectorUpdatesRange prevents the compiler from applying vectorization.
1918  template <typename Number>
1919  struct VectorUpdater
1920  {
1921  VectorUpdater(const Number * rhs,
1922  const Number * matrix_diagonal_inverse,
1923  const unsigned int iteration_index,
1924  const Number factor1,
1925  const Number factor2,
1926  Number * solution_old,
1927  Number * tmp_vector,
1928  Number * solution)
1929  : rhs(rhs)
1930  , matrix_diagonal_inverse(matrix_diagonal_inverse)
1931  , iteration_index(iteration_index)
1932  , factor1(factor1)
1933  , factor2(factor2)
1934  , solution_old(solution_old)
1935  , tmp_vector(tmp_vector)
1936  , solution(solution)
1937  {}
1938 
1939  void
1940  apply_to_subrange(const std::size_t begin, const std::size_t end) const
1941  {
1942  // To circumvent a bug in gcc
1943  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
1944  // copies of the variables factor1 and factor2 and do not check based on
1945  // factor1.
1946  const Number factor1 = this->factor1;
1947  const Number factor1_plus_1 = 1. + this->factor1;
1948  const Number factor2 = this->factor2;
1949  if (iteration_index == 0)
1950  {
1952  for (std::size_t i = begin; i < end; ++i)
1953  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
1954  }
1955  else if (iteration_index == 1)
1956  {
1957  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
1959  for (std::size_t i = begin; i < end; ++i)
1960  // for efficiency reason, write back to temp_vector that is
1961  // already read (avoid read-for-ownership)
1962  tmp_vector[i] =
1963  factor1_plus_1 * solution[i] +
1964  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1965  }
1966  else
1967  {
1968  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
1969  // + f_2 * P^{-1} * (b-A*x^{n})
1971  for (std::size_t i = begin; i < end; ++i)
1972  solution_old[i] =
1973  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
1974  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1975  }
1976  }
1977 
1978  const Number * rhs;
1979  const Number * matrix_diagonal_inverse;
1980  const unsigned int iteration_index;
1981  const Number factor1;
1982  const Number factor2;
1983  mutable Number * solution_old;
1984  mutable Number * tmp_vector;
1985  mutable Number * solution;
1986  };
1987 
1988  template <typename Number>
1989  struct VectorUpdatesRange : public ::parallel::ParallelForInteger
1990  {
1991  VectorUpdatesRange(const VectorUpdater<Number> &updater,
1992  const std::size_t size)
1993  : updater(updater)
1994  {
1996  VectorUpdatesRange::apply_to_subrange(0, size);
1997  else
1998  apply_parallel(
1999  0,
2000  size,
2002  }
2003 
2004  ~VectorUpdatesRange() override = default;
2005 
2006  virtual void
2007  apply_to_subrange(const std::size_t begin,
2008  const std::size_t end) const override
2009  {
2010  updater.apply_to_subrange(begin, end);
2011  }
2012 
2013  const VectorUpdater<Number> &updater;
2014  };
2015 
2016  // selection for diagonal matrix around deal.II vector
2017  template <typename Number>
2018  inline void
2019  vector_updates(const ::Vector<Number> & rhs,
2020  const DiagonalMatrix<::Vector<Number>> &jacobi,
2021  const unsigned int iteration_index,
2022  const double factor1,
2023  const double factor2,
2024  ::Vector<Number> &solution_old,
2025  ::Vector<Number> &temp_vector1,
2026  ::Vector<Number> &,
2027  ::Vector<Number> &solution)
2028  {
2029  VectorUpdater<Number> upd(rhs.begin(),
2030  jacobi.get_vector().begin(),
2031  iteration_index,
2032  factor1,
2033  factor2,
2034  solution_old.begin(),
2035  temp_vector1.begin(),
2036  solution.begin());
2037  VectorUpdatesRange<Number>(upd, rhs.size());
2038 
2039  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2040  if (iteration_index == 0)
2041  {
2042  // nothing to do here because we can immediately write into the
2043  // solution vector without remembering any of the other vectors
2044  }
2045  else if (iteration_index == 1)
2046  {
2047  solution.swap(temp_vector1);
2048  solution_old.swap(temp_vector1);
2049  }
2050  else
2051  solution.swap(solution_old);
2052  }
2053 
2054  // selection for diagonal matrix around parallel deal.II vector
2055  template <typename Number>
2056  inline void
2057  vector_updates(
2059  const DiagonalMatrix<
2061  const unsigned int iteration_index,
2062  const double factor1,
2063  const double factor2,
2065  &solution_old,
2067  &temp_vector1,
2070  {
2071  VectorUpdater<Number> upd(rhs.begin(),
2072  jacobi.get_vector().begin(),
2073  iteration_index,
2074  factor1,
2075  factor2,
2076  solution_old.begin(),
2077  temp_vector1.begin(),
2078  solution.begin());
2079  VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
2080 
2081  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2082  if (iteration_index == 0)
2083  {
2084  // nothing to do here because we can immediately write into the
2085  // solution vector without remembering any of the other vectors
2086  }
2087  else if (iteration_index == 1)
2088  {
2089  solution.swap(temp_vector1);
2090  solution_old.swap(temp_vector1);
2091  }
2092  else
2093  solution.swap(solution_old);
2094  }
2095 
2096  // Detector class to find out whether we have a vmult function that takes
2097  // two additional std::function objects, which we use to run the operation
2098  // on slices of the vector during the matrix-vector product
2099  template <typename MatrixType,
2100  typename VectorType,
2101  typename PreconditionerType>
2102  struct supports_vmult_with_std_functions
2103  {
2104  private:
2105  // this will work always
2106  static bool
2107  detect(...);
2108 
2109  // this detector will work only if we have
2110  // "... MatrixType::vmult(VectorType, const VectorType,
2111  // const std::function<...>&, const std::function<...>&) const"
2112  template <typename MatrixType2>
2113  static decltype(std::declval<MatrixType2 const>().vmult(
2114  std::declval<VectorType &>(),
2115  std::declval<const VectorType &>(),
2116  std::declval<const std::function<void(const unsigned int,
2117  const unsigned int)> &>(),
2118  std::declval<const std::function<void(const unsigned int,
2119  const unsigned int)> &>()))
2120  detect(const MatrixType2 &);
2121 
2122  public:
2123  // finally here we check if our detector has void return type and
2124  // fulfills additional requirements on the vector type and
2125  // preconditioner. This will happen if the compiler can use the second
2126  // detector, otherwise SFINAE let's it work with the more general first
2127  // one that is bool
2128  static const bool value =
2129  !std::is_same<bool,
2130  decltype(detect(std::declval<MatrixType>()))>::value &&
2131  std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value &&
2132  std::is_same<
2133  VectorType,
2134  LinearAlgebra::distributed::Vector<typename VectorType::value_type,
2135  MemorySpace::Host>>::value;
2136  };
2137 
2138  // We need to have a separate declaration for static const members
2139  template <typename T, typename U, typename V>
2140  const bool supports_vmult_with_std_functions<T, U, V>::value;
2141 
2142  template <typename MatrixType,
2143  typename VectorType,
2144  typename PreconditionerType,
2145  typename std::enable_if<
2146  !supports_vmult_with_std_functions<MatrixType,
2147  VectorType,
2148  PreconditionerType>::value,
2149  int>::type * = nullptr>
2150  inline void
2151  vmult_and_update(const MatrixType & matrix,
2152  const PreconditionerType &preconditioner,
2153  const VectorType & rhs,
2154  const unsigned int iteration_index,
2155  const double factor1,
2156  const double factor2,
2157  VectorType & solution,
2158  VectorType & solution_old,
2159  VectorType & temp_vector1,
2160  VectorType & temp_vector2)
2161  {
2162  if (iteration_index > 0)
2163  matrix.vmult(temp_vector1, solution);
2164  vector_updates(rhs,
2165  preconditioner,
2166  iteration_index,
2167  factor1,
2168  factor2,
2169  solution_old,
2170  temp_vector1,
2171  temp_vector2,
2172  solution);
2173  }
2174 
2175  template <typename MatrixType,
2176  typename VectorType,
2177  typename PreconditionerType,
2178  typename std::enable_if<
2179  supports_vmult_with_std_functions<MatrixType,
2180  VectorType,
2181  PreconditionerType>::value,
2182  int>::type * = nullptr>
2183  inline void
2184  vmult_and_update(const MatrixType & matrix,
2185  const PreconditionerType &preconditioner,
2186  const VectorType & rhs,
2187  const unsigned int iteration_index,
2188  const double factor1,
2189  const double factor2,
2190  VectorType & solution,
2191  VectorType & solution_old,
2192  VectorType & temp_vector1,
2193  VectorType &)
2194  {
2195  using Number = typename VectorType::value_type;
2196  VectorUpdater<Number> updater(rhs.begin(),
2197  preconditioner.get_vector().begin(),
2198  iteration_index,
2199  factor1,
2200  factor2,
2201  solution_old.begin(),
2202  temp_vector1.begin(),
2203  solution.begin());
2204  if (iteration_index > 0)
2205  matrix.vmult(
2206  temp_vector1,
2207  solution,
2208  [&](const unsigned int start_range, const unsigned int end_range) {
2209  // zero 'temp_vector1' before running the vmult
2210  // operation
2211  if (end_range > start_range)
2212  std::memset(temp_vector1.begin() + start_range,
2213  0,
2214  sizeof(Number) * (end_range - start_range));
2215  },
2216  [&](const unsigned int start_range, const unsigned int end_range) {
2217  if (end_range > start_range)
2218  updater.apply_to_subrange(start_range, end_range);
2219  });
2220  else
2221  updater.apply_to_subrange(0U, solution.locally_owned_size());
2222 
2223  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2224  if (iteration_index == 0)
2225  {
2226  // nothing to do here because we can immediately write into the
2227  // solution vector without remembering any of the other vectors
2228  }
2229  else if (iteration_index == 1)
2230  {
2231  solution.swap(temp_vector1);
2232  solution_old.swap(temp_vector1);
2233  }
2234  else
2235  solution.swap(solution_old);
2236  }
2237 
2238  template <typename MatrixType, typename PreconditionerType>
2239  inline void
2240  initialize_preconditioner(
2241  const MatrixType & matrix,
2242  std::shared_ptr<PreconditionerType> &preconditioner)
2243  {
2244  (void)matrix;
2245  (void)preconditioner;
2246  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2247  }
2248 
2249  template <typename MatrixType, typename VectorType>
2250  inline void
2251  initialize_preconditioner(
2252  const MatrixType & matrix,
2253  std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
2254  {
2255  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2256  {
2257  if (preconditioner.get() == nullptr)
2258  preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2259 
2260  Assert(
2261  preconditioner->m() == 0,
2262  ExcMessage(
2263  "Preconditioner appears to be initialized but not sized correctly"));
2264 
2265  // This part only works in serial
2266  if (preconditioner->m() != matrix.m())
2267  {
2268  preconditioner->get_vector().reinit(matrix.m());
2269  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2270  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2271  }
2272  }
2273  }
2274 
2275  template <typename VectorType>
2276  void
2277  set_initial_guess(VectorType &vector)
2278  {
2279  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2280  if (vector.locally_owned_elements().is_element(0))
2281  vector(0) = 0.;
2282  }
2283 
2284  template <typename Number>
2285  void
2286  set_initial_guess(::Vector<Number> &vector)
2287  {
2288  // Choose a high-frequency mode consisting of numbers between 0 and 1
2289  // that is cheap to compute (cheaper than random numbers) but avoids
2290  // obviously re-occurring numbers in multi-component systems by choosing
2291  // a period of 11
2292  for (unsigned int i = 0; i < vector.size(); ++i)
2293  vector(i) = i % 11;
2294 
2295  const Number mean_value = vector.mean_value();
2296  vector.add(-mean_value);
2297  }
2298 
2299  template <typename Number>
2300  void
2301  set_initial_guess(
2303  &vector)
2304  {
2305  // Choose a high-frequency mode consisting of numbers between 0 and 1
2306  // that is cheap to compute (cheaper than random numbers) but avoids
2307  // obviously re-occurring numbers in multi-component systems by choosing
2308  // a period of 11.
2309  // Make initial guess robust with respect to number of processors
2310  // by operating on the global index.
2311  types::global_dof_index first_local_range = 0;
2312  if (!vector.locally_owned_elements().is_empty())
2313  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2314  for (unsigned int i = 0; i < vector.locally_owned_size(); ++i)
2315  vector.local_element(i) = (i + first_local_range) % 11;
2316 
2317  const Number mean_value = vector.mean_value();
2318  vector.add(-mean_value);
2319  }
2320 
2321  template <typename Number>
2322  void
2323  set_initial_guess(
2325  {
2326  for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2327  set_initial_guess(vector.block(block));
2328  }
2329 
2330 
2331 # ifdef DEAL_II_COMPILER_CUDA_AWARE
2332  template <typename Number>
2333  __global__ void
2334  set_initial_guess_kernel(const types::global_dof_index offset,
2335  const unsigned int locally_owned_size,
2336  Number * values)
2337 
2338  {
2339  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
2340  if (index < locally_owned_size)
2341  values[index] = (index + offset) % 11;
2342  }
2343 
2344  template <typename Number>
2345  void
2346  set_initial_guess(
2348  &vector)
2349  {
2350  // Choose a high-frequency mode consisting of numbers between 0 and 1
2351  // that is cheap to compute (cheaper than random numbers) but avoids
2352  // obviously re-occurring numbers in multi-component systems by choosing
2353  // a period of 11.
2354  // Make initial guess robust with respect to number of processors
2355  // by operating on the global index.
2356  types::global_dof_index first_local_range = 0;
2357  if (!vector.locally_owned_elements().is_empty())
2358  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2359 
2360  const auto n_local_elements = vector.locally_owned_size();
2361  const int n_blocks =
2362  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
2363  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
2364  first_local_range, n_local_elements, vector.get_values());
2365  AssertCudaKernel();
2366 
2367  const Number mean_value = vector.mean_value();
2368  vector.add(-mean_value);
2369  }
2370 # endif // DEAL_II_COMPILER_CUDA_AWARE
2371 
2372  struct EigenvalueTracker
2373  {
2374  public:
2375  void
2376  slot(const std::vector<double> &eigenvalues)
2377  {
2378  values = eigenvalues;
2379  }
2380 
2381  std::vector<double> values;
2382  };
2383  } // namespace PreconditionChebyshevImplementation
2384 } // namespace internal
2385 
2386 
2387 
2388 template <typename MatrixType, class VectorType, typename PreconditionerType>
2390  AdditionalData::AdditionalData(const unsigned int degree,
2391  const double smoothing_range,
2392  const unsigned int eig_cg_n_iterations,
2393  const double eig_cg_residual,
2394  const double max_eigenvalue)
2395  : degree(degree)
2396  , smoothing_range(smoothing_range)
2397  , eig_cg_n_iterations(eig_cg_n_iterations)
2398  , eig_cg_residual(eig_cg_residual)
2399  , max_eigenvalue(max_eigenvalue)
2400 {}
2401 
2402 
2403 
2404 template <typename MatrixType, class VectorType, typename PreconditionerType>
2405 inline typename PreconditionChebyshev<MatrixType,
2406  VectorType,
2407  PreconditionerType>::AdditionalData &
2409  AdditionalData::operator=(const AdditionalData &other_data)
2410 {
2411  degree = other_data.degree;
2412  smoothing_range = other_data.smoothing_range;
2413  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
2414  eig_cg_residual = other_data.eig_cg_residual;
2415  max_eigenvalue = other_data.max_eigenvalue;
2416  preconditioner = other_data.preconditioner;
2417  constraints.copy_from(other_data.constraints);
2418 
2419  return *this;
2420 }
2421 
2422 
2423 
2424 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2427  : theta(1.)
2428  , delta(1.)
2430 {
2431  static_assert(
2432  std::is_same<size_type, typename VectorType::size_type>::value,
2433  "PreconditionChebyshev and VectorType must have the same size_type.");
2434 }
2435 
2436 
2437 
2438 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2439 inline void
2441  const MatrixType & matrix,
2442  const AdditionalData &additional_data)
2443 {
2444  matrix_ptr = &matrix;
2445  data = additional_data;
2446  Assert(data.degree > 0,
2447  ExcMessage("The degree of the Chebyshev method must be positive."));
2448  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2449  matrix, data.preconditioner);
2451 }
2452 
2453 
2454 
2455 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2456 inline void
2458 {
2460  theta = delta = 1.0;
2461  matrix_ptr = nullptr;
2462  {
2463  VectorType empty_vector;
2464  solution_old.reinit(empty_vector);
2465  temp_vector1.reinit(empty_vector);
2466  temp_vector2.reinit(empty_vector);
2467  }
2468  data.preconditioner.reset();
2469 }
2470 
2471 
2472 
2473 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2474 inline typename PreconditionChebyshev<MatrixType,
2475  VectorType,
2476  PreconditionerType>::EigenvalueInformation
2478  estimate_eigenvalues(const VectorType &src) const
2479 {
2481  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2482 
2484  EigenvalueInformation info{};
2485 
2486  solution_old.reinit(src);
2487  temp_vector1.reinit(src, true);
2488 
2489  if (data.eig_cg_n_iterations > 0)
2490  {
2492  ExcMessage(
2493  "Need to set at least two iterations to find eigenvalues."));
2494 
2495  // set a very strict tolerance to force at least two iterations
2496  ReductionControl control(
2498  std::sqrt(
2500  1e-10,
2501  false,
2502  false);
2503 
2504  internal::PreconditionChebyshevImplementation::EigenvalueTracker
2505  eigenvalue_tracker;
2506  SolverCG<VectorType> solver(control);
2507  solver.connect_eigenvalues_slot(
2508  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2509  eigenvalue_tracker.slot(eigenvalues);
2510  });
2511 
2512  // set an initial guess that contains some high-frequency parts (to the
2513  // extent possible without knowing the discretization and the numbering)
2514  // to trigger high eigenvalues according to the external function
2515  internal::PreconditionChebyshevImplementation::set_initial_guess(
2516  temp_vector1);
2518 
2519  try
2520  {
2521  solver.solve(*matrix_ptr,
2522  solution_old,
2523  temp_vector1,
2524  *data.preconditioner);
2525  }
2527  {}
2528 
2529  // read the eigenvalues from the attached eigenvalue tracker
2530  if (eigenvalue_tracker.values.empty())
2531  info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2532  else
2533  {
2534  info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2535 
2536  // include a safety factor since the CG method will in general not
2537  // be converged
2538  info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
2539  }
2540 
2541  info.cg_iterations = control.last_step();
2542  }
2543  else
2544  {
2545  info.max_eigenvalue_estimate = data.max_eigenvalue;
2546  info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
2547  }
2548 
2549  const double alpha = (data.smoothing_range > 1. ?
2550  info.max_eigenvalue_estimate / data.smoothing_range :
2551  std::min(0.9 * info.max_eigenvalue_estimate,
2552  info.min_eigenvalue_estimate));
2553 
2554  // in case the user set the degree to invalid unsigned int, we have to
2555  // determine the number of necessary iterations from the Chebyshev error
2556  // estimate, given the target tolerance specified by smoothing_range. This
2557  // estimate is based on the error formula given in section 5.1 of
2558  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
2560  {
2561  const double actual_range = info.max_eigenvalue_estimate / alpha;
2562  const double sigma = (1. - std::sqrt(1. / actual_range)) /
2563  (1. + std::sqrt(1. / actual_range));
2564  const double eps = data.smoothing_range;
2565  const_cast<
2567  this)
2568  ->data.degree =
2569  1 + static_cast<unsigned int>(
2570  std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
2571  std::log(1. / sigma));
2572  }
2573 
2574  info.degree = data.degree;
2575 
2576  const_cast<
2578  ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
2579  const_cast<
2581  ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
2582 
2583  // We do not need the second temporary vector in case we have a
2584  // DiagonalMatrix as preconditioner and use deal.II's own vectors
2585  using NumberType = typename VectorType::value_type;
2586  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
2587  false ||
2588  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
2589  ((std::is_same<VectorType,
2591  Vector<NumberType, MemorySpace::Host>>::value ==
2592  false) ||
2593  (std::is_same<VectorType,
2594  LinearAlgebra::distributed::
2595  Vector<NumberType, MemorySpace::CUDA>>::value ==
2596  false))))
2597  temp_vector2.reinit(src, true);
2598  else
2599  {
2600  VectorType empty_vector;
2601  temp_vector2.reinit(empty_vector);
2602  }
2603 
2604  const_cast<
2606  ->eigenvalues_are_initialized = true;
2607 
2608  return info;
2609 }
2610 
2611 
2612 
2613 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2614 inline void
2616  VectorType & solution,
2617  const VectorType &rhs) const
2618 {
2619  std::lock_guard<std::mutex> lock(mutex);
2620  if (eigenvalues_are_initialized == false)
2621  estimate_eigenvalues(rhs);
2622 
2623  internal::PreconditionChebyshevImplementation::vmult_and_update(
2624  *matrix_ptr,
2626  rhs,
2627  0,
2628  0.,
2629  1. / theta,
2630  solution,
2631  solution_old,
2632  temp_vector1,
2633  temp_vector2);
2634 
2635  // if delta is zero, we do not need to iterate because the updates will be
2636  // zero
2637  if (data.degree < 2 || std::abs(delta) < 1e-40)
2638  return;
2639 
2640  double rhok = delta / theta, sigma = theta / delta;
2641  for (unsigned int k = 0; k < data.degree - 1; ++k)
2642  {
2643  const double rhokp = 1. / (2. * sigma - rhok);
2644  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2645  rhok = rhokp;
2646  internal::PreconditionChebyshevImplementation::vmult_and_update(
2647  *matrix_ptr,
2649  rhs,
2650  k + 1,
2651  factor1,
2652  factor2,
2653  solution,
2654  solution_old,
2655  temp_vector1,
2656  temp_vector2);
2657  }
2658 }
2659 
2660 
2661 
2662 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2663 inline void
2665  VectorType & solution,
2666  const VectorType &rhs) const
2667 {
2668  std::lock_guard<std::mutex> lock(mutex);
2669  if (eigenvalues_are_initialized == false)
2670  estimate_eigenvalues(rhs);
2671 
2672  internal::PreconditionChebyshevImplementation::vector_updates(
2673  rhs,
2675  0,
2676  0.,
2677  1. / theta,
2678  solution_old,
2679  temp_vector1,
2680  temp_vector2,
2681  solution);
2682 
2683  if (data.degree < 2 || std::abs(delta) < 1e-40)
2684  return;
2685 
2686  double rhok = delta / theta, sigma = theta / delta;
2687  for (unsigned int k = 0; k < data.degree - 1; ++k)
2688  {
2689  const double rhokp = 1. / (2. * sigma - rhok);
2690  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2691  rhok = rhokp;
2692  matrix_ptr->Tvmult(temp_vector1, solution);
2693  internal::PreconditionChebyshevImplementation::vector_updates(
2694  rhs,
2696  k + 1,
2697  factor1,
2698  factor2,
2699  solution_old,
2700  temp_vector1,
2701  temp_vector2,
2702  solution);
2703  }
2704 }
2705 
2706 
2707 
2708 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2709 inline void
2711  VectorType & solution,
2712  const VectorType &rhs) const
2713 {
2714  std::lock_guard<std::mutex> lock(mutex);
2715  if (eigenvalues_are_initialized == false)
2716  estimate_eigenvalues(rhs);
2717 
2718  internal::PreconditionChebyshevImplementation::vmult_and_update(
2719  *matrix_ptr,
2721  rhs,
2722  1,
2723  0.,
2724  1. / theta,
2725  solution,
2726  solution_old,
2727  temp_vector1,
2728  temp_vector2);
2729 
2730  if (data.degree < 2 || std::abs(delta) < 1e-40)
2731  return;
2732 
2733  double rhok = delta / theta, sigma = theta / delta;
2734  for (unsigned int k = 0; k < data.degree - 1; ++k)
2735  {
2736  const double rhokp = 1. / (2. * sigma - rhok);
2737  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2738  rhok = rhokp;
2739  internal::PreconditionChebyshevImplementation::vmult_and_update(
2740  *matrix_ptr,
2742  rhs,
2743  k + 2,
2744  factor1,
2745  factor2,
2746  solution,
2747  solution_old,
2748  temp_vector1,
2749  temp_vector2);
2750  }
2751 }
2752 
2753 
2754 
2755 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2756 inline void
2758  VectorType & solution,
2759  const VectorType &rhs) const
2760 {
2761  std::lock_guard<std::mutex> lock(mutex);
2762  if (eigenvalues_are_initialized == false)
2763  estimate_eigenvalues(rhs);
2764 
2765  matrix_ptr->Tvmult(temp_vector1, solution);
2766  internal::PreconditionChebyshevImplementation::vector_updates(
2767  rhs,
2769  1,
2770  0.,
2771  1. / theta,
2772  solution_old,
2773  temp_vector1,
2774  temp_vector2,
2775  solution);
2776 
2777  if (data.degree < 2 || std::abs(delta) < 1e-40)
2778  return;
2779 
2780  double rhok = delta / theta, sigma = theta / delta;
2781  for (unsigned int k = 0; k < data.degree - 1; ++k)
2782  {
2783  const double rhokp = 1. / (2. * sigma - rhok);
2784  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2785  rhok = rhokp;
2786  matrix_ptr->Tvmult(temp_vector1, solution);
2787  internal::PreconditionChebyshevImplementation::vector_updates(
2788  rhs,
2790  k + 2,
2791  factor1,
2792  factor2,
2793  solution_old,
2794  temp_vector1,
2795  temp_vector2,
2796  solution);
2797  }
2798 }
2799 
2800 
2801 
2802 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2803 inline typename PreconditionChebyshev<MatrixType,
2804  VectorType,
2805  PreconditionerType>::size_type
2807 {
2808  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2809  return matrix_ptr->m();
2810 }
2811 
2812 
2813 
2814 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2815 inline typename PreconditionChebyshev<MatrixType,
2816  VectorType,
2817  PreconditionerType>::size_type
2819 {
2820  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2821  return matrix_ptr->n();
2822 }
2823 
2824 #endif // DOXYGEN
2825 
2827 
2828 #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:801
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.)
static const char U
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:1571
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)
typename PreconditionRelaxation< MatrixType >::size_type size_type
Definition: precondition.h:765
AdditionalData(const double relaxation=1.)
#define AssertCudaKernel()
Definition: exceptions.h:1859
size_type m() const
size_type locally_owned_size() const
const std::vector< size_type > * permutation
Definition: precondition.h:854
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:797
#define Assert(cond, exc)
Definition: exceptions.h:1461
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:672
void step(VectorType &x, const VectorType &rhs) const
typename PreconditionRelaxation< MatrixType >::size_type size_type
Definition: precondition.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
VectorType::value_type * end(VectorType &V)
void vmult_add(VectorType &, const VectorType &) const
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:50
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
typename PreconditionRelaxation< MatrixType >::size_type size_type
Definition: precondition.h:594
void Tstep(VectorType &x, const VectorType &rhs) const
typename PreconditionRelaxation< MatrixType >::size_type size_type
Definition: precondition.h:666
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:400
VectorType::value_type * begin(VectorType &V)
const std::vector< size_type > * inverse_permutation
Definition: precondition.h:858
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:509
void vmult(VectorType &, const VectorType &) const
const std::vector< size_type > & permutation
Definition: precondition.h:793
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:600
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:725
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1885
bool is_empty() const
Definition: index_set.h:1829
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:138
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())