Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 
1005 template <typename MatrixType = SparseMatrix<double>,
1006  typename VectorType = Vector<double>,
1007  typename PreconditionerType = DiagonalMatrix<VectorType>>
1009 {
1010 public:
1015 
1016  // avoid warning about use of deprecated variables
1017 
1023  {
1027  AdditionalData(const unsigned int degree = 1,
1028  const double smoothing_range = 0.,
1029  const unsigned int eig_cg_n_iterations = 8,
1030  const double eig_cg_residual = 1e-2,
1031  const double max_eigenvalue = 1);
1032 
1036  AdditionalData &
1037  operator=(const AdditionalData &other_data);
1038 
1050  unsigned int degree;
1051 
1064 
1071  unsigned int eig_cg_n_iterations;
1072 
1078 
1085 
1091 
1095  std::shared_ptr<PreconditionerType> preconditioner;
1096  };
1097 
1098 
1100 
1112  void
1113  initialize(const MatrixType & matrix,
1114  const AdditionalData &additional_data = AdditionalData());
1115 
1120  void
1121  vmult(VectorType &dst, const VectorType &src) const;
1122 
1127  void
1128  Tvmult(VectorType &dst, const VectorType &src) const;
1129 
1133  void
1134  step(VectorType &dst, const VectorType &src) const;
1135 
1139  void
1140  Tstep(VectorType &dst, const VectorType &src) const;
1141 
1145  void
1146  clear();
1147 
1152  size_type
1153  m() const;
1154 
1159  size_type
1160  n() const;
1161 
1162 private:
1166  SmartPointer<
1167  const MatrixType,
1170 
1174  mutable VectorType solution_old;
1175 
1179  mutable VectorType temp_vector1;
1180 
1184  mutable VectorType temp_vector2;
1185 
1191 
1195  double theta;
1196 
1201  double delta;
1202 
1208 
1214 
1221  void
1222  estimate_eigenvalues(const VectorType &src) const;
1223 };
1224 
1225 
1226 
1228 /* ---------------------------------- Inline functions ------------------- */
1229 
1230 #ifndef DOXYGEN
1231 
1233  : n_rows(0)
1234  , n_columns(0)
1235 {}
1236 
1237 template <typename MatrixType>
1238 inline void
1239 PreconditionIdentity::initialize(const MatrixType &matrix,
1241 {
1242  n_rows = matrix.m();
1243  n_columns = matrix.n();
1244 }
1245 
1246 
1247 template <class VectorType>
1248 inline void
1249 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1250 {
1251  dst = src;
1252 }
1253 
1254 
1255 
1256 template <class VectorType>
1257 inline void
1258 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1259 {
1260  dst = src;
1261 }
1262 
1263 template <class VectorType>
1264 inline void
1265 PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
1266 {
1267  dst += src;
1268 }
1269 
1270 
1271 
1272 template <class VectorType>
1273 inline void
1274 PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
1275 {
1276  dst += src;
1277 }
1278 
1281 {
1282  Assert(n_rows != 0, ExcNotInitialized());
1283  return n_rows;
1284 }
1285 
1288 {
1289  Assert(n_columns != 0, ExcNotInitialized());
1290  return n_columns;
1291 }
1292 
1293 //---------------------------------------------------------------------------
1294 
1296  const double relaxation)
1297  : relaxation(relaxation)
1298 {}
1299 
1300 
1302  : relaxation(0)
1303  , n_rows(0)
1304  , n_columns(0)
1305 {
1306  AdditionalData add_data;
1307  relaxation = add_data.relaxation;
1308 }
1309 
1310 
1311 
1312 inline void
1314  const PreconditionRichardson::AdditionalData &parameters)
1315 {
1316  relaxation = parameters.relaxation;
1317 }
1318 
1319 
1320 
1321 template <typename MatrixType>
1322 inline void
1324  const MatrixType & matrix,
1325  const PreconditionRichardson::AdditionalData &parameters)
1326 {
1327  relaxation = parameters.relaxation;
1328  n_rows = matrix.m();
1329  n_columns = matrix.n();
1330 }
1331 
1332 
1333 
1334 template <class VectorType>
1335 inline void
1336 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
1337 {
1338  static_assert(
1339  std::is_same<size_type, typename VectorType::size_type>::value,
1340  "PreconditionRichardson and VectorType must have the same size_type.");
1341 
1342  dst.equ(relaxation, src);
1343 }
1344 
1345 
1346 
1347 template <class VectorType>
1348 inline void
1349 PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
1350 {
1351  static_assert(
1352  std::is_same<size_type, typename VectorType::size_type>::value,
1353  "PreconditionRichardson and VectorType must have the same size_type.");
1354 
1355  dst.equ(relaxation, src);
1356 }
1357 
1358 template <class VectorType>
1359 inline void
1360 PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
1361 {
1362  static_assert(
1363  std::is_same<size_type, typename VectorType::size_type>::value,
1364  "PreconditionRichardson and VectorType must have the same size_type.");
1365 
1366  dst.add(relaxation, src);
1367 }
1368 
1369 
1370 
1371 template <class VectorType>
1372 inline void
1373 PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
1374 {
1375  static_assert(
1376  std::is_same<size_type, typename VectorType::size_type>::value,
1377  "PreconditionRichardson and VectorType must have the same size_type.");
1378 
1379  dst.add(relaxation, src);
1380 }
1381 
1384 {
1385  Assert(n_rows != 0, ExcNotInitialized());
1386  return n_rows;
1387 }
1388 
1391 {
1392  Assert(n_columns != 0, ExcNotInitialized());
1393  return n_columns;
1394 }
1395 
1396 //---------------------------------------------------------------------------
1397 
1398 template <typename MatrixType>
1399 inline void
1401  const AdditionalData &parameters)
1402 {
1403  A = &rA;
1404  relaxation = parameters.relaxation;
1405 }
1406 
1407 
1408 template <typename MatrixType>
1409 inline void
1411 {
1412  A = nullptr;
1413 }
1414 
1415 template <typename MatrixType>
1418 {
1419  Assert(A != nullptr, ExcNotInitialized());
1420  return A->m();
1421 }
1422 
1423 template <typename MatrixType>
1426 {
1427  Assert(A != nullptr, ExcNotInitialized());
1428  return A->n();
1429 }
1430 
1431 //---------------------------------------------------------------------------
1432 
1433 template <typename MatrixType>
1434 template <class VectorType>
1435 inline void
1436 PreconditionJacobi<MatrixType>::vmult(VectorType & dst,
1437  const VectorType &src) const
1438 {
1439  static_assert(
1440  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1441  typename VectorType::size_type>::value,
1442  "PreconditionJacobi and VectorType must have the same size_type.");
1443 
1444  Assert(this->A != nullptr, ExcNotInitialized());
1445  this->A->precondition_Jacobi(dst, src, this->relaxation);
1446 }
1447 
1448 
1449 
1450 template <typename MatrixType>
1451 template <class VectorType>
1452 inline void
1454  const VectorType &src) const
1455 {
1456  static_assert(
1457  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1458  typename VectorType::size_type>::value,
1459  "PreconditionJacobi and VectorType must have the same size_type.");
1460 
1461  Assert(this->A != nullptr, ExcNotInitialized());
1462  this->A->precondition_Jacobi(dst, src, this->relaxation);
1463 }
1464 
1465 
1466 
1467 template <typename MatrixType>
1468 template <class VectorType>
1469 inline void
1470 PreconditionJacobi<MatrixType>::step(VectorType & dst,
1471  const VectorType &src) const
1472 {
1473  static_assert(
1474  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1475  typename VectorType::size_type>::value,
1476  "PreconditionJacobi and VectorType must have the same size_type.");
1477 
1478  Assert(this->A != nullptr, ExcNotInitialized());
1479  this->A->Jacobi_step(dst, src, this->relaxation);
1480 }
1481 
1482 
1483 
1484 template <typename MatrixType>
1485 template <class VectorType>
1486 inline void
1487 PreconditionJacobi<MatrixType>::Tstep(VectorType & dst,
1488  const VectorType &src) const
1489 {
1490  static_assert(
1491  std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
1492  typename VectorType::size_type>::value,
1493  "PreconditionJacobi and VectorType must have the same size_type.");
1494 
1495  step(dst, src);
1496 }
1497 
1498 
1499 
1500 //---------------------------------------------------------------------------
1501 
1502 template <typename MatrixType>
1503 template <class VectorType>
1504 inline void
1505 PreconditionSOR<MatrixType>::vmult(VectorType &dst, const VectorType &src) const
1506 {
1507  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1508  typename VectorType::size_type>::value,
1509  "PreconditionSOR and VectorType must have the same size_type.");
1510 
1511  Assert(this->A != nullptr, ExcNotInitialized());
1512  this->A->precondition_SOR(dst, src, this->relaxation);
1513 }
1514 
1515 
1516 
1517 template <typename MatrixType>
1518 template <class VectorType>
1519 inline void
1520 PreconditionSOR<MatrixType>::Tvmult(VectorType & dst,
1521  const VectorType &src) const
1522 {
1523  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1524  typename VectorType::size_type>::value,
1525  "PreconditionSOR and VectorType must have the same size_type.");
1526 
1527  Assert(this->A != nullptr, ExcNotInitialized());
1528  this->A->precondition_TSOR(dst, src, this->relaxation);
1529 }
1530 
1531 
1532 
1533 template <typename MatrixType>
1534 template <class VectorType>
1535 inline void
1536 PreconditionSOR<MatrixType>::step(VectorType &dst, const VectorType &src) const
1537 {
1538  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1539  typename VectorType::size_type>::value,
1540  "PreconditionSOR and VectorType must have the same size_type.");
1541 
1542  Assert(this->A != nullptr, ExcNotInitialized());
1543  this->A->SOR_step(dst, src, this->relaxation);
1544 }
1545 
1546 
1547 
1548 template <typename MatrixType>
1549 template <class VectorType>
1550 inline void
1551 PreconditionSOR<MatrixType>::Tstep(VectorType &dst, const VectorType &src) const
1552 {
1553  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
1554  typename VectorType::size_type>::value,
1555  "PreconditionSOR and VectorType must have the same size_type.");
1556 
1557  Assert(this->A != nullptr, ExcNotInitialized());
1558  this->A->TSOR_step(dst, src, this->relaxation);
1559 }
1560 
1561 
1562 
1563 //---------------------------------------------------------------------------
1564 
1565 template <typename MatrixType>
1566 inline void
1568  const MatrixType & rA,
1569  const typename BaseClass::AdditionalData &parameters)
1570 {
1571  this->PreconditionRelaxation<MatrixType>::initialize(rA, parameters);
1572 
1573  // in case we have a SparseMatrix class, we can extract information about
1574  // the diagonal.
1576  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
1577  &*this->A);
1578 
1579  // calculate the positions first after the diagonal.
1580  if (mat != nullptr)
1581  {
1582  const size_type n = this->A->n();
1583  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1584  for (size_type row = 0; row < n; ++row)
1585  {
1586  // find the first element in this line which is on the right of the
1587  // diagonal. we need to precondition with the elements on the left
1588  // only. note: the first entry in each line denotes the diagonal
1589  // element, which we need not check.
1591  it = mat->begin(row) + 1;
1592  for (; it < mat->end(row); ++it)
1593  if (it->column() > row)
1594  break;
1595  pos_right_of_diagonal[row] = it - mat->begin();
1596  }
1597  }
1598 }
1599 
1600 
1601 template <typename MatrixType>
1602 template <class VectorType>
1603 inline void
1604 PreconditionSSOR<MatrixType>::vmult(VectorType & dst,
1605  const VectorType &src) const
1606 {
1607  static_assert(
1608  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1609  typename VectorType::size_type>::value,
1610  "PreconditionSSOR and VectorType must have the same size_type.");
1611 
1612  Assert(this->A != nullptr, ExcNotInitialized());
1613  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1614 }
1615 
1616 
1617 
1618 template <typename MatrixType>
1619 template <class VectorType>
1620 inline void
1621 PreconditionSSOR<MatrixType>::Tvmult(VectorType & dst,
1622  const VectorType &src) const
1623 {
1624  static_assert(
1625  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1626  typename VectorType::size_type>::value,
1627  "PreconditionSSOR and VectorType must have the same size_type.");
1628 
1629  Assert(this->A != nullptr, ExcNotInitialized());
1630  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
1631 }
1632 
1633 
1634 
1635 template <typename MatrixType>
1636 template <class VectorType>
1637 inline void
1638 PreconditionSSOR<MatrixType>::step(VectorType &dst, const VectorType &src) const
1639 {
1640  static_assert(
1641  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1642  typename VectorType::size_type>::value,
1643  "PreconditionSSOR and VectorType must have the same size_type.");
1644 
1645  Assert(this->A != nullptr, ExcNotInitialized());
1646  this->A->SSOR_step(dst, src, this->relaxation);
1647 }
1648 
1649 
1650 
1651 template <typename MatrixType>
1652 template <class VectorType>
1653 inline void
1654 PreconditionSSOR<MatrixType>::Tstep(VectorType & dst,
1655  const VectorType &src) const
1656 {
1657  static_assert(
1658  std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
1659  typename VectorType::size_type>::value,
1660  "PreconditionSSOR and VectorType must have the same size_type.");
1661 
1662  step(dst, src);
1663 }
1664 
1665 
1666 
1667 //---------------------------------------------------------------------------
1668 
1669 template <typename MatrixType>
1670 inline void
1672  const MatrixType & rA,
1673  const std::vector<size_type> & p,
1674  const std::vector<size_type> & ip,
1675  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1676 {
1677  permutation = &p;
1678  inverse_permutation = &ip;
1680 }
1681 
1682 
1683 template <typename MatrixType>
1684 inline void
1685 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
1686  const AdditionalData &additional_data)
1687 {
1688  initialize(A,
1689  additional_data.permutation,
1690  additional_data.inverse_permutation,
1691  additional_data.parameters);
1692 }
1693 
1694 
1695 template <typename MatrixType>
1696 template <typename VectorType>
1697 inline void
1698 PreconditionPSOR<MatrixType>::vmult(VectorType & dst,
1699  const VectorType &src) const
1700 {
1701  static_assert(
1702  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1703  typename VectorType::size_type>::value,
1704  "PreconditionPSOR and VectorType must have the same size_type.");
1705 
1706  Assert(this->A != nullptr, ExcNotInitialized());
1707  dst = src;
1708  this->A->PSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1709 }
1710 
1711 
1712 
1713 template <typename MatrixType>
1714 template <class VectorType>
1715 inline void
1716 PreconditionPSOR<MatrixType>::Tvmult(VectorType & dst,
1717  const VectorType &src) const
1718 {
1719  static_assert(
1720  std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
1721  typename VectorType::size_type>::value,
1722  "PreconditionPSOR and VectorType must have the same size_type.");
1723 
1724  Assert(this->A != nullptr, ExcNotInitialized());
1725  dst = src;
1726  this->A->TPSOR(dst, *permutation, *inverse_permutation, this->relaxation);
1727 }
1728 
1729 template <typename MatrixType>
1731  const std::vector<size_type> &permutation,
1732  const std::vector<size_type> &inverse_permutation,
1733  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
1734  : permutation(permutation)
1735  , inverse_permutation(inverse_permutation)
1736  , parameters(parameters)
1737 {}
1738 
1739 
1740 //---------------------------------------------------------------------------
1741 
1742 
1743 template <typename MatrixType, class VectorType>
1745  const MatrixType & M,
1746  const function_ptr method)
1747  : matrix(M)
1748  , precondition(method)
1749 {}
1750 
1751 
1752 
1753 template <typename MatrixType, class VectorType>
1754 void
1756  VectorType & dst,
1757  const VectorType &src) const
1758 {
1759  (matrix.*precondition)(dst, src);
1760 }
1761 
1762 //---------------------------------------------------------------------------
1763 
1764 template <typename MatrixType>
1766  const double relaxation)
1767  : relaxation(relaxation)
1768 {}
1769 
1770 
1771 
1772 //---------------------------------------------------------------------------
1773 
1774 namespace internal
1775 {
1776  namespace PreconditionChebyshevImplementation
1777  {
1778  // for deal.II vectors, perform updates for Chebyshev preconditioner all
1779  // at once to reduce memory transfer. Here, we select between general
1780  // vectors and deal.II vectors where we expand the loop over the (local)
1781  // size of the vector
1782 
1783  // generic part for non-deal.II vectors
1784  template <typename VectorType, typename PreconditionerType>
1785  inline void
1786  vector_updates(const VectorType & rhs,
1787  const PreconditionerType &preconditioner,
1788  const unsigned int iteration_index,
1789  const double factor1,
1790  const double factor2,
1791  VectorType & solution_old,
1792  VectorType & temp_vector1,
1793  VectorType & temp_vector2,
1794  VectorType & solution)
1795  {
1796  if (iteration_index == 0)
1797  {
1798  solution.equ(factor2, rhs);
1799  preconditioner.vmult(solution_old, solution);
1800  }
1801  else if (iteration_index == 1)
1802  {
1803  // compute t = P^{-1} * (b-A*x^{n})
1804  temp_vector1.sadd(-1.0, 1.0, rhs);
1805  preconditioner.vmult(solution_old, temp_vector1);
1806 
1807  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
1808  solution_old.sadd(factor2, 1 + factor1, solution);
1809  }
1810  else
1811  {
1812  // compute t = P^{-1} * (b-A*x^{n})
1813  temp_vector1.sadd(-1.0, 1.0, rhs);
1814  preconditioner.vmult(temp_vector2, temp_vector1);
1815 
1816  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
1817  solution_old.sadd(-factor1, factor2, temp_vector2);
1818  solution_old.add(1 + factor1, solution);
1819  }
1820 
1821  solution.swap(solution_old);
1822  }
1823 
1824  // worker routine for deal.II vectors. Because of vectorization, we need
1825  // to put the loop into an extra structure because the virtual function of
1826  // VectorUpdatesRange prevents the compiler from applying vectorization.
1827  template <typename Number>
1828  struct VectorUpdater
1829  {
1830  VectorUpdater(const Number * rhs,
1831  const Number * matrix_diagonal_inverse,
1832  const unsigned int iteration_index,
1833  const Number factor1,
1834  const Number factor2,
1835  Number * solution_old,
1836  Number * tmp_vector,
1837  Number * solution)
1838  : rhs(rhs)
1839  , matrix_diagonal_inverse(matrix_diagonal_inverse)
1840  , iteration_index(iteration_index)
1841  , factor1(factor1)
1842  , factor2(factor2)
1843  , solution_old(solution_old)
1844  , tmp_vector(tmp_vector)
1845  , solution(solution)
1846  {}
1847 
1848  void
1849  apply_to_subrange(const std::size_t begin, const std::size_t end) const
1850  {
1851  // To circumvent a bug in gcc
1852  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
1853  // copies of the variables factor1 and factor2 and do not check based on
1854  // factor1.
1855  const Number factor1 = this->factor1;
1856  const Number factor1_plus_1 = 1. + this->factor1;
1857  const Number factor2 = this->factor2;
1858  if (iteration_index == 0)
1859  {
1860  DEAL_II_OPENMP_SIMD_PRAGMA
1861  for (std::size_t i = begin; i < end; ++i)
1862  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
1863  }
1864  else if (iteration_index == 1)
1865  {
1866  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
1867  DEAL_II_OPENMP_SIMD_PRAGMA
1868  for (std::size_t i = begin; i < end; ++i)
1869  // for efficiency reason, write back to temp_vector that is
1870  // already read (avoid read-for-ownership)
1871  tmp_vector[i] =
1872  factor1_plus_1 * solution[i] +
1873  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1874  }
1875  else
1876  {
1877  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
1878  // + f_2 * P^{-1} * (b-A*x^{n})
1879  DEAL_II_OPENMP_SIMD_PRAGMA
1880  for (std::size_t i = begin; i < end; ++i)
1881  solution_old[i] =
1882  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
1883  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
1884  }
1885  }
1886 
1887  const Number * rhs;
1888  const Number * matrix_diagonal_inverse;
1889  const unsigned int iteration_index;
1890  const Number factor1;
1891  const Number factor2;
1892  mutable Number * solution_old;
1893  mutable Number * tmp_vector;
1894  mutable Number * solution;
1895  };
1896 
1897  template <typename Number>
1898  struct VectorUpdatesRange : public parallel::ParallelForInteger
1899  {
1900  VectorUpdatesRange(const VectorUpdater<Number> &updater,
1901  const std::size_t size)
1902  : updater(updater)
1903  {
1904  if (size < internal::VectorImplementation::minimum_parallel_grain_size)
1905  apply_to_subrange(0, size);
1906  else
1907  apply_parallel(
1908  0,
1909  size,
1910  internal::VectorImplementation::minimum_parallel_grain_size);
1911  }
1912 
1913  ~VectorUpdatesRange() override = default;
1914 
1915  virtual void
1916  apply_to_subrange(const std::size_t begin,
1917  const std::size_t end) const override
1918  {
1919  updater.apply_to_subrange(begin, end);
1920  }
1921 
1922  const VectorUpdater<Number> &updater;
1923  };
1924 
1925  // selection for diagonal matrix around deal.II vector
1926  template <typename Number>
1927  inline void
1928  vector_updates(const ::Vector<Number> & rhs,
1929  const DiagonalMatrix<::Vector<Number>> &jacobi,
1930  const unsigned int iteration_index,
1931  const double factor1,
1932  const double factor2,
1933  ::Vector<Number> &solution_old,
1934  ::Vector<Number> &temp_vector1,
1935  ::Vector<Number> &,
1936  ::Vector<Number> &solution)
1937  {
1938  VectorUpdater<Number> upd(rhs.begin(),
1939  jacobi.get_vector().begin(),
1940  iteration_index,
1941  factor1,
1942  factor2,
1943  solution_old.begin(),
1944  temp_vector1.begin(),
1945  solution.begin());
1946  VectorUpdatesRange<Number>(upd, rhs.size());
1947 
1948  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
1949  if (iteration_index == 0)
1950  {
1951  // nothing to do here because we can immediately write into the
1952  // solution vector without remembering any of the other vectors
1953  }
1954  else if (iteration_index == 1)
1955  {
1956  solution.swap(temp_vector1);
1957  solution_old.swap(temp_vector1);
1958  }
1959  else
1960  solution.swap(solution_old);
1961  }
1962 
1963  // selection for diagonal matrix around parallel deal.II vector
1964  template <typename Number>
1965  inline void
1966  vector_updates(
1968  const DiagonalMatrix<
1970  const unsigned int iteration_index,
1971  const double factor1,
1972  const double factor2,
1974  &solution_old,
1976  &temp_vector1,
1979  {
1980  VectorUpdater<Number> upd(rhs.begin(),
1981  jacobi.get_vector().begin(),
1982  iteration_index,
1983  factor1,
1984  factor2,
1985  solution_old.begin(),
1986  temp_vector1.begin(),
1987  solution.begin());
1988  VectorUpdatesRange<Number>(upd, rhs.local_size());
1989 
1990  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
1991  if (iteration_index == 0)
1992  {
1993  // nothing to do here because we can immediately write into the
1994  // solution vector without remembering any of the other vectors
1995  }
1996  else if (iteration_index == 1)
1997  {
1998  solution.swap(temp_vector1);
1999  solution_old.swap(temp_vector1);
2000  }
2001  else
2002  solution.swap(solution_old);
2003  }
2004 
2005  template <typename MatrixType, typename PreconditionerType>
2006  inline void
2007  initialize_preconditioner(
2008  const MatrixType & matrix,
2009  std::shared_ptr<PreconditionerType> &preconditioner)
2010  {
2011  (void)matrix;
2012  (void)preconditioner;
2013  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2014  }
2015 
2016  template <typename MatrixType, typename VectorType>
2017  inline void
2018  initialize_preconditioner(
2019  const MatrixType & matrix,
2020  std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
2021  {
2022  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2023  {
2024  if (preconditioner.get() == nullptr)
2025  preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2026 
2027  Assert(
2028  preconditioner->m() == 0,
2029  ExcMessage(
2030  "Preconditioner appears to be initialized but not sized correctly"));
2031 
2032  // This part only works in serial
2033  if (preconditioner->m() != matrix.m())
2034  {
2035  preconditioner->get_vector().reinit(matrix.m());
2036  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2037  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2038  }
2039  }
2040  }
2041 
2042  template <typename VectorType>
2043  void
2044  set_initial_guess(VectorType &vector)
2045  {
2046  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2047  if (vector.locally_owned_elements().is_element(0))
2048  vector(0) = 0.;
2049  }
2050 
2051  template <typename Number>
2052  void
2053  set_initial_guess(::Vector<Number> &vector)
2054  {
2055  // Choose a high-frequency mode consisting of numbers between 0 and 1
2056  // that is cheap to compute (cheaper than random numbers) but avoids
2057  // obviously re-occurring numbers in multi-component systems by choosing
2058  // a period of 11
2059  for (unsigned int i = 0; i < vector.size(); ++i)
2060  vector(i) = i % 11;
2061 
2062  const Number mean_value = vector.mean_value();
2063  vector.add(-mean_value);
2064  }
2065 
2066  template <typename Number>
2067  void
2068  set_initial_guess(
2070  &vector)
2071  {
2072  // Choose a high-frequency mode consisting of numbers between 0 and 1
2073  // that is cheap to compute (cheaper than random numbers) but avoids
2074  // obviously re-occurring numbers in multi-component systems by choosing
2075  // a period of 11.
2076  // Make initial guess robust with respect to number of processors
2077  // by operating on the global index.
2078  types::global_dof_index first_local_range = 0;
2079  if (!vector.locally_owned_elements().is_empty())
2080  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2081  for (unsigned int i = 0; i < vector.local_size(); ++i)
2082  vector.local_element(i) = (i + first_local_range) % 11;
2083 
2084  const Number mean_value = vector.mean_value();
2085  vector.add(-mean_value);
2086  }
2087 
2088 
2089 # ifdef DEAL_II_COMPILER_CUDA_AWARE
2090  template <typename Number>
2091  __global__ void
2092  set_initial_guess_kernel(const types::global_dof_index offset,
2093  const unsigned int local_size,
2094  Number * values)
2095 
2096  {
2097  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
2098  if (index < local_size)
2099  values[index] = (index + offset) % 11;
2100  }
2101 
2102  template <typename Number>
2103  void
2104  set_initial_guess(
2106  &vector)
2107  {
2108  // Choose a high-frequency mode consisting of numbers between 0 and 1
2109  // that is cheap to compute (cheaper than random numbers) but avoids
2110  // obviously re-occurring numbers in multi-component systems by choosing
2111  // a period of 11.
2112  // Make initial guess robust with respect to number of processors
2113  // by operating on the global index.
2114  types::global_dof_index first_local_range = 0;
2115  if (!vector.locally_owned_elements().is_empty())
2116  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2117 
2118  const auto n_local_elements = vector.local_size();
2119  const int n_blocks =
2120  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
2121  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
2122  first_local_range, n_local_elements, vector.get_values());
2123 
2124 # ifdef DEBUG
2125  // Check that the kernel was launched correctly
2126  AssertCuda(cudaGetLastError());
2127  // Check that there was no problem during the execution of the kernel
2128  AssertCuda(cudaDeviceSynchronize());
2129 # endif
2130 
2131  const Number mean_value = vector.mean_value();
2132  vector.add(-mean_value);
2133  }
2134 # endif // DEAL_II_COMPILER_CUDA_AWARE
2135 
2136  struct EigenvalueTracker
2137  {
2138  public:
2139  void
2140  slot(const std::vector<double> &eigenvalues)
2141  {
2142  values = eigenvalues;
2143  }
2144 
2145  std::vector<double> values;
2146  };
2147  } // namespace PreconditionChebyshevImplementation
2148 } // namespace internal
2149 
2150 
2151 
2152 template <typename MatrixType, class VectorType, typename PreconditionerType>
2154  AdditionalData::AdditionalData(const unsigned int degree,
2155  const double smoothing_range,
2156  const unsigned int eig_cg_n_iterations,
2157  const double eig_cg_residual,
2158  const double max_eigenvalue)
2159  : degree(degree)
2160  , smoothing_range(smoothing_range)
2161  , eig_cg_n_iterations(eig_cg_n_iterations)
2162  , eig_cg_residual(eig_cg_residual)
2163  , max_eigenvalue(max_eigenvalue)
2164 {}
2165 
2166 
2167 
2168 template <typename MatrixType, class VectorType, typename PreconditionerType>
2169 inline typename PreconditionChebyshev<MatrixType,
2170  VectorType,
2171  PreconditionerType>::AdditionalData &
2173  AdditionalData::operator=(const AdditionalData &other_data)
2174 {
2175  degree = other_data.degree;
2176  smoothing_range = other_data.smoothing_range;
2177  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
2178  eig_cg_residual = other_data.eig_cg_residual;
2179  max_eigenvalue = other_data.max_eigenvalue;
2180  preconditioner = other_data.preconditioner;
2181  constraints.copy_from(other_data.constraints);
2182 
2183  return *this;
2184 }
2185 
2186 
2187 
2188 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2191  : theta(1.)
2192  , delta(1.)
2194 {
2195  static_assert(
2196  std::is_same<size_type, typename VectorType::size_type>::value,
2197  "PreconditionChebyshev and VectorType must have the same size_type.");
2198 }
2199 
2200 
2201 
2202 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2203 inline void
2205  const MatrixType & matrix,
2206  const AdditionalData &additional_data)
2207 {
2208  matrix_ptr = &matrix;
2209  data = additional_data;
2210  Assert(data.degree > 0,
2211  ExcMessage("The degree of the Chebyshev method must be positive."));
2212  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2213  matrix, data.preconditioner);
2215 }
2216 
2217 
2218 
2219 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2220 inline void
2222 {
2224  theta = delta = 1.0;
2225  matrix_ptr = nullptr;
2226  {
2227  VectorType empty_vector;
2228  solution_old.reinit(empty_vector);
2229  temp_vector1.reinit(empty_vector);
2230  temp_vector2.reinit(empty_vector);
2231  }
2232  data.preconditioner.reset();
2233 }
2234 
2235 
2236 
2237 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2238 inline void
2240  estimate_eigenvalues(const VectorType &src) const
2241 {
2243  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2244 
2245  solution_old.reinit(src);
2246  temp_vector1.reinit(src, true);
2247 
2248  // calculate largest eigenvalue using a hand-tuned CG iteration on the
2249  // matrix weighted by its diagonal. we start with a vector that consists of
2250  // ones only, weighted by the length.
2251  double max_eigenvalue, min_eigenvalue;
2252  if (data.eig_cg_n_iterations > 0)
2253  {
2255  ExcMessage(
2256  "Need to set at least two iterations to find eigenvalues."));
2257 
2258  // set a very strict tolerance to force at least two iterations
2259  ReductionControl control(
2261  std::sqrt(
2262  std::numeric_limits<typename VectorType::value_type>::epsilon()),
2263  1e-10,
2264  false,
2265  false);
2266 
2267  internal::PreconditionChebyshevImplementation::EigenvalueTracker
2268  eigenvalue_tracker;
2269  SolverCG<VectorType> solver(control);
2270  solver.connect_eigenvalues_slot(
2271  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2272  eigenvalue_tracker.slot(eigenvalues);
2273  });
2274 
2275  // set an initial guess which is close to the constant vector but where
2276  // one entry is different to trigger high frequencies
2277  internal::PreconditionChebyshevImplementation::set_initial_guess(
2278  temp_vector1);
2279  data.constraints.set_zero(temp_vector1);
2280 
2281  try
2282  {
2283  solver.solve(*matrix_ptr,
2284  solution_old,
2285  temp_vector1,
2286  *data.preconditioner);
2287  }
2289  {}
2290 
2291  // read the eigenvalues from the attached eigenvalue tracker
2292  if (eigenvalue_tracker.values.empty())
2293  min_eigenvalue = max_eigenvalue = 1;
2294  else
2295  {
2296  min_eigenvalue = eigenvalue_tracker.values.front();
2297 
2298  // include a safety factor since the CG method will in general not
2299  // be converged
2300  max_eigenvalue = 1.2 * eigenvalue_tracker.values.back();
2301  }
2302  }
2303  else
2304  {
2305  max_eigenvalue = data.max_eigenvalue;
2306  min_eigenvalue = data.max_eigenvalue / data.smoothing_range;
2307  }
2308 
2309  const double alpha = (data.smoothing_range > 1. ?
2310  max_eigenvalue / data.smoothing_range :
2311  std::min(0.9 * max_eigenvalue, min_eigenvalue));
2312 
2313  // in case the user set the degree to invalid unsigned int, we have to
2314  // determine the number of necessary iterations from the Chebyshev error
2315  // estimate, given the target tolerance specified by smoothing_range. This
2316  // estimate is based on the error formula given in section 5.1 of
2317  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
2319  {
2320  const double actual_range = max_eigenvalue / alpha;
2321  const double sigma = (1. - std::sqrt(1. / actual_range)) /
2322  (1. + std::sqrt(1. / actual_range));
2323  const double eps = data.smoothing_range;
2324  const_cast<
2326  this)
2327  ->data.degree =
2328  1 + static_cast<unsigned int>(
2329  std::log(1. / eps + std::sqrt(1. / eps / eps - 1)) /
2330  std::log(1. / sigma));
2331  }
2332 
2333  const_cast<
2335  ->delta = (max_eigenvalue - alpha) * 0.5;
2336  const_cast<
2338  ->theta = (max_eigenvalue + alpha) * 0.5;
2339 
2340  // We do not need the second temporary vector in case we have a
2341  // DiagonalMatrix as preconditioner and use deal.II's own vectors
2342  using NumberType = typename VectorType::value_type;
2343  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
2344  false ||
2345  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
2346  ((std::is_same<VectorType,
2348  Vector<NumberType, MemorySpace::Host>>::value ==
2349  false) ||
2350  (std::is_same<VectorType,
2351  LinearAlgebra::distributed::
2352  Vector<NumberType, MemorySpace::CUDA>>::value ==
2353  false))))
2354  temp_vector2.reinit(src, true);
2355  else
2356  {
2357  VectorType empty_vector;
2358  temp_vector2.reinit(empty_vector);
2359  }
2360 
2361  const_cast<
2363  ->eigenvalues_are_initialized = true;
2364 }
2365 
2366 
2367 
2368 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2369 inline void
2371  VectorType & solution,
2372  const VectorType &rhs) const
2373 {
2374  std::lock_guard<std::mutex> lock(mutex);
2375  if (eigenvalues_are_initialized == false)
2376  estimate_eigenvalues(rhs);
2377 
2378  internal::PreconditionChebyshevImplementation::vector_updates(
2379  rhs,
2381  0,
2382  0.,
2383  1. / theta,
2384  solution_old,
2385  temp_vector1,
2386  temp_vector2,
2387  solution);
2388 
2389  // if delta is zero, we do not need to iterate because the updates will be
2390  // zero
2391  if (data.degree < 2 || std::abs(delta) < 1e-40)
2392  return;
2393 
2394  double rhok = delta / theta, sigma = theta / delta;
2395  for (unsigned int k = 0; k < data.degree - 1; ++k)
2396  {
2397  matrix_ptr->vmult(temp_vector1, solution);
2398  const double rhokp = 1. / (2. * sigma - rhok);
2399  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2400  rhok = rhokp;
2401  internal::PreconditionChebyshevImplementation::vector_updates(
2402  rhs,
2404  k + 1,
2405  factor1,
2406  factor2,
2407  solution_old,
2408  temp_vector1,
2409  temp_vector2,
2410  solution);
2411  }
2412 }
2413 
2414 
2415 
2416 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2417 inline void
2419  VectorType & solution,
2420  const VectorType &rhs) const
2421 {
2422  std::lock_guard<std::mutex> lock(mutex);
2423  if (eigenvalues_are_initialized == false)
2424  estimate_eigenvalues(rhs);
2425 
2426  internal::PreconditionChebyshevImplementation::vector_updates(
2427  rhs,
2429  0,
2430  0.,
2431  1. / theta,
2432  solution_old,
2433  temp_vector1,
2434  temp_vector2,
2435  solution);
2436 
2437  if (data.degree < 2 || std::abs(delta) < 1e-40)
2438  return;
2439 
2440  double rhok = delta / theta, sigma = theta / delta;
2441  for (unsigned int k = 0; k < data.degree - 1; ++k)
2442  {
2443  matrix_ptr->Tvmult(temp_vector1, solution);
2444  const double rhokp = 1. / (2. * sigma - rhok);
2445  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2446  rhok = rhokp;
2447  internal::PreconditionChebyshevImplementation::vector_updates(
2448  rhs,
2450  k + 1,
2451  factor1,
2452  factor2,
2453  solution_old,
2454  temp_vector1,
2455  temp_vector2,
2456  solution);
2457  }
2458 }
2459 
2460 
2461 
2462 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2463 inline void
2465  VectorType & solution,
2466  const VectorType &rhs) const
2467 {
2468  std::lock_guard<std::mutex> lock(mutex);
2469  if (eigenvalues_are_initialized == false)
2470  estimate_eigenvalues(rhs);
2471 
2472  matrix_ptr->vmult(temp_vector1, solution);
2473  internal::PreconditionChebyshevImplementation::vector_updates(
2474  rhs,
2476  1,
2477  0.,
2478  1. / theta,
2479  solution_old,
2480  temp_vector1,
2481  temp_vector2,
2482  solution);
2483 
2484  if (data.degree < 2 || std::abs(delta) < 1e-40)
2485  return;
2486 
2487  double rhok = delta / theta, sigma = theta / delta;
2488  for (unsigned int k = 0; k < data.degree - 1; ++k)
2489  {
2490  matrix_ptr->vmult(temp_vector1, solution);
2491  const double rhokp = 1. / (2. * sigma - rhok);
2492  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2493  rhok = rhokp;
2494  internal::PreconditionChebyshevImplementation::vector_updates(
2495  rhs,
2497  k + 2,
2498  factor1,
2499  factor2,
2500  solution_old,
2501  temp_vector1,
2502  temp_vector2,
2503  solution);
2504  }
2505 }
2506 
2507 
2508 
2509 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2510 inline void
2512  VectorType & solution,
2513  const VectorType &rhs) const
2514 {
2515  std::lock_guard<std::mutex> lock(mutex);
2516  if (eigenvalues_are_initialized == false)
2517  estimate_eigenvalues(rhs);
2518 
2519  matrix_ptr->Tvmult(temp_vector1, solution);
2520  internal::PreconditionChebyshevImplementation::vector_updates(
2521  rhs,
2523  1,
2524  0.,
2525  1. / theta,
2526  solution_old,
2527  temp_vector1,
2528  temp_vector2,
2529  solution);
2530 
2531  if (data.degree < 2 || std::abs(delta) < 1e-40)
2532  return;
2533 
2534  double rhok = delta / theta, sigma = theta / delta;
2535  for (unsigned int k = 0; k < data.degree - 1; ++k)
2536  {
2537  matrix_ptr->Tvmult(temp_vector1, solution);
2538  const double rhokp = 1. / (2. * sigma - rhok);
2539  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
2540  rhok = rhokp;
2541  internal::PreconditionChebyshevImplementation::vector_updates(
2542  rhs,
2544  k + 2,
2545  factor1,
2546  factor2,
2547  solution_old,
2548  temp_vector1,
2549  temp_vector2,
2550  solution);
2551  }
2552 }
2553 
2554 
2555 
2556 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2557 inline typename PreconditionChebyshev<MatrixType,
2558  VectorType,
2559  PreconditionerType>::size_type
2561 {
2562  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2563  return matrix_ptr->m();
2564 }
2565 
2566 
2567 
2568 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2569 inline typename PreconditionChebyshev<MatrixType,
2570  VectorType,
2571  PreconditionerType>::size_type
2573 {
2574  Assert(matrix_ptr != nullptr, ExcNotInitialized());
2575  return matrix_ptr->n();
2576 }
2577 
2578 #endif // DOXYGEN
2579 
2580 DEAL_II_NAMESPACE_CLOSE
2581 
2582 #endif
void Tvmult(VectorType &, const VectorType &) const
static const unsigned int invalid_unsigned_int
Definition: types.h:187
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
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1847
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
void reinit(const size_type size, const bool omit_zeroing_entries=false)
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:1519
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:1407
typename PreconditionRelaxation< MatrixType >::AdditionalData AdditionalData
Definition: precondition.h:672
#define AssertCuda(error_code)
Definition: exceptions.h:1722
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 estimate_eigenvalues(const VectorType &src) const
void Tstep(VectorType &x, const VectorType &rhs) const
void vmult(VectorType &, const VectorType &) const
void Tstep(VectorType &x, const VectorType &rhs) const
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:89
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())
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
bool is_empty() const
Definition: index_set.h:1791
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())