Ifpack2_Details_Chebyshev_decl.hpp

Go to the documentation of this file.

/*
//@HEADER
// ***********************************************************************
//
//       Ifpack2: Tempated 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
*/

// ***********************************************************************
//
//      Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2004) 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.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************

#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 {
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<const V> D_;

  Teuchos::ArrayRCP<size_t> diagOffsets_;

  bool savedDiagOffsets_;




  Teuchos::RCP<MV> V_;

  Teuchos::RCP<MV> W_;

  ST computedLambdaMax_;

  ST computedLambdaMin_;




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




  void checkConstructorInput () const;

  void checkInputMatrix () const;

  void reset ();

  void
  makeTempMultiVectors (Teuchos::RCP<MV>& V1,
                        Teuchos::RCP<MV>& W,
                        const 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);

  static 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