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_;

  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