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