Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Details_Chebyshev_def.hpp
Go to the documentation of this file.
00001 // ***********************************************************************
00002 // 
00003 //      Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
00004 //                 Copyright (2004) Sandia Corporation
00005 // 
00006 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00007 // license for use of this work by or on behalf of the U.S. Government.
00008 // 
00009 // This library is free software; you can redistribute it and/or modify
00010 // it under the terms of the GNU Lesser General Public License as
00011 // published by the Free Software Foundation; either version 2.1 of the
00012 // License, or (at your option) any later version.
00013 //  
00014 // This library is distributed in the hope that it will be useful, but
00015 // WITHOUT ANY WARRANTY; without even the implied warranty of
00016 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 // Lesser General Public License for more details.
00018 //  
00019 // You should have received a copy of the GNU Lesser General Public
00020 // License along with this library; if not, write to the Free Software
00021 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00022 // USA
00023 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00024 // 
00025 // ***********************************************************************
00026 
00027 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
00028 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
00029 
00036 
00037 #include "Ifpack2_Details_Chebyshev_decl.hpp"
00038 #include <cmath>
00039 
00040 // Uncommit the #define line below if you want Chebyshev to do extra
00041 // debug checking and generate a lot of debugging output to stderr (on
00042 // all processes in the communicator).  Even if you uncomment this
00043 // line, the debugging code will only be enabled if the CMake option
00044 // Teuchos_ENABLE_DEBUG was set to ON when configuring Trilinos.
00045 //#define IFPACK_DETAILS_CHEBYSHEV_DEBUG 1
00046 
00047 namespace Ifpack2 {
00048 namespace Details {
00049 
00050 template<class ScalarType, class MV, class MAT>
00051 void
00052 Chebyshev<ScalarType, MV, MAT>::
00053 checkConstructorInput () const
00054 {
00055   TEUCHOS_TEST_FOR_EXCEPTION(STS::isComplex, std::logic_error,
00056     "Ifpack2::Details::Chebyshev: This class' implementation of Chebyshev "
00057     "iteration only works for real-valued, symmetric positive definite "
00058     "matrices.  However, you instantiated this class for ScalarType=" 
00059     << Teuchos::TypeNameTraits<ScalarType>::name () << ", which is a complex-"
00060     "valued type.  While this may be algorithmically correct if all of the "
00061     "complex numbers in the matrix have zero imaginary part, we forbid using "
00062     "complex ScalarType altogether in order to remind you of the limitations "
00063     "of our implementation (and of the algorithm itself).");
00064   TEUCHOS_TEST_FOR_EXCEPTION(A_.is_null (), std::invalid_argument,
00065     "Ifpack2::Details::Chebyshev: Input matrix to constructor is null.");
00066   TEUCHOS_TEST_FOR_EXCEPTION(
00067     A_->getGlobalNumRows() != A_->getGlobalNumCols(), 
00068     std::invalid_argument,
00069     "Ifpack2::Details::Chebyshev: The input matrix A must be square.  "
00070     "A has " << A_->getGlobalNumRows() << " rows and " 
00071     << A_->getGlobalNumCols() << " columns.");
00072 
00073   // In a debug build, test that the domain and range Maps of the
00074   // matrix are the same.
00075 #ifdef HAVE_TEUCHOS_DEBUG
00076   Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
00077   Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
00078 
00079   // The relation 'isSameAs' is transitive.  It's also a collective,
00080   // so we don't have to do a "shared" test for exception (i.e., a
00081   // global reduction on the test value).  Do the raw pointer
00082   // comparison first, to avoid a collective in the common case.
00083   TEUCHOS_TEST_FOR_EXCEPTION(
00084      domainMap.getRawPtr () != rangeMap.getRawPtr () ||
00085      ! domainMap->isSameAs (*rangeMap),
00086      std::invalid_argument,
00087      "Ifpack2::Details::Chebyshev: The domain Map and range Map of the matrix "
00088      "must be the same (in the sense of isSameAs()).  We only check for this "
00089      "if Trilinos was built with the CMake configuration option Teuchos_ENABLE_"
00090      "DEBUG set to ON.");
00091 
00092 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00093   // mfh 30 Jan 2013: Make sure that A is not bogus.
00094   const MT frobNormA = A_->getFrobeniusNorm ();
00095   TEUCHOS_TEST_FOR_EXCEPTION(
00096     Teuchos::ScalarTraits<MT>::isnaninf (frobNormA),
00097     std::invalid_argument,
00098     "Ifpack2::Details::Chebyshev: Frobenius norm of A is " 
00099     << frobNormA << ", which is Inf or NaN.");
00100 
00101   V X (domainMap);
00102   V Y (rangeMap);
00103   X.putScalar (STS::one ());
00104   A_->apply (X, Y);
00105   const MT infNorm = Y.normInf ();
00106   TEUCHOS_TEST_FOR_EXCEPTION(
00107     Teuchos::ScalarTraits<MT>::isnaninf (infNorm),
00108     std::invalid_argument,
00109     "Ifpack2::Details::Chebyshev: inf-norm of A*[1,...,1]^T is " 
00110     << infNorm << ", which is Inf or NaN.");
00111 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00112 #endif // HAVE_TEUCHOS_DEBUG
00113 }
00114 
00115 template<class ScalarType, class MV, class MAT>
00116 Chebyshev<ScalarType, MV, MAT>::
00117 Chebyshev (Teuchos::RCP<const MAT> A) : 
00118   A_ (A),
00119   computedLambdaMax_ (STS::nan ()),
00120   computedLambdaMin_ (STS::nan ()),
00121   lambdaMax_ (STS::nan ()),
00122   lambdaMin_ (STS::nan ()),
00123   eigRatio_ (Teuchos::as<ST> (30)),
00124   minDiagVal_ (STS::eps ()),
00125   numIters_ (1),
00126   eigMaxIters_ (10),
00127   zeroStartingSolution_ (true),
00128   textbookAlgorithm_ (false),
00129   computeMaxResNorm_ (false)
00130 {
00131   checkConstructorInput ();
00132 }
00133 
00134 template<class ScalarType, class MV, class MAT>
00135 Chebyshev<ScalarType, MV, MAT>::
00136 Chebyshev (Teuchos::RCP<const MAT> A, Teuchos::ParameterList& params) : 
00137   A_ (A), 
00138   computedLambdaMax_ (STS::nan ()),
00139   computedLambdaMin_ (STS::nan ()),
00140   lambdaMax_ (STS::nan ()),
00141   lambdaMin_ (STS::nan ()),
00142   eigRatio_ (Teuchos::as<ST> (30)),
00143   minDiagVal_ (STS::eps ()),
00144   numIters_ (1),
00145   eigMaxIters_ (10),
00146   zeroStartingSolution_ (true),
00147   textbookAlgorithm_ (false),
00148   computeMaxResNorm_ (false)
00149 {
00150   checkConstructorInput ();
00151   setParameters (params);
00152 }
00153 
00154 template<class ScalarType, class MV, class MAT>
00155 void
00156 Chebyshev<ScalarType, MV, MAT>::
00157 setParameters (Teuchos::ParameterList& plist) {
00158   using Teuchos::RCP;
00159   using Teuchos::rcp;
00160   using Teuchos::rcp_const_cast;
00161   using Teuchos::rcpFromRef;
00162 
00163   // Note to developers: The logic for this method is complicated,
00164   // because we want to accept Ifpack and ML parameters whenever
00165   // possible, but we don't want to add their default values to the
00166   // user's ParameterList.  That's why we do all the isParameter()
00167   // checks, instead of using the two-argument version of get()
00168   // everywhere.  The min and max eigenvalue parameters are also a
00169   // special case, because we decide whether or not to do eigenvalue
00170   // analysis based on whether the user supplied the max eigenvalue.
00171 
00172   // Default values of all the parameters.
00173   const ST defaultLambdaMax = STS::nan ();
00174   const ST defaultLambdaMin = STS::nan ();
00175   // 30 is Ifpack::Chebyshev's default.  ML has a different default
00176   // eigRatio for smoothers and the coarse-grid solve (if using
00177   // Chebyshev for that).  The former uses 20; the latter uses 30.
00178   // We're testing smoothers here, so use 20.  (However, if you give
00179   // ML an Epetra matrix, it will use Ifpack for Chebyshev, in which
00180   // case it would defer to Ifpack's default settings.)
00181   const ST defaultEigRatio = Teuchos::as<ST> (30);
00182   const ST defaultMinDiagVal = STS::eps ();
00183   const int defaultNumIters = 1;
00184   const int defaultEigMaxIters = 10;
00185   const bool defaultZeroStartingSolution = true; // Ifpack::Chebyshev default
00186   const bool defaultTextbookAlgorithm = false;
00187   const bool defaultComputeMaxResNorm = false;
00188 
00189   // We'll set the instance data transactionally, after all reads
00190   // from the ParameterList.  That way, if any of the ParameterList
00191   // reads fail (e.g., due to the wrong parameter type), we will not
00192   // have left the instance data in a half-changed state.
00193   RCP<const V> userInvDiag;
00194   ST lambdaMax = defaultLambdaMax;
00195   ST lambdaMin = defaultLambdaMin;
00196   ST eigRatio = defaultEigRatio;
00197   ST minDiagVal = defaultMinDiagVal;
00198   int numIters = defaultNumIters;
00199   int eigMaxIters = defaultEigMaxIters;
00200   bool zeroStartingSolution = defaultZeroStartingSolution;
00201   bool textbookAlgorithm = defaultTextbookAlgorithm;
00202   bool computeMaxResNorm = defaultComputeMaxResNorm;
00203 
00204   // Fetch the parameters from the ParameterList.  Defer all
00205   // externally visible side effects until we have finished all
00206   // ParameterList interaction.  This makes the method satisfy the
00207   // strong exception guarantee.
00208 
00209   // Get the user-supplied inverse diagonal.
00210   //
00211   // Check for a raw pointer (const V* or V*), for Ifpack
00212   // compatibility, as well as for an RCP<const V> or RCP<V>.  We'll
00213   // copy the vector anyway, so it doesn't matter whether it's const
00214   // or nonconst.
00215   if (plist.isParameter ("chebyshev: operator inv diagonal")) {
00216     try { // Could the type be const V*?
00217       const V* rawUserInvDiag = plist.get<const V*> ("chebyshev: operator inv diagonal");
00218       // It's OK to have a nonowning reference; we'll copy the vector anyway.
00219       userInvDiag = rcp (rawUserInvDiag, false);
00220     } catch (Teuchos::Exceptions::InvalidParameterType&) {
00221     }
00222     if (userInvDiag.is_null ()) {
00223       try { // Could the type be V*?
00224   V* rawUserInvDiag = plist.get<V*> ("chebyshev: operator inv diagonal");
00225   // It's OK to have a nonowning reference; we'll copy the vector anyway.
00226   userInvDiag = rcp (const_cast<const V*> (rawUserInvDiag), false);
00227       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00228       }
00229     }
00230     if (userInvDiag.is_null ()) {
00231       try { // Could the type be RCP<const V>?
00232   userInvDiag = plist.get<RCP<const V> > ("chebyshev: operator inv diagonal");
00233       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00234       }
00235     }
00236     if (userInvDiag.is_null ()) {
00237       try { // Could the type be RCP<const V>?
00238   RCP<V> userInvDiagNonconst = plist.get<RCP<V> > ("chebyshev: operator inv diagonal");
00239   userInvDiag = rcp_const_cast<const V> (userInvDiagNonconst);
00240       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00241       }
00242     }
00243     // We've tried all the possible types.  If we got something, then
00244     // make a range Map copy of it.
00245     if (! userInvDiag.is_null ()) {
00246       userInvDiag = makeRangeMapVector (*userInvDiag);
00247     }
00248   }
00249 
00250   // Don't fill in defaults for the max or min eigenvalue, because
00251   // this class uses the existence of those parameters to determine
00252   // whether it should do eigenanalysis.
00253   if (plist.isParameter ("chebyshev: max eigenvalue")) {
00254     lambdaMax = plist.get<ST> ("chebyshev: max eigenvalue");
00255   }
00256   if (plist.isParameter ("chebyshev: min eigenvalue")) {
00257     lambdaMin = plist.get<ST> ("chebyshev: min eigenvalue");
00258   }
00259 
00260   // Only fill in Ifpack2's name for the default parameter, not ML's.
00261   if (plist.isParameter ("smoother: Chebyshev alpha")) { // ML compatibility
00262     eigRatio = plist.get<ST> ("smoother: Chebyshev alpha");
00263   }
00264   // Ifpack2's name overrides ML's name.
00265   eigRatio = plist.get ("chebyshev: ratio eigenvalue", eigRatio);
00266 
00267   // Same name in Ifpack2 and Ifpack.
00268   minDiagVal = plist.get ("chebyshev: min diagonal value", minDiagVal);
00269 
00270   // Only fill in Ifpack2's name, not ML's or Ifpack's.
00271   if (plist.isParameter ("smoother: sweeps")) { // ML compatibility
00272     numIters = plist.get<int> ("smoother: sweeps");
00273   } // Ifpack's name overrides ML's name.
00274   if (plist.isParameter ("relaxation: sweeps")) { // Ifpack compatibility
00275     numIters = plist.get<int> ("relaxation: sweeps");
00276   } // Ifpack2's name overrides Ifpack's name.
00277   numIters = plist.get ("chebyshev: degree", numIters);
00278 
00279   // The last parameter name overrides the first.
00280   if (plist.isParameter ("eigen-analysis: iterations")) { // ML compatibility
00281     eigMaxIters = plist.get<int> ("eigen-analysis: iterations");
00282   } // Ifpack2's name overrides ML's name.
00283   eigMaxIters = plist.get ("chebyshev: eigenvalue max iterations", eigMaxIters);
00284 
00285   zeroStartingSolution = plist.get ("chebyshev: zero starting solution", zeroStartingSolution);
00286 
00287   // We don't want to fill these parameters in, because they shouldn't
00288   // be visible to Ifpack2::Chebyshev users.
00289   if (plist.isParameter ("chebyshev: textbook algorithm")) {
00290     textbookAlgorithm = plist.get<bool> ("chebyshev: textbook algorithm");
00291   }
00292   if (plist.isParameter ("chebyshev: compute max residual norm")) {
00293     computeMaxResNorm = plist.get<bool> ("chebyshev: compute max residual norm");
00294   }
00295 
00296   // Test for Ifpack parameters that we won't ever implement here.
00297   // Be careful to use the one-argument version of get(), since the
00298   // two-argment version adds the parameter if it's not there.
00299   TEUCHOS_TEST_FOR_EXCEPTION(
00300     plist.isParameter ("chebyshev: use block mode") && 
00301     ! plist.get<bool> ("chebyshev: use block mode"),
00302     std::invalid_argument,
00303     "Ifpack2::Chebyshev requires that if you supply the Ifpack parameter "
00304     "\"chebyshev: use block mode\", it must be set to false.  Ifpack2's "
00305     "Chebyshev does not implement Ifpack's block mode.");
00306   TEUCHOS_TEST_FOR_EXCEPTION(
00307     plist.isParameter ("chebyshev: solve normal equations") && 
00308     ! plist.get<bool> ("chebyshev: solve normal equations"),
00309     std::invalid_argument,
00310     "Ifpack2::Chebyshev does not and will never implement the Ifpack "
00311     "parameter \"chebyshev: solve normal equations\".  If you want to solve "
00312     "the normal equations, construct a Tpetra::Operator that implements "
00313     "A^* A, and use Chebyshev to solve A^* A x = A^* b.");
00314 
00315   // Test for Ifpack parameters that we haven't implemented yet.
00316   //
00317   // For now, we only check that this ML parameter, if provided, has
00318   // the one value that we support.  We consider other values "invalid
00319   // arguments" rather than "logic errors," because Ifpack does not
00320   // implement eigenanalyses other than the power method.
00321   std::string eigenAnalysisType ("power-method");
00322   if (plist.isParameter ("eigen-analysis: type")) {
00323     eigenAnalysisType = plist.get<std::string> ("eigen-analysis: type");
00324     TEUCHOS_TEST_FOR_EXCEPTION(
00325       eigenAnalysisType != "power-method" && 
00326       eigenAnalysisType != "power method",
00327       std::invalid_argument,
00328       "Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-"
00329       "analysis: type\" for backwards compatibility.  However, Ifpack2 "
00330       "currently only supports the \"power-method\" option for this "
00331       "parameter.  This imitates Ifpack, which only implements the power "
00332       "method for eigenanalysis.");
00333   }
00334 
00335   userInvDiag_ = userInvDiag;
00336   lambdaMax_ = lambdaMax;
00337   lambdaMin_ = lambdaMin;
00338   eigRatio_ = eigRatio;
00339   minDiagVal_ = minDiagVal;
00340   numIters_ = numIters;
00341   eigMaxIters_ = eigMaxIters;
00342   zeroStartingSolution_ = zeroStartingSolution;
00343   textbookAlgorithm_ = textbookAlgorithm;
00344   computeMaxResNorm_ = computeMaxResNorm;
00345 }
00346 
00347 template<class ScalarType, class MV, class MAT>
00348 void
00349 Chebyshev<ScalarType, MV, MAT>::
00350 compute () {
00351 #ifdef HAVE_TEUCHOS_DEBUG
00352 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00353   {
00354     // mfh 30 Jan 2013: Make sure that A is not bogus.
00355     const MT frobNormA = A_->getFrobeniusNorm ();
00356     TEUCHOS_TEST_FOR_EXCEPTION(
00357       Teuchos::ScalarTraits<MT>::isnaninf (frobNormA),
00358       std::runtime_error,
00359       "Ifpack2::Details::Chebyshev::compute: Frobenius norm of A is " 
00360       << frobNormA << ", which is Inf or NaN.");
00361     V X (A_->getDomainMap ());
00362     V Y (A_->getRangeMap ());
00363     X.putScalar (STS::one ());
00364     A_->apply (X, Y);
00365     const MT infNorm = Y.normInf ();
00366     TEUCHOS_TEST_FOR_EXCEPTION(
00367       Teuchos::ScalarTraits<MT>::isnaninf (infNorm),
00368       std::runtime_error,
00369       "Ifpack2::Details::Chebyshev: inf-norm of A*[1,...,1]^T is " 
00370       << infNorm << ", which is Inf or NaN.");
00371   }
00372 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00373 #endif // HAVE_TEUCHOS_DEBUG
00374 
00375   if (userInvDiag_.is_null ()) {
00376     // The matrix may have changed, and the parameters also affect
00377     // this process, so recompute the inverse diagonal.
00378     D_ = makeInverseDiagonal (*A_);
00379   } else {
00380     D_ = userInvDiag_;
00381   }
00382 
00383 #ifdef HAVE_TEUCHOS_DEBUG
00384 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00385   {
00386     // mfh 30 Jan 2013: Make sure that D_ contains no Inf or Nan entries.
00387     const MT infNormD = D_->normInf ();
00388     TEUCHOS_TEST_FOR_EXCEPTION(
00389       Teuchos::ScalarTraits<MT>::isnaninf (infNormD),
00390       std::runtime_error,
00391       "Ifpack2::Details::Chebyshev::compute: Infinity norm of the "
00392       "vector of inverse diagonal entries of A is " << infNormD 
00393       << ", which is Inf or NaN.");
00394   }
00395 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00396 #endif // HAVE_TEUCHOS_DEBUG
00397 
00398   // The matrix may have changed, so we always compute the
00399   // eigenvalue estimates, even if we computed them before.
00400   // However, if the user gave us lambdaMax already (if it's not
00401   // NaN), then don't (re)compute it.
00402   if (STS::isnaninf (lambdaMax_)) {
00403     const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
00404     const ST computedLambdaMin = computedLambdaMax / eigRatio_;
00405 
00406     // Defer "committing" results until all computations succeeded.
00407     computedLambdaMax_ = computedLambdaMax;
00408     computedLambdaMin_ = computedLambdaMin;
00409   }
00410 }
00411 
00412 template<class ScalarType, class MV, class MAT>
00413 typename Chebyshev<ScalarType, MV, MAT>::MT
00414 Chebyshev<ScalarType, MV, MAT>::
00415 apply (const MV& B, MV& X) {
00416   ST eigRatio = eigRatio_;
00417   ST lambdaMax = lambdaMax_;
00418   ST lambdaMin = lambdaMax / eigRatio_; // lambdaMin_; (???)
00419 
00420   // FIXME (mfh 22 Jan 2013) We really only want to check if it's
00421   // NaN.  Inf means something different.  However,
00422   // Teuchos::ScalarTraits doesn't distinguish the two cases.
00423   if (STS::isnaninf (lambdaMax)) {
00424     lambdaMax = computedLambdaMax_;
00425   }
00426   TEUCHOS_TEST_FOR_EXCEPTION(
00427     STS::isnaninf (lambdaMax), 
00428     std::runtime_error, 
00429     "Chebyshev::apply: Both lambdaMax_ and computedLambdaMax_ are NaN or Inf.  "
00430     "This means one of the following: (a) you didn't call compute() before "
00431     "calling apply(), "
00432     "(b) you didn't call compute() after calling setParameters(), or "
00433     "(c) there is a bug in this class (compute() has not done the "
00434     "eigenanalysis that it should have done).");
00435   if (STS::isnaninf (lambdaMin)) {
00436     lambdaMin = computedLambdaMin_;
00437   }
00438   TEUCHOS_TEST_FOR_EXCEPTION(
00439     STS::isnaninf (lambdaMin), 
00440     std::runtime_error, 
00441     "Chebyshev::apply: Both lambdaMin_ and computedLambdaMin_ are NaN or Inf.  "
00442     "This means one of the following: (a) you didn't call compute() before "
00443     "calling apply(), "
00444     "(b) you didn't call compute() after calling setParameters(), or "
00445     "(c) there is a bug in this class (compute() has not done the "
00446     "eigenanalysis that it should have done).");
00447   TEUCHOS_TEST_FOR_EXCEPTION(
00448     STS::isnaninf (eigRatio), 
00449     std::logic_error, 
00450     "Chebyshev::apply: eigRatio is NaN or Inf.");
00451   TEUCHOS_TEST_FOR_EXCEPTION(
00452     D_.is_null (),
00453     std::runtime_error, 
00454     "Chebyshev::apply: The vector of inverse diagonal entries of the matrix is null.  "
00455     "This means either that you didn't call compute() before calling apply(), "
00456     "or that you didn't call compute() after calling setParameters().");
00457 
00458 #ifdef HAVE_TEUCHOS_DEBUG
00459 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00460   {
00461     const MT infNormD = D_->normInf ();
00462     TEUCHOS_TEST_FOR_EXCEPTION(
00463       STS::isnaninf (infNormD),
00464       std::runtime_error,
00465       "Chebyshev::apply: The vector of inverse diagonal entries of the matrix "
00466       "has infinity norm " << infNormD << ", which is Inf or NaN.  This will "
00467       "make left-scaled Chebyshev iteration produce invalid results.");
00468   }
00469 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00470 #endif // HAVE_TEUCHOS_DEBUG
00471 
00472   if (! textbookAlgorithm_) {
00473     const ST one = Teuchos::as<ST> (1);
00474     if (STS::magnitude (lambdaMax - one) < Teuchos::as<MT> (1.0e-6)) {
00475       lambdaMin = one;
00476       lambdaMax = lambdaMin;
00477       eigRatio = one; // Ifpack doesn't include this line.
00478     }
00479     ifpackApplyImpl (*A_, B, X, numIters_, lambdaMax, lambdaMin, eigRatio, *D_);
00480   } else {
00481     textbookApplyImpl (*A_, B, X, numIters_, lambdaMax, lambdaMin, eigRatio, *D_);
00482   }
00483 
00484   if (computeMaxResNorm_) {
00485     MV R (B.getMap (), B.getNumVectors ());
00486     computeResidual (R, B, *A_, X);
00487     Teuchos::Array<MT> norms (B.getNumVectors ());
00488     R.norm2 (norms ());
00489     return *std::max_element (norms.begin (), norms.end ());
00490   } else {
00491     return Teuchos::as<MT> (0);
00492   }
00493 }
00494 
00495 template<class ScalarType, class MV, class MAT>
00496 void
00497 Chebyshev<ScalarType, MV, MAT>::
00498 print (std::ostream& out) {
00499   using std::endl;
00500   // YAML 1.2 (see http://www.yaml.org/spec/1.2/spec.html) block
00501   // mapping (see Section 8.2.2).  It should parse as follows:
00502   // { 
00503   //   Chebyshev: { 
00504   //     computedLambdaMax_: <value>,
00505   //     computedLambdaMin_: <value>,
00506   //     ...
00507   //     computeMaxResNorm_: <value>
00508   //   }
00509   // }
00510   // The spec says that the first line ("Chebyshev:") is limited to
00511   // 1024 characters, because this particular format requires
00512   // lookahead.
00513   out << "Chebyshev:" << endl
00514       << "  computedLambdaMax_: " << computedLambdaMax_ << endl
00515       << "  computedLambdaMin_: " << computedLambdaMin_ << endl
00516       << "  lambdaMax_: " << lambdaMax_ << endl
00517       << "  lambdaMin_: " << lambdaMin_ << endl
00518       << "  eigRatio_: " << eigRatio_ << endl
00519       << "  numIters_: " << numIters_ << endl
00520       << "  eigMaxIters_: " << eigMaxIters_ << endl
00521       << "  zeroStartingSolution_: " << zeroStartingSolution_ << endl
00522       << "  textbookAlgorithm_: " << textbookAlgorithm_ << endl
00523       << "  computeMaxResNorm_: " << computeMaxResNorm_ << endl;
00524 }
00525 
00526 template<class ScalarType, class MV, class MAT>
00527 void
00528 Chebyshev<ScalarType, MV, MAT>::
00529 computeResidual (MV& R, const MV& B, const MAT& A, const MV& X, 
00530      const Teuchos::ETransp mode) 
00531 {
00532 #ifdef HAVE_TEUCHOS_DEBUG
00533 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00534   using std::cerr;
00535   using std::endl;
00536   cerr << "- computeResidual:" << endl
00537        << "--- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl
00538        << "--- \\|B\\|_{\\infty} = " << maxNormInf (X) << endl;
00539 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00540 #endif // HAVE_TEUCHOS_DEBUG
00541 
00542   R = B;
00543   A.apply (X, R, mode, -STS::one(), STS::one());
00544 
00545 #ifdef HAVE_TEUCHOS_DEBUG
00546 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00547   cerr << "--- \\|B - A*X\\|_{\\infty} = " << maxNormInf (R) << endl;
00548 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00549 #endif // HAVE_TEUCHOS_DEBUG
00550 }
00551 
00552 template<class ScalarType, class MV, class MAT>
00553 void
00554 Chebyshev<ScalarType, MV, MAT>::
00555 solve (MV& Z, const V& D_inv, const MV& R) {
00556   Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
00557 }
00558 
00559 template<class ScalarType, class MV, class MAT>
00560 void
00561 Chebyshev<ScalarType, MV, MAT>::
00562 solve (MV& Z, const ST alpha, const V& D_inv, const MV& R) {
00563   Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
00564 }
00565 
00566 template<class ScalarType, class MV, class MAT>
00567 Teuchos::RCP<typename Chebyshev<ScalarType, MV, MAT>::V>
00568 Chebyshev<ScalarType, MV, MAT>::
00569 makeInverseDiagonal (const MAT& A) const {
00570   using Teuchos::RCP;
00571 
00572   RCP<V> D_rowMap (new V (A.getGraph ()->getRowMap ()));
00573   A.getLocalDiagCopy (*D_rowMap);
00574   RCP<V> D_rangeMap = makeRangeMapVector (*D_rowMap);
00575 
00576   // Invert the diagonal entries, replacing entries less than
00577   // machine precision with machine precision.
00578   typedef Kokkos::MultiVector<ST, typename MAT::node_type> KMV;
00579   KMV& localDiag = D_rangeMap->getLocalMVNonConst ();
00580   typedef Kokkos::DefaultArithmetic<KMV> KMVT;
00581   KMVT::ReciprocalThreshold (localDiag, minDiagVal_);
00582   return D_rangeMap;
00583 }
00584 
00585 template<class ScalarType, class MV, class MAT>
00586 Teuchos::RCP<typename Chebyshev<ScalarType, MV, MAT>::V>
00587 Chebyshev<ScalarType, MV, MAT>::
00588 makeRangeMapVector (const V& D) const {
00589   using Teuchos::RCP;
00590   using Teuchos::rcp;
00591   typedef Tpetra::Export<typename MV::local_ordinal_type,
00592        typename MV::global_ordinal_type,
00593        typename MV::node_type> export_type;
00594   RCP<const map_type> sourceMap = D.getMap ();
00595   RCP<const map_type> rangeMap = A_->getRangeMap ();
00596   RCP<const map_type> rowMap = A_->getRowMap ();
00597   RCP<V> D_out;
00598 
00599   if (rangeMap.getRawPtr () == sourceMap.getRawPtr ()) {
00600     // The given vector's Map is identical to the matrix's range Map.
00601     // That means we don't need to Export.  (We don't call isSameAs()
00602     // here, because that takes at least one global reduction.  This
00603     // is the fast path.  We may call isSameAs() in the slow path
00604     // below, to avoid an even slower path.
00605     D_out = rcp (new V (D));
00606   }
00607   else { // We need to Export.
00608     RCP<const export_type> exporter;
00609     // Making an Export object from scratch is expensive enough that
00610     // it's worth the O(1) global reductions to call isSameAs(), to
00611     // see if we can avoid that cost.
00612     if (sourceMap.getRawPtr () == rowMap.getRawPtr () || 
00613   sourceMap->isSameAs (*rowMap)) {
00614       // We can reuse the matrix's Export object, if there is one.
00615       exporter = A_->getGraph ()->getExporter ();
00616     }
00617     else { // We have to make a new Export object.
00618       exporter = rcp (new export_type (sourceMap, rangeMap));
00619     }
00620      
00621     D_out = rcp (new V (D));
00622     if (exporter.is_null ()) { 
00623       // Row Map and range Map are the same; no need to Export.
00624       *D_out = D;
00625     }
00626     else {
00627       D_out->doExport (D, *exporter, Tpetra::ADD);
00628     }
00629   } // if we don't need to Export, or if we do
00630 
00631   return D_out;
00632 }
00633 
00634 
00635 template<class ScalarType, class MV, class MAT>
00636 void
00637 Chebyshev<ScalarType, MV, MAT>::
00638 textbookApplyImpl (const MAT& A,
00639        const MV& B,
00640        MV& X,
00641        const int numIters,
00642        const ST lambdaMax,
00643        const ST lambdaMin,
00644        const ST eigRatio,
00645        const V& D_inv) const
00646 {
00647   (void) lambdaMin; // Forestall compiler warning.
00648   const ST myLambdaMin = lambdaMax / eigRatio;
00649 
00650   const ST zero = Teuchos::as<ST> (0);
00651   const ST one = Teuchos::as<ST> (1);
00652   const ST two = Teuchos::as<ST> (2);
00653   const ST d = (lambdaMax + myLambdaMin) / two; // Ifpack2 calls this theta
00654   const ST c = (lambdaMax - myLambdaMin) / two; // Ifpack2 calls this 1/delta
00655 
00656   if (zeroStartingSolution_ && numIters > 0) {
00657     // If zero iterations, then input X is output X.
00658     X.putScalar (zero);
00659   }
00660   MV R (B.getMap (), B.getNumVectors (), false);
00661   MV P (B.getMap (), B.getNumVectors (), false);
00662   MV Z (B.getMap (), B.getNumVectors (), false);
00663   ST alpha, beta;
00664   for (int i = 0; i < numIters; ++i) {
00665     computeResidual (R, B, A, X); // R = B - A*X
00666     solve (Z, D_inv, R); // z = D_inv * R, that is, D \ R.
00667     if (i == 0) {
00668       P = Z;
00669       alpha = two / d;
00670     } else {
00671       //beta = (c * alpha / two)^2;
00672       //const ST sqrtBeta = c * alpha / two;
00673       //beta = sqrtBeta * sqrtBeta;
00674       beta = alpha * (c/two) * (c/two);
00675       alpha = one / (d - beta);
00676       P.update (one, Z, beta); // P = Z + beta*P
00677     }
00678     X.update (alpha, P, one); // X = X + alpha*P
00679     // If we compute the residual here, we could either do R = B -
00680     // A*X, or R = R - alpha*A*P.  Since we choose the former, we
00681     // can move the computeResidual call to the top of the loop.
00682   }
00683 }
00684 
00685 template<class ScalarType, class MV, class MAT>
00686 typename Chebyshev<ScalarType, MV, MAT>::MT
00687 Chebyshev<ScalarType, MV, MAT>::maxNormInf (const MV& X) {
00688   std::vector<MT> norms (X.getNumVectors ());
00689   X.normInf (norms);
00690   return *std::max_element (norms.begin (), norms.end ());
00691 }
00692 
00693 template<class ScalarType, class MV, class MAT>
00694 void
00695 Chebyshev<ScalarType, MV, MAT>::
00696 ifpackApplyImpl (const MAT& A,
00697      const MV& B,
00698      MV& X,
00699      const int numIters,
00700      const ST lambdaMax,
00701      const ST lambdaMin,
00702      const ST eigRatio,
00703      const V& D_inv) 
00704 {
00705 #ifdef HAVE_TEUCHOS_DEBUG
00706 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00707   using std::cerr;
00708   using std::endl;
00709   cerr << "\\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
00710   cerr << "\\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
00711 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00712 #endif // HAVE_TEUCHOS_DEBUG
00713 
00714   if (numIters <= 0) {
00715     return;
00716   }
00717   const ST one = Teuchos::as<ST> (1);
00718   const ST two = Teuchos::as<ST> (2);
00719 
00720   // Quick solve when the matrix A is the identity.
00721   if (lambdaMin == one && lambdaMax == lambdaMin) {
00722     solve (X, D_inv, B);
00723     return;
00724   }
00725 
00726   // Initialize coefficients
00727   const ST alpha = lambdaMax / eigRatio;
00728   const ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
00729   const ST delta = two / (beta - alpha);
00730   const ST theta = (beta + alpha) / two;
00731   const ST s1 = theta * delta;
00732 
00733 #ifdef HAVE_TEUCHOS_DEBUG
00734 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00735   cerr << "alpha = " << alpha << endl
00736        << "beta = " << beta << endl
00737        << "delta = " << delta << endl
00738        << "theta = " << theta << endl
00739        << "s1 = " << s1 << endl;
00740 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00741 #endif // HAVE_TEUCHOS_DEBUG
00742 
00743   // Fetch cached temporary vectors.
00744   Teuchos::RCP<MV> V_ptr, W_ptr;
00745   makeTempMultiVectors (V_ptr, W_ptr, B);
00746 
00747   // mfh 28 Jan 2013: We write V1 instead of V, so as not to confuse
00748   // the multivector V with the typedef V (for Tpetra::Vector).
00749   //MV V1 (B.getMap (), B.getNumVectors (), false);
00750   //MV W (B.getMap (), B.getNumVectors (), false);
00751   MV& V1 = *V_ptr;
00752   MV& W = *W_ptr;
00753 
00754 #ifdef HAVE_TEUCHOS_DEBUG
00755 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00756   cerr << "Iteration " << 1 << ":" << endl
00757        << "- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
00758 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00759 #endif // HAVE_TEUCHOS_DEBUG
00760 
00761   // Special case for the first iteration.
00762   if (! zeroStartingSolution_) {
00763     computeResidual (V1, B, A, X); // V1 = B - A*X
00764 
00765 #ifdef HAVE_TEUCHOS_DEBUG
00766 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00767     cerr << "- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
00768 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00769 #endif // HAVE_TEUCHOS_DEBUG
00770 
00771     solve (W, one/theta, D_inv, V1); // W = (1/theta)*D_inv*(B-A*X)
00772 
00773 #ifdef HAVE_TEUCHOS_DEBUG
00774 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00775     cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
00776 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00777 #endif // HAVE_TEUCHOS_DEBUG
00778 
00779     X.update (one, W, one); // X = X + W
00780   } else {
00781     solve (W, one/theta, D_inv, B); // W = (1/theta)*D_inv*B
00782 
00783 #ifdef HAVE_TEUCHOS_DEBUG
00784 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00785     cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
00786 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00787 #endif // HAVE_TEUCHOS_DEBUG
00788 
00789     X = W; // X = 0 + W
00790   }
00791 #ifdef HAVE_TEUCHOS_DEBUG
00792 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00793   cerr << "- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
00794 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00795 #endif // HAVE_TEUCHOS_DEBUG
00796 
00797   // The rest of the iterations.
00798   ST rhok = one / s1;
00799   ST rhokp1, dtemp1, dtemp2;
00800   for (int deg = 1; deg < numIters; ++deg) {
00801 
00802 #ifdef HAVE_TEUCHOS_DEBUG
00803 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00804     cerr << "Iteration " << deg+1 << ":" << endl;
00805     cerr << "- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
00806     cerr << "- \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
00807     cerr << "- \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm () << endl;
00808     cerr << "- rhok = " << rhok << endl;
00809     V1.putScalar (STS::zero ());
00810 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00811 #endif // HAVE_TEUCHOS_DEBUG
00812 
00813     computeResidual (V1, B, A, X); // V1 = B - A*X
00814 
00815 #ifdef HAVE_TEUCHOS_DEBUG
00816 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00817     cerr << "- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
00818 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00819 #endif // HAVE_TEUCHOS_DEBUG
00820 
00821     rhokp1 = one / (two * s1 - rhok);
00822     dtemp1 = rhokp1 * rhok;
00823     dtemp2 = two * rhokp1 * delta;
00824     rhok = rhokp1;
00825 
00826 #ifdef HAVE_TEUCHOS_DEBUG
00827 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00828     cerr << "- dtemp1 = " << dtemp1 << endl
00829    << "- dtemp2 = " << dtemp1 << endl;
00830 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00831 #endif // HAVE_TEUCHOS_DEBUG
00832 
00833     W.scale (dtemp1);
00834     W.elementWiseMultiply (dtemp2, D_inv, V1, one);
00835     X.update (one, W, one);
00836 
00837 #ifdef HAVE_TEUCHOS_DEBUG
00838 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00839     cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
00840     cerr << "- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
00841 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00842 #endif // HAVE_TEUCHOS_DEBUG
00843   }
00844 }
00845 
00846 template<class ScalarType, class MV, class MAT>
00847 typename Chebyshev<ScalarType, MV, MAT>::ST
00848 Chebyshev<ScalarType, MV, MAT>::
00849 powerMethod (const MAT& A, const V& D_inv, const int numIters) {
00850   const ST zero = Teuchos::as<ST> (0);
00851   const ST one = Teuchos::as<ST> (1);
00852   ST lambdaMax = zero;
00853   ST RQ_top, RQ_bottom, norm;
00854 
00855   V x (A.getDomainMap ());
00856   V y (A.getRangeMap ());
00857   x.randomize ();
00858   norm = x.norm2 ();
00859   TEUCHOS_TEST_FOR_EXCEPTION(norm == zero, std::runtime_error, 
00860     "Chebyshev::powerMethod: Tpetra::Vector::randomize filled the vector "
00861     "with zeros.  This is not impossible, but is unlikely.");
00862 
00863   x.scale (one / norm);
00864   for (int iter = 0; iter < numIters; ++iter) {
00865     A.apply (x, y);
00866     solve (y, D_inv, y);
00867     RQ_top = y.dot (x);
00868     RQ_bottom = x.dot (x);
00869     lambdaMax = RQ_top / RQ_bottom;
00870     norm = y.norm2 ();
00871     if (norm == zero) { // Return something reasonable.
00872       return zero;
00873     }
00874     x.update (one / norm, y, zero);
00875   }
00876   return lambdaMax;
00877 }
00878 
00879 template<class ScalarType, class MV, class MAT>
00880 Teuchos::RCP<const MAT>
00881 Chebyshev<ScalarType, MV, MAT>::
00882 getMatrix() const {
00883   return A_;
00884 }
00885 
00886 template<class ScalarType, class MV, class MAT>
00887 bool
00888 Chebyshev<ScalarType, MV, MAT>::
00889 hasTransposeApply() const {
00890   // Technically, this is true, because the matrix must be symmetric.
00891   return true;
00892 }
00893 
00894 template<class ScalarType, class MV, class MAT>
00895 void
00896 Chebyshev<ScalarType, MV, MAT>::
00897 makeTempMultiVectors (Teuchos::RCP<MV>& V1,
00898           Teuchos::RCP<MV>& W,
00899           const MV& X)
00900 {
00901   // Don't fill the vectors with zeros.
00902   if (V_.is_null ()) {
00903     V_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), false));
00904   }
00905   if (W_.is_null ()) {
00906     W_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), false));
00907   }
00908   V1 = V_;
00909   W = W_;
00910 }
00911 
00912 } // namespace Details
00913 } // namespace Ifpack2
00914 
00915 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends