Ifpack2_Details_Chebyshev_decl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//
// Copyright 2009 NTESS and the Ifpack2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
#define IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
 
 
#include "Ifpack2_ConfigDefs.hpp"
#include "Teuchos_VerbosityLevel.hpp"
#include "Teuchos_Describable.hpp"
#include "Tpetra_CrsMatrix.hpp"
 
namespace Ifpack2 {
namespace Details {
 
#ifndef DOXYGEN_SHOULD_SKIP_THIS
template <class TpetraOperatorType>
class ChebyshevKernel;  // forward declaration
#endif                  // DOXYGEN_SHOULD_SKIP_THIS
 
template <class ScalarType, class MV>
class Chebyshev : public Teuchos::Describable {
 public:
 
 
  typedef ScalarType ST;
  typedef Teuchos::ScalarTraits<ScalarType> STS;
  typedef typename STS::magnitudeType MT;
  typedef Tpetra::Operator<typename MV::scalar_type,
                           typename MV::local_ordinal_type,
                           typename MV::global_ordinal_type,
                           typename MV::node_type>
      op_type;
  typedef Tpetra::RowMatrix<typename MV::scalar_type,
                            typename MV::local_ordinal_type,
                            typename MV::global_ordinal_type,
                            typename MV::node_type>
      row_matrix_type;
  typedef Tpetra::Vector<typename MV::scalar_type,
                         typename MV::local_ordinal_type,
                         typename MV::global_ordinal_type,
                         typename MV::node_type>
      V;
  typedef Tpetra::Map<typename MV::local_ordinal_type,
                      typename MV::global_ordinal_type,
                      typename MV::node_type>
      map_type;
 
  Chebyshev(Teuchos::RCP<const row_matrix_type> A);
 
  Chebyshev(Teuchos::RCP<const row_matrix_type> A, Teuchos::ParameterList& params);
 
  void setParameters(Teuchos::ParameterList& plist);
 
  void setZeroStartingSolution(bool zeroStartingSolution) { zeroStartingSolution_ = zeroStartingSolution; }
 
  void compute();
 
  MT apply(const MV& B, MV& X);
 
  ST getLambdaMaxForApply() const;
 
  Teuchos::RCP<const row_matrix_type> getMatrix() const;
 
  void setMatrix(const Teuchos::RCP<const row_matrix_type>& A);
 
  bool hasTransposeApply() const;
 
  void print(std::ostream& out);
 
 
 
  std::string description() const;
 
  void
  describe(Teuchos::FancyOStream& out,
           const Teuchos::EVerbosityLevel verbLevel =
               Teuchos::Describable::verbLevel_default) const;
 private:
 
 
  Teuchos::RCP<const row_matrix_type> A_;
 
  Teuchos::RCP<ChebyshevKernel<op_type> > ck_;
 
  Teuchos::RCP<const V> D_;
 
  typedef Kokkos::View<size_t*, typename MV::node_type::device_type> offsets_type;
 
  offsets_type diagOffsets_;
 
  bool savedDiagOffsets_;
 
 
 
  Teuchos::RCP<MV> W_;
  Teuchos::RCP<MV> W2_;
 
  ST computedLambdaMax_;
 
  ST computedLambdaMin_;
 
 
 
  ST lambdaMaxForApply_;
 
  MT boostFactor_;
  ST lambdaMinForApply_;
  ST eigRatioForApply_;
 
 
 
  Teuchos::RCP<const V> userInvDiag_;
 
  ST userLambdaMax_;
 
  ST userLambdaMin_;
 
  ST userEigRatio_;
 
  ST minDiagVal_;
 
  int numIters_;
 
  int eigMaxIters_;
 
  MT eigRelTolerance_;
 
  bool eigKeepVectors_;
 
  std::string eigenAnalysisType_;
 
  Teuchos::RCP<V> eigVector_;
  Teuchos::RCP<V> eigVector2_;
 
  int eigNormalizationFreq_;
 
  bool zeroStartingSolution_;
 
  bool assumeMatrixUnchanged_;
 
  std::string chebyshevAlgorithm_;
 
  bool computeMaxResNorm_;
 
  bool computeSpectralRadius_;
 
  bool ckUseNativeSpMV_;
 
  bool preAllocateTempVector_;
 
  Teuchos::RCP<Teuchos::FancyOStream> out_;
 
  bool debug_;
 
 
 
  void checkConstructorInput() const;
 
  void checkInputMatrix() const;
 
  void reset();
 
  Teuchos::RCP<MV> makeTempMultiVector(const MV& B);
 
  Teuchos::RCP<MV> makeSecondTempMultiVector(const MV& B);
 
  void
  firstIterationWithZeroStartingSolution(MV& W,
                                         const ScalarType& alpha,
                                         const V& D_inv,
                                         const MV& B,
                                         MV& X);
 
