Ifpack2_Details_Chebyshev_decl.hpp

Go to the documentation of this file.

/*
//@HEADER
// ***********************************************************************
//
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2009) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@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 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_;

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

  bool zeroStartingSolution_;

  bool assumeMatrixUnchanged_;

  bool textbookAlgorithm_;

  bool computeMaxResNorm_;

  Teuchos::RCP<Teuchos::FancyOStream> out_;

  bool debug_;



  void checkConstructorInput () const;

  void checkInputMatrix () const;

  void reset ();

  Teuchos::RCP<MV> makeTempMultiVector (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
  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);

  void computeInitialGuessForPowerMethod (V& x, const bool nonnegativeRealParts = false) const;

  ST
  powerMethodWithInitGuess (const op_type& A, const V& D_inv, const int numIters, V& x);

  ST
  powerMethod (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