Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
precondition.h
1 // ---------------------------------------------------------------------
2 // @f$Id: precondition.h 31481 2013-10-29 12:20:59Z kronbichler @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__precondition_h
18 #define __deal2__precondition_h
19 
20 // This file contains simple preconditioners.
21 
22 #include <deal.II/base/config.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/base/utilities.h>
25 #include <deal.II/base/parallel.h>
26 #include <deal.II/base/template_constraints.h>
27 #include <deal.II/lac/tridiagonal_matrix.h>
28 #include <deal.II/lac/solver_cg.h>
29 #include <deal.II/lac/vector_memory.h>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 // forward declarations
34 
35 template <typename number> class Vector;
36 template <typename number> class SparseMatrix;
37 namespace parallel
38 {
39  namespace distributed
40  {
41  template <typename number> class Vector;
42  }
43 }
44 
45 
46 
74 {
75 public:
85  {
90  };
91 
98  template <class MATRIX>
99  void initialize (const MATRIX &matrix,
100  const AdditionalData &additional_data = AdditionalData());
101 
105  template<class VECTOR>
106  void vmult (VECTOR &, const VECTOR &) const;
107 
115  template<class VECTOR>
116  void Tvmult (VECTOR &, const VECTOR &) const;
117 
121  template<class VECTOR>
122  void vmult_add (VECTOR &, const VECTOR &) const;
123 
131  template<class VECTOR>
132  void Tvmult_add (VECTOR &, const VECTOR &) const;
133 
142  void clear () {}
143 };
144 
145 
146 
160 {
161 public:
167  {
168  public:
175  AdditionalData (const double relaxation = 1.);
176 
180  double relaxation;
181  };
182 
189 
193  void initialize (const AdditionalData &parameters);
194 
203  template <class MATRIX>
204  void initialize (const MATRIX &,
205  const AdditionalData &parameters);
206 
210  template<class VECTOR>
211  void vmult (VECTOR &, const VECTOR &) const;
212 
220  template<class VECTOR>
221  void Tvmult (VECTOR &, const VECTOR &) const;
225  template<class VECTOR>
226  void vmult_add (VECTOR &, const VECTOR &) const;
227 
235  template<class VECTOR>
236  void Tvmult_add (VECTOR &, const VECTOR &) const;
237 
246  void clear () {}
247 
248 private:
253  double relaxation;
254 };
255 
256 
257 
299 template<class MATRIX = SparseMatrix<double>, class VECTOR = Vector<double> >
301 {
302 public:
307  typedef void ( MATRIX::* function_ptr)(VECTOR &, const VECTOR &) const;
308 
318  PreconditionUseMatrix(const MATRIX &M,
319  const function_ptr method);
320 
327  void vmult (VECTOR &dst,
328  const VECTOR &src) const;
329 
330 private:
334  const MATRIX &matrix;
335 
341 };
342 
343 
344 
352 template<class MATRIX = SparseMatrix<double> >
354 {
355 public:
360  {
361  public:
365  AdditionalData (const double relaxation = 1.);
366 
370  double relaxation;
371  };
372 
383  void initialize (const MATRIX &A,
384  const AdditionalData &parameters = AdditionalData());
385 
390  void clear();
391 
392 protected:
397 
401  double relaxation;
402 };
403 
404 
405 
431 template <class MATRIX = SparseMatrix<double> >
433 {
434 public:
438  template<class VECTOR>
439  void vmult (VECTOR &, const VECTOR &) const;
440 
448  template<class VECTOR>
449  void Tvmult (VECTOR &, const VECTOR &) const;
450 
456  template<class VECTOR>
457  void step (VECTOR &x, const VECTOR &rhs) const;
458 
464  template<class VECTOR>
465  void Tstep (VECTOR &x, const VECTOR &rhs) const;
466 };
467 
468 
515 template <class MATRIX = SparseMatrix<double> >
517 {
518 public:
522  template<class VECTOR>
523  void vmult (VECTOR &, const VECTOR &) const;
524 
529  template<class VECTOR>
530  void Tvmult (VECTOR &, const VECTOR &) const;
531 
537  template<class VECTOR>
538  void step (VECTOR &x, const VECTOR &rhs) const;
539 
545  template<class VECTOR>
546  void Tstep (VECTOR &x, const VECTOR &rhs) const;
547 };
548 
549 
550 
576 template <class MATRIX = SparseMatrix<double> >
578 {
579 public:
584 
589 
590 
601  void initialize (const MATRIX &A,
602  const typename BaseClass::AdditionalData &parameters = typename BaseClass::AdditionalData());
603 
607  template<class VECTOR>
608  void vmult (VECTOR &, const VECTOR &) const;
609 
617  template<class VECTOR>
618  void Tvmult (VECTOR &, const VECTOR &) const;
619 
620 
626  template<class VECTOR>
627  void step (VECTOR &x, const VECTOR &rhs) const;
628 
634  template<class VECTOR>
635  void Tstep (VECTOR &x, const VECTOR &rhs) const;
636 
637 private:
643  std::vector<std::size_t> pos_right_of_diagonal;
644 };
645 
646 
678 template <class MATRIX = SparseMatrix<double> >
680 {
681 public:
686 
706  void initialize (const MATRIX &A,
707  const std::vector<size_type> &permutation,
708  const std::vector<size_type> &inverse_permutation,
710  parameters = typename PreconditionRelaxation<MATRIX>::AdditionalData());
711 
715  template<class VECTOR>
716  void vmult (VECTOR &, const VECTOR &) const;
717 
722  template<class VECTOR>
723  void Tvmult (VECTOR &, const VECTOR &) const;
724 private:
728  const std::vector<size_type> *permutation;
733  const std::vector<size_type> *inverse_permutation;
734 };
735 
736 
737 
789 template<class SOLVER, class MATRIX = SparseMatrix<double>, class PRECONDITION = PreconditionIdentity>
791 {
792 public:
798 
805  void initialize (SOLVER &,
806  const MATRIX &,
807  const PRECONDITION &);
808 
812  template<class VECTOR>
813  void vmult (VECTOR &, const VECTOR &) const;
814 
815 private:
820 
825 
831 
832 
833 
847 template<class MATRIX, class PRECOND, class VECTOR>
849 {
850 public:
857  PreconditionedMatrix (const MATRIX &A,
858  const PRECOND &P,
860 
865  void vmult (VECTOR &dst, const VECTOR &src) const;
866 
871  void Tvmult (VECTOR &dst, const VECTOR &src) const;
872 
876  double residual (VECTOR &dst, const VECTOR &src, const VECTOR &rhs) const;
877 
878 private:
882  const MATRIX &A;
886  const PRECOND &P;
892 
893 
894 
909 template <class MATRIX=SparseMatrix<double>, class VECTOR=Vector<double> >
911 {
912 public:
917 
923  {
927  AdditionalData (const unsigned int degree = 0,
928  const double smoothing_range = 0.,
929  const bool nonzero_starting = false,
930  const unsigned int eig_cg_n_iterations = 8,
931  const double eig_cg_residual = 1e-2,
932  const double max_eigenvalue = 1);
933 
940  unsigned int degree;
941 
954 
965 
971  unsigned int eig_cg_n_iterations;
972 
978 
985 
990  };
991 
993 
1005  void initialize (const MATRIX &matrix,
1006  const AdditionalData &additional_data = AdditionalData());
1007 
1012  void vmult (VECTOR &dst,
1013  const VECTOR &src) const;
1014 
1019  void Tvmult (VECTOR &dst,
1020  const VECTOR &src) const;
1021 
1025  void clear ();
1026 
1027 private:
1028 
1033 
1037  mutable VECTOR update1;
1038 
1042  mutable VECTOR update2;
1043 
1048 
1052  double theta;
1053 
1058  double delta;
1059 
1064 };
1065 
1066 
1067 
1069 /* ---------------------------------- Inline functions ------------------- */
1070 
1071 #ifndef DOXYGEN
1072 
1073 template <class MATRIX>
1074 inline void
1076  const MATRIX &,
1078 {}
1079 
1080 
1081 template<class VECTOR>
1082 inline void
1083 PreconditionIdentity::vmult (VECTOR &dst, const VECTOR &src) const
1084 {
1085  dst = src;
1086 }
1087 
1088 
1089 
1090 template<class VECTOR>
1091 inline void
1092 PreconditionIdentity::Tvmult (VECTOR &dst, const VECTOR &src) const
1093 {
1094  dst = src;
1095 }
1096 
1097 template<class VECTOR>
1098 inline void
1099 PreconditionIdentity::vmult_add (VECTOR &dst, const VECTOR &src) const
1100 {
1101  dst.add(src);
1102 }
1103 
1104 
1105 
1106 template<class VECTOR>
1107 inline void
1108 PreconditionIdentity::Tvmult_add (VECTOR &dst, const VECTOR &src) const
1109 {
1110  dst.add(src);
1111 }
1112 
1113 //---------------------------------------------------------------------------
1114 
1115 inline
1117  const double relaxation)
1118  :
1119  relaxation(relaxation)
1120 {}
1121 
1122 
1123 inline
1125  :
1126  relaxation(0)
1127 {
1128  AdditionalData add_data;
1129  relaxation=add_data.relaxation;
1130 }
1131 
1132 
1133 
1134 inline void
1136  const PreconditionRichardson::AdditionalData &parameters)
1137 {
1138  relaxation = parameters.relaxation;
1139 }
1140 
1141 
1142 
1143 template <class MATRIX>
1144 inline void
1146  const MATRIX &,
1147  const PreconditionRichardson::AdditionalData &parameters)
1148 {
1149  relaxation = parameters.relaxation;
1150 }
1151 
1152 
1153 
1154 template<class VECTOR>
1155 inline void
1156 PreconditionRichardson::vmult (VECTOR &dst, const VECTOR &src) const
1157 {
1158  dst.equ(relaxation,src);
1159 }
1160 
1161 
1162 
1163 template<class VECTOR>
1164 inline void
1165 PreconditionRichardson::Tvmult (VECTOR &dst, const VECTOR &src) const
1166 {
1167  dst.equ(relaxation,src);
1168 }
1169 
1170 template<class VECTOR>
1171 inline void
1172 PreconditionRichardson::vmult_add (VECTOR &dst, const VECTOR &src) const
1173 {
1174  dst.add(relaxation,src);
1175 }
1176 
1177 
1178 
1179 template<class VECTOR>
1180 inline void
1181 PreconditionRichardson::Tvmult_add (VECTOR &dst, const VECTOR &src) const
1182 {
1183  dst.add(relaxation,src);
1184 }
1185 
1186 //---------------------------------------------------------------------------
1187 
1188 template <class MATRIX>
1189 inline void
1191  const AdditionalData &parameters)
1192 {
1193  A = &rA;
1194  relaxation = parameters.relaxation;
1195 }
1196 
1197 
1198 template <class MATRIX>
1199 inline void
1201 {
1202  A = 0;
1203 }
1204 
1205 
1206 //---------------------------------------------------------------------------
1207 
1208 template <class MATRIX>
1209 template<class VECTOR>
1210 inline void
1211 PreconditionJacobi<MATRIX>::vmult (VECTOR &dst, const VECTOR &src) const
1212 {
1213  Assert (this->A!=0, ExcNotInitialized());
1214  this->A->precondition_Jacobi (dst, src, this->relaxation);
1215 }
1216 
1217 
1218 
1219 template <class MATRIX>
1220 template<class VECTOR>
1221 inline void
1222 PreconditionJacobi<MATRIX>::Tvmult (VECTOR &dst, const VECTOR &src) const
1223 {
1224  Assert (this->A!=0, ExcNotInitialized());
1225  this->A->precondition_Jacobi (dst, src, this->relaxation);
1226 }
1227 
1228 
1229 
1230 template <class MATRIX>
1231 template<class VECTOR>
1232 inline void
1233 PreconditionJacobi<MATRIX>::step (VECTOR &dst, const VECTOR &src) const
1234 {
1235  Assert (this->A!=0, ExcNotInitialized());
1236  this->A->Jacobi_step (dst, src, this->relaxation);
1237 }
1238 
1239 
1240 
1241 template <class MATRIX>
1242 template<class VECTOR>
1243 inline void
1244 PreconditionJacobi<MATRIX>::Tstep (VECTOR &dst, const VECTOR &src) const
1245 {
1246  step (dst, src);
1247 }
1248 
1249 
1250 
1251 //---------------------------------------------------------------------------
1252 
1253 template <class MATRIX>
1254 template<class VECTOR>
1255 inline void
1256 PreconditionSOR<MATRIX>::vmult (VECTOR &dst, const VECTOR &src) const
1257 {
1258  Assert (this->A!=0, ExcNotInitialized());
1259  this->A->precondition_SOR (dst, src, this->relaxation);
1260 }
1261 
1262 
1263 
1264 template <class MATRIX>
1265 template<class VECTOR>
1266 inline void
1267 PreconditionSOR<MATRIX>::Tvmult (VECTOR &dst, const VECTOR &src) const
1268 {
1269  Assert (this->A!=0, ExcNotInitialized());
1270  this->A->precondition_TSOR (dst, src, this->relaxation);
1271 }
1272 
1273 
1274 
1275 template <class MATRIX>
1276 template<class VECTOR>
1277 inline void
1278 PreconditionSOR<MATRIX>::step (VECTOR &dst, const VECTOR &src) const
1279 {
1280  Assert (this->A!=0, ExcNotInitialized());
1281  this->A->SOR_step (dst, src, this->relaxation);
1282 }
1283 
1284 
1285 
1286 template <class MATRIX>
1287 template<class VECTOR>
1288 inline void
1289 PreconditionSOR<MATRIX>::Tstep (VECTOR &dst, const VECTOR &src) const
1290 {
1291  Assert (this->A!=0, ExcNotInitialized());
1292  this->A->TSOR_step (dst, src, this->relaxation);
1293 }
1294 
1295 
1296 
1297 //---------------------------------------------------------------------------
1298 
1299 template <class MATRIX>
1300 inline void
1302  const typename BaseClass::AdditionalData &parameters)
1303 {
1304  this->PreconditionRelaxation<MATRIX>::initialize (rA, parameters);
1305 
1306  // in case we have a SparseMatrix class, we can extract information about
1307  // the diagonal.
1309  dynamic_cast<const SparseMatrix<typename MATRIX::value_type> *>(&*this->A);
1310 
1311  // calculate the positions first after the diagonal.
1312  if (mat != 0)
1313  {
1314  const size_type n = this->A->n();
1315  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
1316  for (size_type row=0; row<n; ++row)
1317  {
1318  // find the first element in this line which is on the right of the
1319  // diagonal. we need to precondition with the elements on the left
1320  // only. note: the first entry in each line denotes the diagonal
1321  // element, which we need not check.
1323  it = mat->begin(row)+1;
1324  for ( ; it < mat->end(row); ++it)
1325  if (it->column() > row)
1326  break;
1327  pos_right_of_diagonal[row] = it - mat->begin();
1328  }
1329  }
1330 }
1331 
1332 
1333 template <class MATRIX>
1334 template<class VECTOR>
1335 inline void
1336 PreconditionSSOR<MATRIX>::vmult (VECTOR &dst, const VECTOR &src) const
1337 {
1338  Assert (this->A!=0, ExcNotInitialized());
1339  this->A->precondition_SSOR (dst, src, this->relaxation, pos_right_of_diagonal);
1340 }
1341 
1342 
1343 
1344 template <class MATRIX>
1345 template<class VECTOR>
1346 inline void
1347 PreconditionSSOR<MATRIX>::Tvmult (VECTOR &dst, const VECTOR &src) const
1348 {
1349  Assert (this->A!=0, ExcNotInitialized());
1350  this->A->precondition_SSOR (dst, src, this->relaxation, pos_right_of_diagonal);
1351 }
1352 
1353 
1354 
1355 template <class MATRIX>
1356 template<class VECTOR>
1357 inline void
1358 PreconditionSSOR<MATRIX>::step (VECTOR &dst, const VECTOR &src) const
1359 {
1360  Assert (this->A!=0, ExcNotInitialized());
1361  this->A->SSOR_step (dst, src, this->relaxation);
1362 }
1363 
1364 
1365 
1366 template <class MATRIX>
1367 template<class VECTOR>
1368 inline void
1369 PreconditionSSOR<MATRIX>::Tstep (VECTOR &dst, const VECTOR &src) const
1370 {
1371  step (dst, src);
1372 }
1373 
1374 
1375 
1376 //---------------------------------------------------------------------------
1377 
1378 template <class MATRIX>
1379 inline void
1381  const MATRIX &rA,
1382  const std::vector<size_type> &p,
1383  const std::vector<size_type> &ip,
1384  const typename PreconditionRelaxation<MATRIX>::AdditionalData &parameters)
1385 {
1386  permutation = &p;
1387  inverse_permutation = &ip;
1389 }
1390 
1391 
1392 template <class MATRIX>
1393 template<class VECTOR>
1394 inline void
1395 PreconditionPSOR<MATRIX>::vmult (VECTOR &dst, const VECTOR &src) const
1396 {
1397  Assert (this->A!=0, ExcNotInitialized());
1398  dst = src;
1399  this->A->PSOR (dst, *permutation, *inverse_permutation, this->relaxation);
1400 }
1401 
1402 
1403 
1404 template <class MATRIX>
1405 template<class VECTOR>
1406 inline void
1407 PreconditionPSOR<MATRIX>::Tvmult (VECTOR &dst, const VECTOR &src) const
1408 {
1409  Assert (this->A!=0, ExcNotInitialized());
1410  dst = src;
1411  this->A->TPSOR (dst, *permutation, *inverse_permutation, this->relaxation);
1412 }
1413 
1414 
1415 //---------------------------------------------------------------------------
1416 
1417 
1418 template<class MATRIX, class VECTOR>
1420  const function_ptr method)
1421  :
1422  matrix(M), precondition(method)
1423 {}
1424 
1425 
1426 
1427 template<class MATRIX, class VECTOR>
1428 void
1430  const VECTOR &src) const
1431 {
1432  (matrix.*precondition)(dst, src);
1433 }
1434 
1435 //---------------------------------------------------------------------------
1436 
1437 template<class MATRIX>
1438 inline
1440 AdditionalData (const double relaxation)
1441  :
1442  relaxation (relaxation)
1443 {}
1444 
1445 
1446 
1448 
1449 template<class SOLVER, class MATRIX, class PRECONDITION>
1452  :
1453  solver(0), matrix(0), precondition(0)
1454 {}
1455 
1456 
1457 template<class SOLVER, class MATRIX, class PRECONDITION>
1458 void
1460 ::initialize (SOLVER &s,
1461  const MATRIX &m,
1462  const PRECONDITION &p)
1463 {
1464  solver = &s;
1465  matrix = &m;
1466  precondition = &p;
1467 }
1468 
1469 
1470 template<class SOLVER, class MATRIX, class PRECONDITION>
1471 template<class VECTOR>
1472 void
1474  const VECTOR &src) const
1475 {
1476  Assert (solver !=0 && matrix != 0 && precondition != 0,
1477  ExcNotInitialized());
1478 
1479  solver->solve(*matrix, dst, src, *precondition);
1480 }
1481 
1483 
1484 
1485 template<class MATRIX, class PRECOND, class VECTOR>
1486 inline
1488 ::PreconditionedMatrix (const MATRIX &A,
1489  const PRECOND &P,
1490  VectorMemory<VECTOR> &mem):
1491  A(A), P(P), mem(mem)
1492 {}
1493 
1494 
1495 template<class MATRIX, class PRECOND, class VECTOR>
1496 inline void
1498 ::vmult (VECTOR &dst,
1499  const VECTOR &src) const
1500 {
1501  VECTOR *h = mem.alloc();
1502  h->reinit(src);
1503  A.vmult(*h, src);
1504  P.vmult(dst, *h);
1505  mem.free(h);
1506 }
1507 
1508 
1509 
1510 template<class MATRIX, class PRECOND, class VECTOR>
1511 inline void
1513 ::Tvmult (VECTOR &dst,
1514  const VECTOR &src) const
1515 {
1516  VECTOR *h = mem.alloc();
1517  h->reinit(src);
1518  A.Tvmult(*h, src);
1519  P.Tvmult(dst, *h);
1520  mem.free(h);
1521 }
1522 
1523 
1524 
1525 template<class MATRIX, class PRECOND, class VECTOR>
1526 inline double
1528 ::residual (VECTOR &dst,
1529  const VECTOR &src,
1530  const VECTOR &rhs) const
1531 {
1532  VECTOR *h = mem.alloc();
1533  h->reinit(src);
1534  A.vmult(*h, src);
1535  P.vmult(dst, *h);
1536  mem.free(h);
1537  dst.sadd(-1.,1.,rhs);
1538  return dst.l2_norm ();
1539 }
1540 
1541 
1542 
1543 //---------------------------------------------------------------------------
1544 
1545 namespace internal
1546 {
1547  namespace PreconditionChebyshev
1548  {
1549  // for deal.II vectors, perform updates for Chebyshev preconditioner all
1550  // at once to reduce memory transfer. Here, we select between general
1551  // vectors and deal.II vectors where we expand the loop over the (local)
1552  // size of the vector
1553 
1554  // generic part for non-deal.II vectors
1555  template <typename VECTOR>
1556  inline
1557  void
1558  vector_updates (const VECTOR &src,
1559  const VECTOR &matrix_diagonal_inverse,
1560  const bool start_zero,
1561  const double factor1,
1562  const double factor2,
1563  VECTOR &update1,
1564  VECTOR &update2,
1565  VECTOR &dst)
1566  {
1567  if (start_zero)
1568  {
1569  dst.equ (factor2, src);
1570  dst.scale (matrix_diagonal_inverse);
1571  update1.equ(-1.,dst);
1572  }
1573  else
1574  {
1575  update2 -= src;
1576  update2.scale (matrix_diagonal_inverse);
1577  if (factor1 == 0.)
1578  update1.equ(factor2, update2);
1579  else
1580  update1.sadd(factor1, factor2, update2);
1581  dst -= update1;
1582  }
1583  }
1584 
1585  // worker loop for deal.II vectors
1586  template <typename Number>
1587  struct VectorUpdatesRange : public parallel::ParallelForInteger
1588  {
1590 
1591  VectorUpdatesRange (const size_t size,
1592  const Number *src,
1593  const Number *matrix_diagonal_inverse,
1594  const bool start_zero,
1595  const Number factor1,
1596  const Number factor2,
1597  Number *update1,
1598  Number *update2,
1599  Number *dst)
1600  :
1601  src (src),
1602  matrix_diagonal_inverse (matrix_diagonal_inverse),
1603  start_zero (start_zero),
1604  factor1 (factor1),
1605  factor2 (factor2),
1606  update1 (update1),
1607  update2 (update2),
1608  dst (dst)
1609  {
1610  if (size < internal::Vector::minimum_parallel_grain_size)
1611  apply_to_subrange (0, size);
1612  else
1613  apply_parallel (0, size,
1614  internal::Vector::minimum_parallel_grain_size);
1615  }
1616 
1617  ~VectorUpdatesRange()
1618  {}
1619 
1620  virtual void
1621  apply_to_subrange (const size_t begin,
1622  const size_t end) const
1623  {
1624  if (factor1 == Number())
1625  {
1626  if (start_zero)
1627  for (size_type i=begin; i<end; ++i)
1628  {
1629  dst[i] = factor2 * src[i] * matrix_diagonal_inverse[i];
1630  update1[i] = -dst[i];
1631  }
1632  else
1633  for (size_type i=begin; i<end; ++i)
1634  {
1635  update1[i] = ((update2[i]-src[i]) *
1636  factor2*matrix_diagonal_inverse[i]);
1637  dst[i] -= update1[i];
1638  }
1639  }
1640  else
1641  for (size_type i=begin; i<end; ++i)
1642  {
1643  const Number update2i = ((update2[i] - src[i]) *
1644  matrix_diagonal_inverse[i]);
1645  update1[i] = factor1 * update1[i] + factor2 * update2i;
1646  dst[i] -= update1[i];
1647  }
1648  }
1649 
1650  const Number *src;
1651  const Number *matrix_diagonal_inverse;
1652  const bool start_zero;
1653  const Number factor1;
1654  const Number factor2;
1655  mutable Number *update1;
1656  mutable Number *update2;
1657  mutable Number *dst;
1658  };
1659 
1660  // selection for deal.II vector
1661  template <typename Number>
1662  inline
1663  void
1664  vector_updates (const ::Vector<Number> &src,
1665  const ::Vector<Number> &matrix_diagonal_inverse,
1666  const bool start_zero,
1667  const double factor1,
1668  const double factor2,
1669  ::Vector<Number> &update1,
1670  ::Vector<Number> &update2,
1671  ::Vector<Number> &dst)
1672  {
1673  VectorUpdatesRange<Number>(src.size(), src.begin(),
1674  matrix_diagonal_inverse.begin(),
1675  start_zero, factor1, factor2,
1676  update1.begin(), update2.begin(), dst.begin());
1677  }
1678 
1679  // selection for parallel deal.II vector
1680  template <typename Number>
1681  inline
1682  void
1683  vector_updates (const parallel::distributed::Vector<Number> &src,
1684  const parallel::distributed::Vector<Number> &matrix_diagonal_inverse,
1685  const bool start_zero,
1686  const double factor1,
1687  const double factor2,
1691  {
1692  VectorUpdatesRange<Number>(src.local_size(), src.begin(),
1693  matrix_diagonal_inverse.begin(),
1694  start_zero, factor1, factor2,
1695  update1.begin(), update2.begin(), dst.begin());
1696  }
1697 
1698  template <typename VECTOR>
1699  struct DiagonalPreconditioner
1700  {
1701  DiagonalPreconditioner (const VECTOR &vector)
1702  :
1703  diagonal_vector(vector)
1704  {}
1705 
1706  void vmult (VECTOR &dst,
1707  const VECTOR &src) const
1708  {
1709  dst = src;
1710  dst.scale(diagonal_vector);
1711  }
1712 
1713  const VECTOR &diagonal_vector;
1714  };
1715  }
1716 }
1717 
1718 
1719 
1720 template <class MATRIX, class VECTOR>
1721 inline
1723 AdditionalData (const unsigned int degree,
1724  const double smoothing_range,
1725  const bool nonzero_starting,
1726  const unsigned int eig_cg_n_iterations,
1727  const double eig_cg_residual,
1728  const double max_eigenvalue)
1729  :
1730  degree (degree),
1731  smoothing_range (smoothing_range),
1732  nonzero_starting (nonzero_starting),
1733  eig_cg_n_iterations (eig_cg_n_iterations),
1734  eig_cg_residual (eig_cg_residual),
1735  max_eigenvalue (max_eigenvalue)
1736 {}
1737 
1738 
1739 
1740 template <class MATRIX, class VECTOR>
1741 inline
1743  :
1744  is_initialized (false)
1745 {}
1746 
1747 
1748 
1749 template <class MATRIX, class VECTOR>
1750 inline
1751 void
1753  const AdditionalData &additional_data)
1754 {
1755  matrix_ptr = &matrix;
1756  data = additional_data;
1757  if (data.matrix_diagonal_inverse.size() != matrix.m())
1758  {
1759  Assert(data.matrix_diagonal_inverse.size() == 0,
1760  ExcMessage("Matrix diagonal vector set but not sized correctly"));
1761  data.matrix_diagonal_inverse.reinit(matrix.m());
1762  for (unsigned int i=0; i<matrix.m(); ++i)
1763  data.matrix_diagonal_inverse(i) = 1./matrix.el(i,i);
1764  }
1765 
1766 
1767  // calculate largest eigenvalue using a hand-tuned CG iteration on the
1768  // matrix weighted by its diagonal. we start with a vector that consists of
1769  // ones only, weighted by the length.
1770  double max_eigenvalue, min_eigenvalue;
1771  if (data.eig_cg_n_iterations > 0)
1772  {
1773  Assert (additional_data.eig_cg_n_iterations > 2,
1774  ExcMessage ("Need to set at least two iterations to find eigenvalues."));
1775 
1776  // attach stream to SolverCG, run it with log report for eigenvalues
1777  std::ostream *old_stream = deallog.has_file() ? &deallog.get_file_stream() :
1778  static_cast<std::ostream *>(0);
1779  if (old_stream)
1780  deallog.detach();
1781 
1782  std::ostringstream log_msg;
1783  deallog.attach(log_msg);
1784 
1785  // set a very strict tolerance to force at least two iterations
1786  ReductionControl control (data.eig_cg_n_iterations, 1e-20, 1e-20);
1788  VECTOR *rhs = memory.alloc();
1789  VECTOR *dummy = memory.alloc();
1790  rhs->reinit(data.matrix_diagonal_inverse, true);
1791  dummy->reinit(data.matrix_diagonal_inverse);
1792  *rhs = 1./std::sqrt(static_cast<double>(matrix.m()));
1793 
1794  typename SolverCG<VECTOR>::AdditionalData cg_data;
1795  cg_data.compute_eigenvalues = true;
1796  SolverCG<VECTOR> solver (control, memory, cg_data);
1797  internal::PreconditionChebyshev::DiagonalPreconditioner<VECTOR>
1798  preconditioner(data.matrix_diagonal_inverse);
1799  try
1800  {
1801  solver.solve(matrix, *dummy, *rhs, preconditioner);
1802  }
1804  {
1805  }
1806  Assert(control.last_step() >= 2,
1807  ExcMessage("Could not find eigenvalues"));
1808 
1809  memory.free(dummy);
1810  memory.free(rhs);
1811 
1812  // read the log stream: grab the first and last eigenvalue
1813  std::string cg_message = log_msg.str();
1814  const std::size_t pos = cg_message.find("cg:: ");
1815  if (pos != std::string::npos)
1816  {
1817  cg_message.erase(0, pos+5);
1818  std::string first = cg_message;
1819 
1820  if (cg_message.find_first_of(" ") != std::string::npos)
1821  first.erase(cg_message.find_first_of(" "), std::string::npos);
1822  std::istringstream(first) >> min_eigenvalue;
1823 
1824  if (cg_message.find_last_of(" ") != std::string::npos)
1825  {
1826  cg_message.erase(0, cg_message.find_last_of(" ")+1);
1827  std::istringstream(cg_message) >> max_eigenvalue;
1828  }
1829  else max_eigenvalue = min_eigenvalue;
1830  }
1831  else
1832  min_eigenvalue = max_eigenvalue = 1;
1833 
1834  // reset deal.II stream
1835  deallog.detach();
1836  if (old_stream)
1837  deallog.attach(*old_stream, false);
1838 
1839  // include a safety factor since the CG method will in general not be
1840  // converged
1841  max_eigenvalue *= 1.2;
1842  }
1843  else
1844  {
1845  max_eigenvalue = data.max_eigenvalue;
1846  min_eigenvalue = data.max_eigenvalue/data.smoothing_range;
1847  }
1848 
1849  const double alpha = (data.smoothing_range > 1. ?
1850  max_eigenvalue / data.smoothing_range :
1851  std::min(0.9*max_eigenvalue, min_eigenvalue));
1852  delta = (max_eigenvalue-alpha)*0.5;
1853  theta = (max_eigenvalue+alpha)*0.5;
1854 
1855  update1.reinit (data.matrix_diagonal_inverse, true);
1856  update2.reinit (data.matrix_diagonal_inverse, true);
1857 
1858  is_initialized = true;
1859 }
1860 
1861 
1862 
1863 template <class MATRIX, class VECTOR>
1864 inline
1865 void
1867  const VECTOR &src) const
1868 {
1869  Assert (is_initialized, ExcMessage("Preconditioner not initialized"));
1870  double rhok = delta / theta, sigma = theta / delta;
1871  if (data.nonzero_starting && !dst.all_zero())
1872  {
1873  matrix_ptr->vmult (update2, dst);
1874  internal::PreconditionChebyshev::vector_updates
1875  (src, data.matrix_diagonal_inverse, false, 0., 1./theta, update1,
1876  update2, dst);
1877  }
1878  else
1879  internal::PreconditionChebyshev::vector_updates
1880  (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1,
1881  update2, dst);
1882 
1883  for (unsigned int k=0; k<data.degree; ++k)
1884  {
1885  matrix_ptr->vmult (update2, dst);
1886  const double rhokp = 1./(2.*sigma-rhok);
1887  const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta;
1888  rhok = rhokp;
1889  internal::PreconditionChebyshev::vector_updates
1890  (src, data.matrix_diagonal_inverse, false, factor1, factor2, update1,
1891  update2, dst);
1892  }
1893 }
1894 
1895 
1896 
1897 template <class MATRIX, class VECTOR>
1898 inline
1899 void
1901  const VECTOR &src) const
1902 {
1903  Assert (is_initialized, ExcMessage("Preconditioner not initialized"));
1904  double rhok = delta / theta, sigma = theta / delta;
1905  if (data.nonzero_starting && !dst.all_zero())
1906  {
1907  matrix_ptr->Tvmult (update2, dst);
1908  internal::PreconditionChebyshev::vector_updates
1909  (src, data.matrix_diagonal_inverse, false, 0., 1./theta, update1,
1910  update2, dst);
1911  }
1912  else
1913  internal::PreconditionChebyshev::vector_updates
1914  (src, data.matrix_diagonal_inverse, true, 0., 1./theta, update1,
1915  update2, dst);
1916 
1917  for (unsigned int k=0; k<data.degree; ++k)
1918  {
1919  matrix_ptr->Tvmult (update2, dst);
1920  const double rhokp = 1./(2.*sigma-rhok);
1921  const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta;
1922  rhok = rhokp;
1923  internal::PreconditionChebyshev::vector_updates
1924  (src, data.matrix_diagonal_inverse, false, factor1, factor2, update1,
1925  update2, dst);
1926  }
1927 }
1928 
1929 
1930 
1931 template <class MATRIX, class VECTOR>
1932 inline
1934 {
1935  is_initialized = false;
1936  matrix_ptr = 0;
1937  data.matrix_diagonal_inverse.reinit(0);
1938  update1.reinit(0);
1939  update2.reinit(0);
1940 }
1941 
1942 
1943 
1944 #endif // DOXYGEN
1945 
1946 DEAL_II_NAMESPACE_CLOSE
1947 
1948 #endif
void initialize(const MATRIX &matrix, const AdditionalData &additional_data=AdditionalData())
AdditionalData data
const std::vector< size_type > * inverse_permutation
Definition: precondition.h:733
void Tvmult(VECTOR &, const VECTOR &) const
void initialize(const MATRIX &A, const AdditionalData &parameters=AdditionalData())
void Tvmult_add(VECTOR &, const VECTOR &) const
void Tstep(VECTOR &x, const VECTOR &rhs) const
void Tvmult(VECTOR &, const VECTOR &) const
const function_ptr precondition
Definition: precondition.h:340
void vmult(VECTOR &dst, const VECTOR &src) const
void vmult(VECTOR &dst, const VECTOR &src) const
void Tvmult(VECTOR &dst, const VECTOR &src) const
void step(VECTOR &x, const VECTOR &rhs) const
SmartPointer< const MATRIX, PreconditionChebyshev< MATRIX, VECTOR > > matrix_ptr
::ExceptionBase & ExcMessage(std::string arg1)
void reinit(const size_type size, const bool fast=false)
void Tvmult(VECTOR &dst, const VECTOR &src) const
types::global_dof_index size_type
Definition: precondition.h:583
void Tvmult(VECTOR &, const VECTOR &) const
void vmult(VECTOR &dst, const VECTOR &src) const
void Tvmult(VECTOR &, const VECTOR &) const
double residual(VECTOR &dst, const VECTOR &src, const VECTOR &rhs) const
AdditionalData(const double relaxation=1.)
VectorMemory< VECTOR > & mem
Definition: precondition.h:890
void vmult_add(VECTOR &, const VECTOR &) const
bool has_file() const
std::ostream & get_file_stream()
virtual VECTOR * alloc()
const_iterator begin() const
SmartPointer< SOLVER, PreconditionLACSolver< SOLVER, MATRIX, PRECONDITION > > solver
Definition: precondition.h:819
void initialize(const AdditionalData &parameters)
void vmult(VECTOR &, const VECTOR &) const
const MATRIX & A
Definition: precondition.h:882
unsigned int global_dof_index
Definition: types.h:100
void Tvmult(VECTOR &, const VECTOR &) const
#define Assert(cond, exc)
Definition: exceptions.h:299
void Tvmult(VECTOR &, const VECTOR &) const
void Tvmult_add(VECTOR &, const VECTOR &) const
PreconditionedMatrix(const MATRIX &A, const PRECOND &P, VectorMemory< VECTOR > &mem)
AdditionalData(const unsigned int degree=0, const double smoothing_range=0., const bool nonzero_starting=false, const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1)
void step(VECTOR &x, const VECTOR &rhs) const
void vmult(VECTOR &, const VECTOR &) const
void(MATRIX::* function_ptr)(VECTOR &, const VECTOR &) const
Definition: precondition.h:307
size_type local_size() const
SmartPointer< const PRECONDITION, PreconditionLACSolver< SOLVER, MATRIX, PRECONDITION > > precondition
Definition: precondition.h:829
iterator begin()
void initialize(const MATRIX &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename PreconditionRelaxation< MATRIX >::AdditionalData &parameters=typename PreconditionRelaxation< MATRIX >::AdditionalData())
types::global_dof_index size_type
Definition: precondition.h:916
SmartPointer< const MATRIX, PreconditionRelaxation< MATRIX > > A
Definition: precondition.h:396
void initialize(SOLVER &, const MATRIX &, const PRECONDITION &)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const_iterator end() const
void attach(std::ostream &o, const bool print_job_id=true)
types::global_dof_index size_type
Definition: precondition.h:685
void step(VECTOR &x, const VECTOR &rhs) const
void vmult(VECTOR &, const VECTOR &) const
const MATRIX & matrix
Definition: precondition.h:334
const PRECOND & P
Definition: precondition.h:886
std::vector< std::size_t > pos_right_of_diagonal
Definition: precondition.h:643
const std::vector< size_type > * permutation
Definition: precondition.h:728
PreconditionUseMatrix(const MATRIX &M, const function_ptr method)
void initialize(const MATRIX &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
void Tstep(VECTOR &x, const VECTOR &rhs) const
void vmult_add(VECTOR &, const VECTOR &) const
::ExceptionBase & ExcNotInitialized()
void vmult(VECTOR &, const VECTOR &) const
AdditionalData(const double relaxation=1.)
void vmult(VECTOR &, const VECTOR &) const
void detach()
void Tstep(VECTOR &x, const VECTOR &rhs) const
SmartPointer< const MATRIX, PreconditionLACSolver< SOLVER, MATRIX, PRECONDITION > > matrix
Definition: precondition.h:824
virtual void free(const VECTOR *const)
void vmult(VECTOR &, const VECTOR &) const
void vmult(VECTOR &, const VECTOR &) const
void initialize(const MATRIX &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionRelaxation< MATRIX > BaseClass
Definition: precondition.h:588