  static void
  computeResidual(MV& R, const MV& B, const op_type& A, const MV& X,
                  const Teuchos::ETransp mode = Teuchos::NO_TRANS);
 
  static void solve(MV& Z, const V& D_inv, const MV& R);
 
  static void solve(MV& Z, const ST alpha, const V& D_inv, const MV& R);
 
  Teuchos::RCP<const V>
  makeInverseDiagonal(const row_matrix_type& A, const bool useDiagOffsets = false) const;
 
  Teuchos::RCP<V> makeRangeMapVector(const Teuchos::RCP<V>& D) const;
 
  Teuchos::RCP<const V>
  makeRangeMapVectorConst(const Teuchos::RCP<const V>& D) const;
 
  void
  textbookApplyImpl(const op_type& A,
                    const MV& B,
                    MV& X,
                    const int numIters,
                    const ST lambdaMax,
                    const ST lambdaMin,
                    const ST eigRatio,
                    const V& D_inv) const;
 
  void
  fourthKindApplyImpl(const op_type& A,
                      const MV& B,
                      MV& X,
                      const int numIters,
                      const ST lambdaMax,
                      const V& D_inv);
 
  void
  ifpackApplyImpl(const op_type& A,
                  const MV& B,
                  MV& X,
                  const int numIters,
                  const ST lambdaMax,
                  const ST lambdaMin,
                  const ST eigRatio,
                  const V& D_inv);
 
  ST cgMethodWithInitGuess(const op_type& A, const V& D_inv, const int numIters, V& x);
 
  ST cgMethod(const op_type& A, const V& D_inv, const int numIters);
 
  static MT maxNormInf(const MV& X);
 
  // mfh 22 Jan 2013: This code builds correctly, but does not
  // converge.  I'm leaving it in place in case someone else wants to
  // work on it.
#if 0
  void
  mlApplyImpl (const op_type& A,
               const MV& B,
               MV& X,
               const int numIters,
               const ST lambdaMax,
               const ST lambdaMin,
               const ST eigRatio,
               const V& D_inv)
  {
    const ST zero = Teuchos::as<ST> (0);
    const ST one = Teuchos::as<ST> (1);
    const ST two = Teuchos::as<ST> (2);
 
    MV pAux (B.getMap (), B.getNumVectors ()); // Result of A*X
    MV dk (B.getMap (), B.getNumVectors ()); // Solution update
    MV R (B.getMap (), B.getNumVectors ()); // Not in original ML; need for B - pAux
 
    ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
    ST alpha = lambdaMax / eigRatio;
 
    ST delta = (beta - alpha) / two;
    ST theta = (beta + alpha) / two;
    ST s1 = theta / delta;
    ST rhok = one / s1;
 
    // Diagonal: ML replaces entries containing 0 with 1.  We
    // shouldn't have any entries like that in typical test problems,
    // so it's OK not to do that here.
 
    // The (scaled) matrix is the identity: set X = D_inv * B.  (If A
    // is the identity, then certainly D_inv is too.  D_inv comes from
    // A, so if D_inv * A is the identity, then we still need to apply
    // the "preconditioner" D_inv to B as well, to get X.)
    if (lambdaMin == one && lambdaMin == lambdaMax) {
      solve (X, D_inv, B);
      return;
    }
 
    // The next bit of code is a direct translation of code from ML's
    // ML_Cheby function, in the "normal point scaling" section, which
    // is in lines 7365-7392 of ml_smoother.c.
 
    if (! zeroStartingSolution_) {
      // dk = (1/theta) * D_inv * (B - (A*X))
      A.apply (X, pAux); // pAux = A * X
      R = B;
      R.update (-one, pAux, one); // R = B - pAux
      dk.elementWiseMultiply (one/theta, D_inv, R, zero); // dk = (1/theta)*D_inv*R
      X.update (one, dk, one); // X = X + dk
    } else {
      dk.elementWiseMultiply (one/theta, D_inv, B, zero); // dk = (1/theta)*D_inv*B
      X = dk;
    }
 
    ST rhokp1, dtemp1, dtemp2;
    for (int k = 0; k < numIters-1; ++k) {
      A.apply (X, pAux);
      rhokp1 = one / (two*s1 - rhok);
      dtemp1 = rhokp1*rhok;
      dtemp2 = two*rhokp1/delta;
      rhok = rhokp1;
 
      R = B;
      R.update (-one, pAux, one); // R = B - pAux
      // dk = dtemp1 * dk + dtemp2 * D_inv * (B - pAux)
      dk.elementWiseMultiply (dtemp2, D_inv, B, dtemp1);
      X.update (one, dk, one); // X = X + dk
    }
  }
#endif  // 0
};
 
}  // namespace Details
}  // namespace Ifpack2
 
#endif  // IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP