Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_Details_Chebyshev_decl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4 //
5 // Copyright 2009 NTESS and the Ifpack2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
11 #define IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
12 
19 
20 #include "Ifpack2_ConfigDefs.hpp"
22 #include "Teuchos_Describable.hpp"
23 #include "Tpetra_CrsMatrix.hpp"
24 
25 namespace Ifpack2 {
26 namespace Details {
27 
28 #ifndef DOXYGEN_SHOULD_SKIP_THIS
29 template <class TpetraOperatorType>
30 class ChebyshevKernel; // forward declaration
31 #endif // DOXYGEN_SHOULD_SKIP_THIS
32 
74 template <class ScalarType, class MV>
76  public:
78 
79 
81  typedef ScalarType ST;
85  typedef typename STS::magnitudeType MT;
87  typedef Tpetra::Operator<typename MV::scalar_type,
88  typename MV::local_ordinal_type,
89  typename MV::global_ordinal_type,
90  typename MV::node_type>
93  typedef Tpetra::RowMatrix<typename MV::scalar_type,
94  typename MV::local_ordinal_type,
95  typename MV::global_ordinal_type,
96  typename MV::node_type>
99  typedef Tpetra::Vector<typename MV::scalar_type,
100  typename MV::local_ordinal_type,
101  typename MV::global_ordinal_type,
102  typename MV::node_type>
103  V;
105  typedef Tpetra::Map<typename MV::local_ordinal_type,
106  typename MV::global_ordinal_type,
107  typename MV::node_type>
110 
119 
130 
240 
241  void setZeroStartingSolution(bool zeroStartingSolution) { zeroStartingSolution_ = zeroStartingSolution; }
242 
276  void compute();
277 
294  MT apply(const MV& B, MV& X);
295 
296  ST getLambdaMaxForApply() const;
297 
300 
307 
309  bool hasTransposeApply() const;
310 
312  void print(std::ostream& out);
313 
315 
317 
319  std::string description() const;
320 
322  void
324  const Teuchos::EVerbosityLevel verbLevel =
327  private:
329 
330 
338 
341 
353 
355  typedef Kokkos::View<size_t*, typename MV::node_type::device_type> offsets_type;
356 
362  offsets_type diagOffsets_;
363 
371  bool savedDiagOffsets_;
372 
374 
376 
380  Teuchos::RCP<MV> W_;
381  Teuchos::RCP<MV> W2_;
382 
388  ST computedLambdaMax_;
389 
395  ST computedLambdaMin_;
396 
398 
400 
403  ST lambdaMaxForApply_;
404 
411  MT boostFactor_;
414  ST lambdaMinForApply_;
417  ST eigRatioForApply_;
418 
420 
422 
427  Teuchos::RCP<const V> userInvDiag_;
428 
432  ST userLambdaMax_;
433 
437  ST userLambdaMin_;
438 
442  ST userEigRatio_;
443 
448  ST minDiagVal_;
449 
451  int numIters_;
452 
454  int eigMaxIters_;
455 
457  MT eigRelTolerance_;
458 
460  bool eigKeepVectors_;
461 
463  std::string eigenAnalysisType_;
464 
466  Teuchos::RCP<V> eigVector_;
467  Teuchos::RCP<V> eigVector2_;
468 
470  int eigNormalizationFreq_;
471 
473  bool zeroStartingSolution_;
474 
481  bool assumeMatrixUnchanged_;
482 
484  std::string chebyshevAlgorithm_;
485 
487  bool computeMaxResNorm_;
488 
490  bool computeSpectralRadius_;
491 
494  bool ckUseNativeSpMV_;
495 
499  bool preAllocateTempVector_;
500 
505 
507  bool debug_;
508 
510 
512 
514  void checkConstructorInput() const;
515 
517  void checkInputMatrix() const;
518 
526  void reset();
527 
533  Teuchos::RCP<MV> makeTempMultiVector(const MV& B);
534 
541  Teuchos::RCP<MV> makeSecondTempMultiVector(const MV& B);
542 
544  void
545  firstIterationWithZeroStartingSolution(MV& W,
546  const ScalarType& alpha,
547  const V& D_inv,
548  const MV& B,
549  MV& X);
550 
552  static void
553  computeResidual(MV& R, const MV& B, const op_type& A, const MV& X,
554  const Teuchos::ETransp mode = Teuchos::NO_TRANS);
555 
561  static void solve(MV& Z, const V& D_inv, const MV& R);
562 
568  static void solve(MV& Z, const ST alpha, const V& D_inv, const MV& R);
569 
578  makeInverseDiagonal(const row_matrix_type& A, const bool useDiagOffsets = false) const;
579 
589  Teuchos::RCP<V> makeRangeMapVector(const Teuchos::RCP<V>& D) const;
590 
593  makeRangeMapVectorConst(const Teuchos::RCP<const V>& D) const;
594 
613  void
614  textbookApplyImpl(const op_type& A,
615  const MV& B,
616  MV& X,
617  const int numIters,
618  const ST lambdaMax,
619  const ST lambdaMin,
620  const ST eigRatio,
621  const V& D_inv) const;
622 
644  void
645  fourthKindApplyImpl(const op_type& A,
646  const MV& B,
647  MV& X,
648  const int numIters,
649  const ST lambdaMax,
650  const V& D_inv);
651 
674  void
675  ifpackApplyImpl(const op_type& A,
676  const MV& B,
677  MV& X,
678  const int numIters,
679  const ST lambdaMax,
680  const ST lambdaMin,
681  const ST eigRatio,
682  const V& D_inv);
683 
696  ST cgMethodWithInitGuess(const op_type& A, const V& D_inv, const int numIters, V& x);
697 
707  ST cgMethod(const op_type& A, const V& D_inv, const int numIters);
708 
710  static MT maxNormInf(const MV& X);
711 
712  // mfh 22 Jan 2013: This code builds correctly, but does not
713  // converge. I'm leaving it in place in case someone else wants to
714  // work on it.
715 #if 0
716  void
739  mlApplyImpl (const op_type& A,
740  const MV& B,
741  MV& X,
742  const int numIters,
743  const ST lambdaMax,
744  const ST lambdaMin,
745  const ST eigRatio,
746  const V& D_inv)
747  {
748  const ST zero = Teuchos::as<ST> (0);
749  const ST one = Teuchos::as<ST> (1);
750  const ST two = Teuchos::as<ST> (2);
751 
752  MV pAux (B.getMap (), B.getNumVectors ()); // Result of A*X
753  MV dk (B.getMap (), B.getNumVectors ()); // Solution update
754  MV R (B.getMap (), B.getNumVectors ()); // Not in original ML; need for B - pAux
755 
756  ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
757  ST alpha = lambdaMax / eigRatio;
758 
759  ST delta = (beta - alpha) / two;
760  ST theta = (beta + alpha) / two;
761  ST s1 = theta / delta;
762  ST rhok = one / s1;
763 
764  // Diagonal: ML replaces entries containing 0 with 1. We
765  // shouldn't have any entries like that in typical test problems,
766  // so it's OK not to do that here.
767 
768  // The (scaled) matrix is the identity: set X = D_inv * B. (If A
769  // is the identity, then certainly D_inv is too. D_inv comes from
770  // A, so if D_inv * A is the identity, then we still need to apply
771  // the "preconditioner" D_inv to B as well, to get X.)
772  if (lambdaMin == one && lambdaMin == lambdaMax) {
773  solve (X, D_inv, B);
774  return;
775  }
776 
777  // The next bit of code is a direct translation of code from ML's
778  // ML_Cheby function, in the "normal point scaling" section, which
779  // is in lines 7365-7392 of ml_smoother.c.
780 
781  if (! zeroStartingSolution_) {
782  // dk = (1/theta) * D_inv * (B - (A*X))
783  A.apply (X, pAux); // pAux = A * X
784  R = B;
785  R.update (-one, pAux, one); // R = B - pAux
786  dk.elementWiseMultiply (one/theta, D_inv, R, zero); // dk = (1/theta)*D_inv*R
787  X.update (one, dk, one); // X = X + dk
788  } else {
789  dk.elementWiseMultiply (one/theta, D_inv, B, zero); // dk = (1/theta)*D_inv*B
790  X = dk;
791  }
792 
793  ST rhokp1, dtemp1, dtemp2;
794  for (int k = 0; k < numIters-1; ++k) {
795  A.apply (X, pAux);
796  rhokp1 = one / (two*s1 - rhok);
797  dtemp1 = rhokp1*rhok;
798  dtemp2 = two*rhokp1/delta;
799  rhok = rhokp1;
800 
801  R = B;
802  R.update (-one, pAux, one); // R = B - pAux
803  // dk = dtemp1 * dk + dtemp2 * D_inv * (B - pAux)
804  dk.elementWiseMultiply (dtemp2, D_inv, B, dtemp1);
805  X.update (one, dk, one); // X = X + dk
806  }
807  }
808 #endif // 0
809 
810 };
811 
812 } // namespace Details
813 } // namespace Ifpack2
814 
815 #endif // IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:340
MT apply(const MV &B, MV &X)
Solve Ax=b for x with Chebyshev iteration with left diagonal scaling.
Definition: Ifpack2_Details_Chebyshev_def.hpp:999
Tpetra::Operator< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > op_type
Specialization of Tpetra::Operator.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:91
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print a description of the Chebyshev solver to out.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1661
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:755
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1597
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:275
Tpetra::Map< typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > map_type
Specialization of Tpetra::Map.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:108
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1049
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:85
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:81
Tpetra::RowMatrix< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > row_matrix_type
Specialization of Tpetra::RowMatrix.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:97
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:75
static const EVerbosityLevel verbLevel_default
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1641
Teuchos::ScalarTraits< ScalarType > STS
Traits class for ST.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:83
bool hasTransposeApply() const
Whether it&#39;s possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1603
void compute()
(Re)compute the left scaling D_inv, and estimate min and max eigenvalues of D_inv * A...
Definition: Ifpack2_Details_Chebyshev_def.hpp:791
Tpetra::Vector< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > V
Specialization of Tpetra::Vector.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:103