Reference documentation for deal.II version Git 040c6ad7d4 2020-09-26 18:01:03 +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 
78 {
79 public:
84 
90  {
94  AdditionalData() = default;
95  };
96 
101 
106  template <typename MatrixType>
107  void
108  initialize(const MatrixType & matrix,
109  const AdditionalData &additional_data = AdditionalData());
110 
114  template <class VectorType>
115  void
116  vmult(VectorType &, const VectorType &) const;
117 
122  template <class VectorType>
123  void
124  Tvmult(VectorType &, const VectorType &) const;
125 
129  template <class VectorType>
130  void
131  vmult_add(VectorType &, const VectorType &) const;
132 
137  template <class VectorType>
138  void
139  Tvmult_add(VectorType &, const VectorType &) const;
140 
145  void
147  {}
148 
156  size_type
157  m() const;
158 
166  size_type
167  n() const;
168 
169 private:
174 
179 };
180 
181 
182 
194 {
195 public:
200 
205  {
206  public:
211  AdditionalData(const double relaxation = 1.);
212 
216  double relaxation;
217  };
218 
224 
228  void
229  initialize(const AdditionalData &parameters);
230 
236  template <typename MatrixType>
237  void
238  initialize(const MatrixType &matrix, const AdditionalData &parameters);
239 
243  template <class VectorType>
244  void
245  vmult(VectorType &, const VectorType &) const;
246 
251  template <class VectorType>
252  void
253  Tvmult(VectorType &, const VectorType &) const;
257  template <class VectorType>
258  void
259  vmult_add(VectorType &, const VectorType &) const;
260 
265  template <class VectorType>
266  void
267  Tvmult_add(VectorType &, const VectorType &) const;
268 
273  void
275  {}
276 
284  size_type
285  m() const;
286 
294  size_type
295  n() const;
296 
297 private:
301  double relaxation;
302 
307 
312 };
313 
314 
315 
355 template <typename MatrixType = SparseMatrix<double>,
356  class VectorType = Vector<double>>
358 {
359 public:
363  using function_ptr = void (MatrixType::*)(VectorType &,
364  const VectorType &) const;
365 
371  PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
372 
377  void
378  vmult(VectorType &dst, const VectorType &src) const;
379 
380 private:
384  const MatrixType &matrix;
385 
390 };
391 
392 
393 
399 template <typename MatrixType = SparseMatrix<double>>
401 {
402 public:
407 
412  {
413  public:
417  AdditionalData(const double relaxation = 1.);
418 
422  double relaxation;
423  };
424 
430  void
431  initialize(const MatrixType & A,
432  const AdditionalData &parameters = AdditionalData());
433 
437  void
438  clear();
439 
444  size_type
445  m() const;
446 
451  size_type
452  n() const;
453 
454 protected:
459 
463  double relaxation;
464 };
465 
466 
467 
494 template <typename MatrixType = SparseMatrix<double>>
495 class PreconditionJacobi : public PreconditionRelaxation<MatrixType>
496 {
497 public:
501  using AdditionalData =
503 
507  template <class VectorType>
508  void
509  vmult(VectorType &, const VectorType &) const;
510 
515  template <class VectorType>
516  void
517  Tvmult(VectorType &, const VectorType &) const;
518 
522  template <class VectorType>
523  void
524  step(VectorType &x, const VectorType &rhs) const;
525 
529  template <class VectorType>
530  void
531  Tstep(VectorType &x, const VectorType &rhs) const;
532 };
533 
534 
580 template <typename MatrixType = SparseMatrix<double>>
581 class PreconditionSOR : public PreconditionRelaxation<MatrixType>
582 {
583 public:
587  using AdditionalData =
589 
593  template <class VectorType>
594  void
595  vmult(VectorType &, const VectorType &) const;
596 
600  template <class VectorType>
601  void
602  Tvmult(VectorType &, const VectorType &) const;
603 
607  template <class VectorType>
608  void
609  step(VectorType &x, const VectorType &rhs) const;
610 
614  template <class VectorType>
615  void
616  Tstep(VectorType &x, const VectorType &rhs) const;
617 };
618 
619 
620 
647 template <typename MatrixType = SparseMatrix<double>>
648 class PreconditionSSOR : public PreconditionRelaxation<MatrixType>
649 {
650 public:
654  using AdditionalData =
656 
661 
666 
667 
673  void
674  initialize(const MatrixType & A,
675  const typename BaseClass::AdditionalData &parameters =
676  typename BaseClass::AdditionalData());
677 
681  template <class VectorType>
682  void
683  vmult(VectorType &, const VectorType &) const;
684 
689  template <class VectorType>
690  void
691  Tvmult(VectorType &, const VectorType &) const;
692 
693 
697  template <class VectorType>
698  void
699  step(VectorType &x, const VectorType &rhs) const;
700 
704  template <class VectorType>
705  void
706  Tstep(VectorType &x, const VectorType &rhs) const;
707 
708 private:
713  std::vector<std::size_t> pos_right_of_diagonal;
714 };
715 
716 
746 template <typename MatrixType = SparseMatrix<double>>
747 class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
748 {
749 public:
754 
759  {
760  public:
772  const std::vector<size_type> &permutation,
773  const std::vector<size_type> &inverse_permutation,
775  &parameters =
777 
781  const std::vector<size_type> &permutation;
785  const std::vector<size_type> &inverse_permutation;
790  };
791 
803  void
804  initialize(const MatrixType & A,
805  const std::vector<size_type> &permutation,
806  const std::vector<size_type> &inverse_permutation,
808  &parameters =
810 
821  void
822  initialize(const MatrixType &A, const AdditionalData &additional_data);
823 
827  template <class VectorType>
828  void
829  vmult(VectorType &, const VectorType &) const;
830 
834  template <class VectorType>
835  void
836  Tvmult(VectorType &, const VectorType &) const;
837 
838 private:
842  const std::vector<size_type> *permutation;
846  const std::vector<size_type> *inverse_permutation;
847 };
848 
849 
850 
986 template <typename MatrixType = SparseMatrix<double>,
987  typename VectorType = Vector<double>,
988  typename PreconditionerType = DiagonalMatrix<VectorType>>
990 {
991 public:
996 
1002  {
1006  AdditionalData(const unsigned int degree = 1,
1007  const double smoothing_range = 0.,
1008  const unsigned int eig_cg_n_iterations = 8,
1009  const double eig_cg_residual = 1e-2,
1010  const double max_eigenvalue = 1);
1011 
1015  AdditionalData &
1016  operator=(const AdditionalData &other_data);
1017 
1029  unsigned int degree;
1030 
1043 
1050  unsigned int eig_cg_n_iterations;
1051 
1057 
1064 
1070 
1074  std::shared_ptr<PreconditionerType> preconditioner;
1075  };
1076 
1077 
1082 
1094  void
1095  initialize(const MatrixType & matrix,
1096  const AdditionalData &additional_data = AdditionalData());
1097 
1102  void
1103  vmult(VectorType &dst, const VectorType &src) const;
1104 
1109  void
1110  Tvmult(VectorType &dst, const VectorType &src) const;
1111 
1115  void
1116  step(VectorType &dst, const VectorType &src) const;
1117 
1121  void
1122  Tstep(VectorType &dst, const VectorType &src) const;
1123 
1127  void
1128  clear();
1129 
1134  size_type
1135  m() const;
1136 
1141  size_type
1142  n() const;
1143 
1149  {
1161  unsigned int cg_iterations;
1166  unsigned int degree;
1171  : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
1172  , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
1173  , cg_iterations{0}
1174  , degree{0}
1175  {}
1176  };
1177 
1191  estimate_eigenvalues(const VectorType &src) const;
1192 
1193 private:
1197  SmartPointer<
1198  const MatrixType,
1201 
1206 
1211 
1216 
1222 
1226  double theta;
1227 
1232  double delta;
1233 
1239 
1245 };
1246 
1247 
1248 
1250 /* ---------------------------------- Inline functions ------------------- */
1251 
1252 #ifndef DOXYGEN
1253 
1255  : n_rows(0)
1256  , n_columns(0)
1257 {}
1258 
1259 template <typename MatrixType>
1260 inline void
1261 PreconditionIdentity::initialize(const MatrixType &matrix,
1263 {
1264  n_rows = matrix.m();
1265  n_columns = matrix.n();
1266 }
1267 
1268 
1269 template <class VectorType>
1270 inline void
1271 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1272 {
1273  dst = src;
1274 }
1275 
1276 
1277 
1278 template <class VectorType>
1279 inline void
1280 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1281 {
1282  dst = src;
1283 }
1284 
1285 template <class VectorType>
1286 inline void
1288 {
1289  dst += src;
1290 }
1291 
1292 
1293 
1294 template <class VectorType>
1295 inline void
1297 {
1298  dst += src;
1299 }
1300 
1303 {
1304  Assert(n_rows != 0, ExcNotInitialized());
1305  return n_rows;
1306 }
1307 
1310 {
1311  Assert(n_columns != 0, ExcNotInitialized());
1312  return n_columns;
1313 }
1314 
1315 //---------------------------------------------------------------------------
1316 
1318  const double relaxation)
1319  : relaxation(relaxation)
1320 {}
1321 
1322 
1324  : relaxation(0)
1325  , n_rows(0)
1326  , n_columns(0)
1327 {
1328  AdditionalData add_data;
1329  relaxation = add_data.relaxation;
1330 }
1331 
1332 
1333 
1334 inline void
1336  const PreconditionRichardson::AdditionalData &parameters)
1337 {
1338  relaxation = parameters.relaxation;
1339 }
1340 
1341 
1342 
1343 template <typename MatrixType>
1344 inline void
1346  const MatrixType & matrix,
1347  const PreconditionRichardson::AdditionalData &parameters)
1348 {
1349  relaxation = parameters.relaxation;
1350  n_rows = matrix.m();
1351  n_columns = matrix.n();
1352 }
1353 
1354 
1355 
1356 template <class VectorType>
1357 inline void
1358 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
1359 {
1360  static_assert(
1361  std::is_same<size_type, typename VectorType::size_type>::value,
1362  "PreconditionRichardson and VectorType must have the same size_type.");
1363 
1364  dst.equ(relaxation, src);
1365 }
1366 
1367 
1368 
1369 template <class VectorType>
1370 inline void
1372 {
1373  static_assert(
1374  std::is_same<size_type, typename VectorType::size_type>::value,
1375  "PreconditionRichardson and VectorType must have the same size_type.");
1376 
1377  dst.equ(relaxation, src);
1378 }
1379 
1380 template <class VectorType>
1381 inline void
1383 {
1384  static_assert(
1385  std::is_same<size_type, typename VectorType::size_type>::value,
1386  "PreconditionRichardson and VectorType must have the same size_type.");
1387 
1388  dst.add(relaxation, src);
1389 }
1390 
1391 
1392 
1393 template <class VectorType>
1394 inline void
1396 {
1397  static_assert(
1398  std::is_same<size_type, typename VectorType::size_type>::value,
1399  "PreconditionRichardson and VectorType must have the same size_type.");
1400 
1401  dst.add(relaxation, src);
1402 }
1403 
1406 {
1407  Assert(n_rows != 0, ExcNotInitialized());
1408  return n_rows;
1409 }
1410 
1413 {
1414  Assert(n_columns != 0, ExcNotInitialized());
1415  return n_columns;
1416 }
1417 
1418 //---------------------------------------------------------------------------
1419 
1420 template <typename MatrixType>
1421 inline void
1423  const AdditionalData &parameters)
1424 {
1425  A = &rA;
1426  relaxation = parameters.relaxation;
1427 }
1428 
1429 
1430 template <typename MatrixType>
1431 inline void
1433 {
1434  A = nullptr;
1435 }
1436 
1437 template <typename MatrixType>
1440 {
1441  Assert(A != nullptr, ExcNotInitialized());
1442  return A->m();
1443 }
1444 
1445 template <typename MatrixType>
1448 {
1449  Assert(A != nullptr, ExcNotInitialized());
1450  return A->n();
1451 }
1452 
1453 //---------------------------------------------------------------------------
1454 
1455 template <typename MatrixType>
1456 template <class VectorType>
1457 inline void
1459  const VectorType &src) const
1460 {
1461  static_assert(
1462  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1463  typename VectorType::size_type>::value,
1464  "PreconditionJacobi and VectorType must have the same size_type.");
1465 
1466  Assert(this->A != nullptr, ExcNotInitialized());
1467  this->A->precondition_Jacobi(dst, src, this->relaxation);
1468 }
1469 
1470 
1471 
1472 template <typename MatrixType>
1473 template <class VectorType>
1474 inline void
1476  const VectorType &src) const
1477 {
1478  static_assert(
1479  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1480  typename VectorType::size_type>::value,
1481  "PreconditionJacobi and VectorType must have the same size_type.");
1482 
1483  Assert(this->A != nullptr, ExcNotInitialized());
1484  this->A->precondition_Jacobi(dst, src, this->relaxation);
1485 }
1486 
1487 
1488 
1489 template <typename MatrixType>
1490 template <class VectorType>
1491 inline void
1493  const VectorType &src) const
1494 {
1495  static_assert(
1496  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1497  typename VectorType::size_type>::value,
1498  "PreconditionJacobi and VectorType must have the same size_type.");
1499 
1500  Assert(this->A != nullptr, ExcNotInitialized());
1501  this->A->Jacobi_step(dst, src, this->relaxation);
1502 }
1503 
1504 
1505 
1506 template <typename MatrixType>
1507 template <class VectorType>
1508 inline void
1510  const VectorType &src) const
1511 {
1512  static_assert(
1513  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1514  typename VectorType::size_type>::value,
1515  "PreconditionJacobi and VectorType must have the same size_type.");
1516 
1517  step(dst, src);
1518 }
1519 
1520 
1521 
1522 //---------------------------------------------------------------------------
1523 
1524 template <typename MatrixType>
1525 template <class VectorType>
1526 inline void
1528 {
1529  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1530  typename VectorType::size_type>::value,
1531  "PreconditionSOR and VectorType must have the same size_type.");
1532 
1533  Assert(this->A != nullptr, ExcNotInitialized());
1534  this->A->precondition_SOR(dst, src, this->relaxation);
1535 }
1536 
1537 
1538 
1539 template <typename MatrixType>
1540 template <class VectorType>
1541 inline void
1543  const VectorType &src) const
1544 {
1545  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1546  typename VectorType::size_type>::value,
1547  "PreconditionSOR and VectorType must have the same size_type.");
1548 
1549  Assert(this->A != nullptr, ExcNotInitialized());
1550  this->A->precondition_TSOR(dst, src, this->relaxation);
1551 }
1552 
1553 
1554 
1555 template <typename MatrixType>
1556 template <class VectorType>
1557 inline void
1559 {
1560  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1561  typename VectorType::size_type>::value,
1562  "PreconditionSOR and VectorType must have the same size_type.");
1563 
1564  Assert(this->A != nullptr, ExcNotInitialized());
1565  this->A->SOR_step(dst, src, this->relaxation);
1566 }
1567 
1568 
1569 
1570 template <typename MatrixType>
1571 template <class VectorType>
1572 inline void
1574 {
1575  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1576  typename VectorType::size_type>::value,
1577  "PreconditionSOR and VectorType must have the same size_type.");
1578 
1579  Assert(this->A != nullptr, ExcNotInitialized());
1580  this->A->TSOR_step(dst, src, this->relaxation);
1581 }
1582 
1583 
1584 
1585 //---------------------------------------------------------------------------
1586 
1587 template <typename MatrixType>
1588 inline void
1590  const MatrixType & rA,
1591  const typename BaseClass::AdditionalData &parameters)
1592 {
1593  this->PreconditionRelaxation<MatrixType>::initialize(rA, parameters);
1594 
1595  // in case we have a SparseMatrix class, we can extract information about
1596  // the diagonal.
1598  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
1599  &*this->A);
1600 
1601  // calculate the positions first after the diagonal.
1602  if (mat != nullptr)
1603  {
1604  const size_type n = this->A->n();
1605  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1606  for (size_type row = 0; row < n; ++row)
1607  {
1608  // find the first element in this line which is on the right of the
1609  // diagonal. we need to precondition with the elements on the left
1610  // only. note: the first entry in each line denotes the diagonal
1611  // element, which we need not check.
1613  it = mat->begin(row) + 1;
1614  for (; it < mat->end(row); ++it)
1615  if (it->column() > row)
1616  break;
1617  pos_right_of_diagonal[row] = it - mat->begin();
1618  }
1619  }
1620 }
1621 
1622 
1623 template <typename MatrixType>
1624 template <class VectorType>
1625 inline void
1627  const VectorType &src) const
1628 {
1629  static_assert(
1630  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1631  typename VectorType::size_type>::value,
1632  "PreconditionSSOR and VectorType must have the same size_type.");
1633 
1634  Assert(this->A != nullptr, ExcNotInitialized());
1635  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1636 }
1637 
1638 
1639 
1640 template <typename MatrixType>
1641 template <class VectorType>
1642 inline void
1644  const VectorType &src) const
1645 {
1646  static_assert(
1647  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1648  typename VectorType::size_type>::value,
1649  "PreconditionSSOR and VectorType must have the same size_type.");
1650 
1651  Assert(this->A != nullptr, ExcNotInitialized());
1652  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1653 }
1654 
1655 
1656 
1657 template <typename MatrixType>
1658 template <class VectorType>
1659 inline void
1661 {
1662  static_assert(
1663  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1664  typename VectorType::size_type>::value,
1665  "PreconditionSSOR and VectorType must have the same size_type.");
1666 
1667  Assert(this->A != nullptr, ExcNotInitialized());
1668  this->A->SSOR_step(dst, src, this->relaxation);
1669 }
1670 
1671 
1672 
1673 template <typename MatrixType>
1674 template <class VectorType>
1675 inline void
1677  const VectorType &src) const
1678 {
1679  static_assert(
1680  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1681  typename VectorType::size_type>::value,
1682  "PreconditionSSOR and VectorType must have the same size_type.");
1683 
1684  step(dst, src);
1685 }
1686 
1687 
1688 
1689 //---------------------------------------------------------------------------
1690 
1691 template <typename MatrixType>
1692 inline void
1694  const MatrixType & rA,
1695  const std::vector<size_type> & p,
1696  const std::vector<size_type> & ip,
1697  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1698 {
1699  permutation = &p;
1700  inverse_permutation = &ip;
1702 }
1703 
1704 
1705 template <typename MatrixType>
1706 inline void
1707 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
1708  const AdditionalData &additional_data)
1709 {
1710  initialize(A,
1711  additional_data.permutation,
1712  additional_data.inverse_permutation,
1713  additional_data.parameters);
1714 }
1715 
1716 
1717 template <typename MatrixType>
1718 template <typename VectorType>
1719 inline void
1721  const VectorType &src) const
1722 {
1723  static_assert(
1724  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1725  typename VectorType::size_type>::value,
1726  "PreconditionPSOR and VectorType must have the same size_type.");
1727 
1728  Assert(this->A != nullptr, ExcNotInitialized());
1729  dst = src;
1730  this->A->PSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1731 }
1732 
1733 
1734 
1735 template <typename MatrixType>
1736 template <class VectorType>
1737 inline void
1739  const VectorType &src) const
1740 {
1741  static_assert(
1742  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1743  typename VectorType::size_type>::value,
1744  "PreconditionPSOR and VectorType must have the same size_type.");
1745 
1746  Assert(this->A != nullptr, ExcNotInitialized());
1747  dst = src;
1748  this->A->TPSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1749 }
1750 
1751 template <typename MatrixType>
1753  const std::vector<size_type> &permutation,
1754  const std::vector<size_type> &inverse_permutation,
1755  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1756  : permutation(permutation)
1757  , inverse_permutation(inverse_permutation)
1758  , parameters(parameters)
1759 {}
1760 
1761 
1762 //---------------------------------------------------------------------------
1763 
1764 
1765 template <typename MatrixType, class VectorType>
1767  const MatrixType & M,
1768  const function_ptr method)
1769  : matrix(M)
1770  , precondition(method)
1771 {}
1772 
1773 
1774 
1775 template <typename MatrixType, class VectorType>
1776 void
1778  VectorType & dst,
1779  const VectorType &src) const
1780 {
1781  (matrix.*precondition)(dst, src);
1782 }
1783 
1784 //---------------------------------------------------------------------------
1785 
1786 template <typename MatrixType>
1788  const double relaxation)
1789  : relaxation(relaxation)
1790 {}
1791 
1792 
1793 
1794 //---------------------------------------------------------------------------
1795 
1796 namespace internal
1797 {
1798  namespace PreconditionChebyshevImplementation
1799  {
1800  // for deal.II vectors, perform updates for Chebyshev preconditioner all
1801  // at once to reduce memory transfer. Here, we select between general
1802  // vectors and deal.II vectors where we expand the loop over the (local)
1803  // size of the vector
1804 
1805  // generic part for non-deal.II vectors
1806  template <typename VectorType, typename PreconditionerType>
1807  inline void
1808  vector_updates(const VectorType & rhs,
1809  const PreconditionerType &preconditioner,
1810  const unsigned int iteration_index,
1811  const double factor1,
1812  const double factor2,
1816  VectorType & solution)
1817  {
1818  if (iteration_index == 0)
1819  {
1820  solution.equ(factor2, rhs);
1821  preconditioner.vmult(solution_old, solution);
1822  }
1823  else if (iteration_index == 1)
1824  {
1825  // compute t = P^{-1} * (b-A*x^{n})
1826  temp_vector1.sadd(-1.0, 1.0, rhs);
1827  preconditioner.vmult(solution_old, temp_vector1);
1828 
1829  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
1830  solution_old.sadd(factor2, 1 + factor1, solution);
1831  }
1832  else
1833  {
1834  // compute t = P^{-1} * (b-A*x^{n})
1835  temp_vector1.sadd(-1.0, 1.0, rhs);
1836  preconditioner.vmult(temp_vector2, temp_vector1);
1837 
1838  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
1839  solution_old.sadd(-factor1, factor2, temp_vector2);
1840  solution_old.add(1 + factor1, solution);
1841  }
1842 
1843  solution.swap(solution_old);
1844  }
1845 
1846  // worker routine for deal.II vectors. Because of vectorization, we need
1847  // to put the loop into an extra structure because the virtual function of
1848  // VectorUpdatesRange prevents the compiler from applying vectorization.
1849  template <typename Number>
1850  struct VectorUpdater
1851  {
1852  VectorUpdater(const Number * rhs,
1853  const Number * matrix_diagonal_inverse,
1854  const unsigned int iteration_index,
1855  const Number factor1,
1856  const Number factor2,
1857  Number * solution_old,
1858  Number * tmp_vector,
1859  Number * solution)
1860  : rhs(rhs)
1861  , matrix_diagonal_inverse(matrix_diagonal_inverse)
1862  , iteration_index(iteration_index)
1863  , factor1(factor1)
1864  , factor2(factor2)
1865  , solution_old(solution_old)
1866  , tmp_vector(tmp_vector)
1867  , solution(solution)
1868  {}
1869 
1870  void
1871  apply_to_subrange(const std::size_t begin, const std::size_t end) const
1872  {
1873  // To circumvent a bug in gcc
1874  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
1875  // copies of the variables factor1 and factor2 and do not check based on
1876  // factor1.
1877  const Number factor1 = this->factor1;
1878  const Number factor1_plus_1 = 1. + this->factor1;
1879  const Number factor2 = this->factor2;
1880  if (iteration_index == 0)
1881  {
1883  for (std::size_t i = begin; i < end; ++i)
1884  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
1885  }
1886  else if (iteration_index == 1)
1887  {
1888  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
1890  for (std::size_t i = begin; i < end; ++i)
1891  // for efficiency reason, write back to temp_vector that is
1892  // already read (avoid read-for-ownership)
1893  tmp_vector[i] =
1894  factor1_plus_1 * solution[i] +
1895  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1896  }
1897  else
1898  {
1899  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
1900  // + f_2 * P^{-1} * (b-A*x^{n})
1902  for (std::size_t i = begin; i < end; ++i)
1903  solution_old[i] =
1904  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
1905  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1906  }
1907  }
1908 
1909  const Number * rhs;
1910  const Number * matrix_diagonal_inverse;
1911  const unsigned int iteration_index;
1912  const Number factor1;
1913  const Number factor2;
1914  mutable Number * solution_old;
1915  mutable Number * tmp_vector;
1916  mutable Number * solution;
1917  };
1918 
1919  template <typename Number>
1920  struct VectorUpdatesRange : public ::parallel::ParallelForInteger
1921  {
1922  VectorUpdatesRange(const VectorUpdater<Number> &updater,
1923  const std::size_t size)
1924  : updater(updater)
1925  {
1927  VectorUpdatesRange::apply_to_subrange(0, size);
1928  else
1929  apply_parallel(
1930  0,
1931  size,
1933  }
1934 
1935  ~VectorUpdatesRange() override = default;
1936 
1937  virtual void
1938  apply_to_subrange(const std::size_t begin,
1939  const std::size_t end) const override
1940  {
1941  updater.apply_to_subrange(begin, end);
1942  }
1943 
1944  const VectorUpdater<Number> &updater;
1945  };
1946 
1947  // selection for diagonal matrix around deal.II vector
1948  template <typename Number>
1949  inline void
1950  vector_updates(const ::Vector<Number> & rhs,
1951  const DiagonalMatrix<::Vector<Number>> &jacobi,
1952  const unsigned int iteration_index,
1953  const double factor1,
1954  const double factor2,
1955  ::Vector<Number> &solution_old,
1956  ::Vector<Number> &temp_vector1,
1957  ::Vector<Number> &,
1958  ::Vector<Number> &solution)
1959  {
1960  VectorUpdater<Number> upd(rhs.begin(),
1961  jacobi.get_vector().begin(),
1962  iteration_index,
1963  factor1,
1964  factor2,
1965  solution_old.begin(),
1966  temp_vector1.begin(),
1967  solution.begin());
1968  VectorUpdatesRange<Number>(upd, rhs.size());
1969 
1970  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
1971  if (iteration_index == 0)
1972  {
1973  // nothing to do here because we can immediately write into the
1974  // solution vector without remembering any of the other vectors
1975  }
1976  else if (iteration_index == 1)
1977  {
1978  solution.swap(temp_vector1);
1979  solution_old.swap(temp_vector1);
1980  }
1981  else
1982  solution.swap(solution_old);
1983  }
1984 
1985  // selection for diagonal matrix around parallel deal.II vector
1986  template <typename Number>
1987  inline void
1988  vector_updates(
1990  const DiagonalMatrix<
1992  const unsigned int iteration_index,
1993  const double factor1,
1994  const double factor2,
1996  &solution_old,
1998  &temp_vector1,
2001  {
2002  VectorUpdater<Number> upd(rhs.begin(),
2003  jacobi.get_vector().begin(),
2004  iteration_index,
2005  factor1,
2006  factor2,
2007  solution_old.begin(),
2008  temp_vector1.begin(),
2009  solution.begin());
2010  VectorUpdatesRange<Number>(upd, rhs.local_size());
2011 
2012  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2013  if (iteration_index == 0)
2014  {
2015  // nothing to do here because we can immediately write into the
2016  // solution vector without remembering any of the other vectors
2017  }
2018  else if (iteration_index == 1)
2019  {
2020  solution.swap(temp_vector1);
2021  solution_old.swap(temp_vector1);
2022  }
2023  else
2024  solution.swap(solution_old);
2025  }
2026 
2027  template <typename MatrixType, typename PreconditionerType>
2028  inline void
2029  initialize_preconditioner(
2030  const MatrixType & matrix,
2031  std::shared_ptr<PreconditionerType> &preconditioner)
2032  {
2033  (void)matrix;
2034  (void)preconditioner;
2035  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2036  }
2037 
2038  template <typename MatrixType, typename VectorType>
2039  inline void
2040  initialize_preconditioner(
2041  const MatrixType & matrix,
2042  std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
2043  {
2044  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2045  {
2046  if (preconditioner.get() == nullptr)
2047  preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2048 
2049  Assert(
2050  preconditioner->m() == 0,
2051  ExcMessage(
2052  "Preconditioner appears to be initialized but not sized correctly"));
2053 
2054  // This part only works in serial
2055  if (preconditioner->m() != matrix.m())
2056  {
2057  preconditioner->get_vector().reinit(matrix.m());
2058  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2059  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2060  }
2061  }
2062  }
2063 
2064  template <typename VectorType>
2065  void
2066  set_initial_guess(VectorType &vector)
2067  {
2068  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2069  if (vector.locally_owned_elements().is_element(0))
2070  vector(0) = 0.;
2071  }
2072 
2073  template <typename Number>
2074  void
2075  set_initial_guess(::Vector<Number> &vector)
2076  {
2077  // Choose a high-frequency mode consisting of numbers between 0 and 1
2078  // that is cheap to compute (cheaper than random numbers) but avoids
2079  // obviously re-occurring numbers in multi-component systems by choosing
2080  // a period of 11
2081  for (unsigned int i = 0; i < vector.size(); ++i)
2082  vector(i) = i % 11;
2083 
2084  const Number mean_value = vector.mean_value();
2085  vector.add(-mean_value);
2086  }
2087 
2088  template <typename Number>
2089  void
2090  set_initial_guess(
2092  &vector)
2093  {
2094  // Choose a high-frequency mode consisting of numbers between 0 and 1
2095  // that is cheap to compute (cheaper than random numbers) but avoids
2096  // obviously re-occurring numbers in multi-component systems by choosing
2097  // a period of 11.
2098  // Make initial guess robust with respect to number of processors
2099  // by operating on the global index.
2100  types::global_dof_index first_local_range = 0;
2101  if (!vector.locally_owned_elements().is_empty())
2102  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2103  for (unsigned int i = 0; i < vector.local_size(); ++i)
2104  vector.local_element(i) = (i + first_local_range) % 11;
2105 
2106  const Number mean_value = vector.mean_value();
2107  vector.add(-mean_value);
2108  }
2109 
2110 
2111 # ifdef DEAL_II_COMPILER_CUDA_AWARE
2112  template <typename Number>
2113  __global__ void
2114  set_initial_guess_kernel(const types::global_dof_index offset,
2115  const unsigned int local_size,
2116  Number * values)
2117 
2118  {
2119  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
2120  if (index < local_size)
2121  values[index] = (index + offset) % 11;
2122  }
2123 
2124  template <typename Number>
2125  void
2126  set_initial_guess(
2128  &vector)
2129  {
2130  // Choose a high-frequency mode consisting of numbers between 0 and 1
2131  // that is cheap to compute (cheaper than random numbers) but avoids
2132  // obviously re-occurring numbers in multi-component systems by choosing
2133  // a period of 11.
2134  // Make initial guess robust with respect to number of processors
2135  // by operating on the global index.
2136  types::global_dof_index first_local_range = 0;
2137  if (!vector.locally_owned_elements().is_empty())
2138  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2139 
2140  const auto n_local_elements = vector.local_size();
2141  const int n_blocks =
2142  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
2143  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
2144  first_local_range, n_local_elements, vector.get_values());
2145  AssertCudaKernel();
2146 
2147  const Number mean_value = vector.mean_value();
2148  vector.add(-mean_value);
2149  }
2150 # endif // DEAL_II_COMPILER_CUDA_AWARE
2151 
2152  struct EigenvalueTracker
2153  {
2154  public:
2155  void
2156  slot(const std::vector<double> &eigenvalues)
2157  {
2158  values = eigenvalues;
2159  }
2160 
2161  std::vector<double> values;
2162  };
2163  } // namespace PreconditionChebyshevImplementation
2164 } // namespace internal
2165 
2166 
2167 
2168 template <typename MatrixType, class VectorType, typename PreconditionerType>
2170  AdditionalData::AdditionalData(const unsigned int degree,
2171  const double smoothing_range,
2172  const unsigned int eig_cg_n_iterations,
2173  const double eig_cg_residual,
2174  const double max_eigenvalue)
2175  : degree(degree)
2176  , smoothing_range(smoothing_range)
2177  , eig_cg_n_iterations(eig_cg_n_iterations)
2178  , eig_cg_residual(eig_cg_residual)
2179  , max_eigenvalue(max_eigenvalue)
2180 {}
2181 
2182 
2183 
2184 template <typename MatrixType, class VectorType, typename PreconditionerType>
2185 inline typename PreconditionChebyshev<MatrixType,
2186  VectorType,
2187  PreconditionerType>::AdditionalData &
2189  AdditionalData::operator=(const AdditionalData &other_data)
2190 {
2191  degree = other_data.degree;
2192  smoothing_range = other_data.smoothing_range;
2193  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
2194  eig_cg_residual = other_data.eig_cg_residual;
2195  max_eigenvalue = other_data.max_eigenvalue;
2196  preconditioner = other_data.preconditioner;
2197  constraints.copy_from(other_data.constraints);
2198 
2199  return *this;
2200 }
2201 
2202 
2203 
2204 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2207  : theta(1.)
2208  , delta(1.)
2210 {
2211  static_assert(
2212  std::is_same<size_type, typename VectorType::size_type>::value,
2213  "PreconditionChebyshev and VectorType must have the same size_type.");
2214 }
2215 
2216 
2217 
2218 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2219 inline void
2221  const MatrixType & matrix,
2222  const AdditionalData &additional_data)
2223 {
2224  matrix_ptr = &matrix;
2225  data = additional_data;
2226  Assert(data.degree > 0,
2227  ExcMessage("The degree of the Chebyshev method must be positive."));
2228  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2229  matrix, data.preconditioner);
2231 }
2232 
2233 
2234 
2235 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2236 inline void
2238 {
2240  theta = delta = 1.0;
2241  matrix_ptr = nullptr;
2242  {
2243  VectorType empty_vector;
2244  solution_old.reinit(empty_vector);
2245  temp_vector1.reinit(empty_vector);
2246  temp_vector2.reinit(empty_vector);
2247  }
2248  data.preconditioner.reset();
2249 }
2250 
2251 
2252 
2253 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2254 inline typename PreconditionChebyshev<MatrixType,
2255  VectorType,
2256  PreconditionerType>::EigenvalueInformation
2258  estimate_eigenvalues(const VectorType &src) const
2259 {
2261  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2262 
2264  EigenvalueInformation info{};
2265 
2266  solution_old.reinit(src);
2267  temp_vector1.reinit(src, true);
2268 
2269  // calculate largest eigenvalue using a hand-tuned CG iteration on the
2270  // matrix weighted by its diagonal. we start with a vector that consists of
2271  // ones only, weighted by the length.
2272  if (data.eig_cg_n_iterations > 0)
2273  {
2275  ExcMessage(
2276  "Need to set at least two iterations to find eigenvalues."));
2277 
2278  // set a very strict tolerance to force at least two iterations
2279  ReductionControl control(
2281  std::sqrt(
2283  1e-10,
2284  false,
2285  false);
2286 
2287  internal::PreconditionChebyshevImplementation::EigenvalueTracker
2288  eigenvalue_tracker;
2289  SolverCG<VectorType> solver(control);
2290  solver.connect_eigenvalues_slot(
2291  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2292  eigenvalue_tracker.slot(eigenvalues);
2293  });
2294 
2295  // set an initial guess which is close to the constant vector but where
2296  // one entry is different to trigger high frequencies
2297  internal::PreconditionChebyshevImplementation::set_initial_guess(
2298  temp_vector1);
2300 
2301  try
2302  {
2303  solver.solve(*matrix_ptr,
2304  solution_old,
2305  temp_vector1,
2306  *data.preconditioner);
2307  }
2309  {}
2310 
2311  // read the eigenvalues from the attached eigenvalue tracker
2312  if (eigenvalue_tracker.values.empty())
2313  info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
2314  else
2315  {
2316  info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
2317 
2318  // include a safety factor since the CG method will in general not
2319  // be converged
2320  info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
2321  }
2322 
2323  info.cg_iterations = control.last_step();
2324  }
2325  else
2326  {
2327  info.max_eigenvalue_estimate = data.max_eigenvalue;
2328  info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
2329  }
2330 
2331  const double alpha = (data.smoothing_range > 1. ?
2332  info.max_eigenvalue_estimate / data.smoothing_range :
2333  std::min(0.9 * info.max_eigenvalue_estimate,
2334  info.min_eigenvalue_estimate));
2335 
2336  // in case the user set the degree to invalid unsigned int, we have to
2337  // determine the number of necessary iterations from the Chebyshev error
2338  // estimate, given the target tolerance specified by smoothing_range. This
2339  // estimate is based on the error formula given in section 5.1 of
2340  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
2342  {
2343  const double actual_range = info.max_eigenvalue_estimate / alpha;
2344  const double sigma = (1. - std::sqrt(1. / actual_range)) /
2345  (1. + std::sqrt(1. / actual_range));
2346  const double eps = data.smoothing_range;
2347  const_cast<
2349  this)
2350  ->data.degree =
2351  1 + static_cast<unsigned int>(
2352  std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
2353  std::log(1. / sigma));
2354 
2355  info.degree = data.degree;
2356  }
2357 
2358  const_cast<
2360  ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
2361  const_cast<
2363  ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
2364 
2365  // We do not need the second temporary vector in case we have a
2366  // DiagonalMatrix as preconditioner and use deal.II's own vectors
2367  using NumberType = typename VectorType::value_type;
2368  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
2369  false ||
2370  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
2371  ((std::is_same<VectorType,
2373  Vector<NumberType, MemorySpace::Host>>::value ==
2374  false) ||
2375  (std::is_same<VectorType,
2376  LinearAlgebra::distributed::
2377  Vector<NumberType, MemorySpace::CUDA>>::value ==
2378  false))))
2379  temp_vector2.reinit(src, true);
2380  else
2381  {
2382  VectorType empty_vector;
2383  temp_vector2.reinit(empty_vector);
2384  }
2385 
2386  const_cast<
2388  ->eigenvalues_are_initialized = true;
2389 
2390  return info;
2391 }
2392 
2393 
2394 
2395 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2396 inline void
2398  VectorType & solution,
2399  const VectorType &rhs) const
2400 {
2401  std::lock_guard<std::mutex> lock(mutex);
2402  if (eigenvalues_are_initialized == false)
2403  estimate_eigenvalues(rhs);
2404 
2405  internal::PreconditionChebyshevImplementation::vector_updates(
2406  rhs,
2408  0,
2409  0.,
2410  1. / theta,
2411  solution_old,
2412  temp_vector1,
2413  temp_vector2,
2414  solution);
2415 
2416  // if delta is zero, we do not need to iterate because the updates will be
2417  // zero
2418  if (data.degree < 2 || std::abs(delta) < 1e-40)
2419  return;
2420 
2421  double rhok = delta / theta, sigma = theta / delta;
2422  for (unsigned int k = 0; k < data.degree - 1; ++k)
2423  {
2424  matrix_ptr->vmult(temp_vector1, solution);
2425  const double rhokp = 1. / (2. * sigma - rhok);
2426  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2427  rhok = rhokp;
2428  internal::PreconditionChebyshevImplementation::vector_updates(
2429  rhs,
2431  k + 1,
2432  factor1,
2433  factor2,
2434  solution_old,
2435  temp_vector1,
2436  temp_vector2,
2437  solution);
2438  }
2439 }
2440 
2441 
2442 
2443 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2444 inline void
2446  VectorType & solution,
2447  const VectorType &rhs) const
2448 {
2449  std::lock_guard<std::mutex> lock(mutex);
2450  if (eigenvalues_are_initialized == false)
2451  estimate_eigenvalues(rhs);
2452 
2453  internal::PreconditionChebyshevImplementation::vector_updates(
2454  rhs,
2456  0,
2457  0.,
2458  1. / theta,
2459  solution_old,
2460  temp_vector1,
2461  temp_vector2,
2462  solution);
2463 
2464  if (data.degree < 2 || std::abs(delta) < 1e-40)
2465  return;
2466 
2467  double rhok = delta / theta, sigma = theta / delta;
2468  for (unsigned int k = 0; k < data.degree - 1; ++k)
2469  {
2470  matrix_ptr->Tvmult(temp_vector1, solution);
2471  const double rhokp = 1. / (2. * sigma - rhok);
2472  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2473  rhok = rhokp;
2474  internal::PreconditionChebyshevImplementation::vector_updates(
2475  rhs,
2477  k + 1,
2478  factor1,
2479  factor2,
2480  solution_old,
2481  temp_vector1,
2482  temp_vector2,
2483  solution);
2484  }
2485 }
2486 
2487 
2488 
2489 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2490 inline void
2492  VectorType & solution,
2493  const VectorType &rhs) const
2494 {
2495  std::lock_guard<std::mutex> lock(mutex);
2496  if (eigenvalues_are_initialized == false)
2497  estimate_eigenvalues(rhs);
2498 
2499  matrix_ptr->vmult(temp_vector1, solution);
2500  internal::PreconditionChebyshevImplementation::vector_updates(
2501  rhs,
2503  1,
2504  0.,
2505  1. / theta,
2506  solution_old,
2507  temp_vector1,
2508  temp_vector2,
2509  solution);
2510 
2511  if (data.degree < 2 || std::abs(delta) < 1e-40)
2512  return;
2513 
2514  double rhok = delta / theta, sigma = theta / delta;
2515  for (unsigned int k = 0; k < data.degree - 1; ++k)
2516  {
2517  matrix_ptr->vmult(temp_vector1, solution);
2518  const double rhokp = 1. / (2. * sigma - rhok);
2519  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2520  rhok = rhokp;
2521  internal::PreconditionChebyshevImplementation::vector_updates(
2522  rhs,
2524  k + 2,
2525  factor1,
2526  factor2,
2527  solution_old,
2528  temp_vector1,
2529  temp_vector2,
2530  solution);
2531  }
2532 }
2533 
2534 
2535 
2536 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2537 inline void
2539  VectorType & solution,
2540  const VectorType &rhs) const
2541 {
2542  std::lock_guard<std::mutex> lock(mutex);
2543  if (eigenvalues_are_initialized == false)
2544  estimate_eigenvalues(rhs);
2545 
2546  matrix_ptr->Tvmult(temp_vector1, solution);
2547  internal::PreconditionChebyshevImplementation::vector_updates(
2548  rhs,
2550  1,
2551  0.,
2552  1. / theta,
2553  solution_old,
2554  temp_vector1,
2555  temp_vector2,
2556  solution);
2557 
2558  if (data.degree < 2 || std::abs(delta) < 1e-40)
2559  return;
2560 
2561  double rhok = delta / theta, sigma = theta / delta;
2562  for (unsigned int k = 0; k < data.degree - 1; ++k)
2563  {
2564  matrix_ptr->Tvmult(temp_vector1, solution);
2565  const double rhokp = 1. / (2. * sigma - rhok);
2566  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2567  rhok = rhokp;
2568  internal::PreconditionChebyshevImplementation::vector_updates(
2569  rhs,
2571  k + 2,
2572  factor1,
2573  factor2,
2574  solution_old,
2575  temp_vector1,
2576  temp_vector2,
2577  solution);
2578  }
2579 }
2580 
2581 
2582 
2583 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2584 inline typename PreconditionChebyshev<MatrixType,
2585  VectorType,
2586  PreconditionerType>::size_type
2588 {
2589  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2590  return matrix_ptr->m();
2591 }
2592 
2593 
2594 
2595 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2596 inline typename PreconditionChebyshev<MatrixType,
2597  VectorType,
2598  PreconditionerType>::size_type
2600 {
2601  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2602  return matrix_ptr->n();
2603 }
2604 
2605 #endif // DOXYGEN
2606 
2608 
2609 #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:789
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:384
AdditionalData data
AdditionalData(const double relaxation=1.)
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition: precondition.h:364
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:1521
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:1773
size_type m() const
const std::vector< size_type > * permutation
Definition: precondition.h:842
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
const function_ptr precondition
Definition: precondition.h:389
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:406
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:785
#define Assert(cond, exc)
Definition: exceptions.h:1411
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:655
void step(VectorType &x, const VectorType &rhs) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
VectorType::value_type * end(VectorType &V)
typename MatrixType::size_type size_type
Definition: precondition.h:753
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:458
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:362
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:660
const std::vector< size_type > * inverse_permutation
Definition: precondition.h:846
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:502
void vmult(VectorType &, const VectorType &) const
const std::vector< size_type > & permutation
Definition: precondition.h:781
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:588
Eigenvalue vector is filled.
virtual Number mean_value() const override
void vmult_add(VectorType &, const VectorType &) const
std::vector< std::size_t > pos_right_of_diagonal
Definition: precondition.h:713
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1878
bool is_empty() const
Definition: index_set.h:1822
T max(const T &t, const MPI_Comm &mpi_communicator)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:134
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())