Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Details_Chebyshev_def.hpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00006 //                 Copyright (2009) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 // ***********************************************************************
00045 //
00046 //      Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
00047 //                 Copyright (2004) Sandia Corporation
00048 //
00049 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00050 // license for use of this work by or on behalf of the U.S. Government.
00051 //
00052 // This library is free software; you can redistribute it and/or modify
00053 // it under the terms of the GNU Lesser General Public License as
00054 // published by the Free Software Foundation; either version 2.1 of the
00055 // License, or (at your option) any later version.
00056 //
00057 // This library is distributed in the hope that it will be useful, but
00058 // WITHOUT ANY WARRANTY; without even the implied warranty of
00059 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00060 // Lesser General Public License for more details.
00061 //
00062 // You should have received a copy of the GNU Lesser General Public
00063 // License along with this library; if not, write to the Free Software
00064 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00065 // USA
00066 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00067 //
00068 // ***********************************************************************
00069 
00070 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
00071 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
00072 
00079 
00080 #include "Ifpack2_Details_Chebyshev_decl.hpp"
00081 #include <cmath>
00082 
00083 // Uncommit the #define line below if you want Chebyshev to do extra
00084 // debug checking and generate a lot of debugging output to stderr (on
00085 // all processes in the communicator).  Even if you uncomment this
00086 // line, the debugging code will only be enabled if the CMake option
00087 // Teuchos_ENABLE_DEBUG was set to ON when configuring Trilinos.
00088 //#define IFPACK_DETAILS_CHEBYSHEV_DEBUG 1
00089 
00090 namespace Ifpack2 {
00091 namespace Details {
00092 
00093 namespace {
00094   // We use this text a lot in error messages.
00095   const char computeBeforeApplyReminder[] =
00096     "This means one of the following:\n"
00097     "  - you have not yet called compute() on this instance, or \n"
00098     "  - you didn't call compute() after calling setParameters().\n\n"
00099     "After creating an Ifpack2::Chebyshev instance,\n"
00100     "you must _always_ call compute() at least once before calling apply().\n"
00101     "After calling compute() once, you do not need to call it again,\n"
00102     "unless the matrix has changed or you have changed parameters\n"
00103     "(by calling setParameters()).";
00104 }
00105 
00106 template<class ScalarType, class MV>
00107 void Chebyshev<ScalarType, MV>::checkInputMatrix () const
00108 {
00109   TEUCHOS_TEST_FOR_EXCEPTION(
00110     ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
00111     std::invalid_argument,
00112     "Ifpack2::Chebyshev: The input matrix A must be square.  "
00113     "A has " << A_->getGlobalNumRows () << " rows and "
00114     << A_->getGlobalNumCols () << " columns.");
00115 
00116   // In a debug build, test that the domain and range Maps of the
00117   // matrix are the same.
00118 #ifdef HAVE_TEUCHOS_DEBUG
00119   if (! A_.is_null ()) {
00120     Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
00121     Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
00122 
00123     // isSameAs is a collective, but if the two pointers are the same,
00124     // isSameAs will assume that they are the same on all processes, and
00125     // return true without an all-reduce.
00126     TEUCHOS_TEST_FOR_EXCEPTION(
00127       ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
00128       "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
00129       "the same (in the sense of isSameAs())." << std::endl << "We only check "
00130       "for this if Trilinos was built with the CMake configuration option "
00131       "Teuchos_ENABLE_DEBUG set to ON.");
00132   }
00133 #endif // HAVE_TEUCHOS_DEBUG
00134 }
00135 
00136 
00137 template<class ScalarType, class MV>
00138 void
00139 Chebyshev<ScalarType, MV>::
00140 checkConstructorInput () const
00141 {
00142   TEUCHOS_TEST_FOR_EXCEPTION(STS::isComplex, std::logic_error,
00143     "Ifpack2::Chebyshev: This class' implementation of Chebyshev iteration "
00144     "only works for real-valued, symmetric positive definite matrices.  "
00145     "However, you instantiated this class for ScalarType = "
00146     << Teuchos::TypeNameTraits<ScalarType>::name () << ", which is a complex-"
00147     "valued type.  While this may be algorithmically correct if all of the "
00148     "complex numbers in the matrix have zero imaginary part, we forbid using "
00149     "complex ScalarType altogether in order to remind you of the limitations "
00150     "of our implementation (and of the algorithm itself).");
00151 
00152   checkInputMatrix ();
00153 }
00154 
00155 template<class ScalarType, class MV>
00156 Chebyshev<ScalarType, MV>::
00157 Chebyshev (Teuchos::RCP<const row_matrix_type> A) :
00158   A_ (A),
00159   savedDiagOffsets_ (false),
00160   computedLambdaMax_ (STS::nan ()),
00161   computedLambdaMin_ (STS::nan ()),
00162   lambdaMaxForApply_ (STS::nan ()),
00163   lambdaMinForApply_ (STS::nan ()),
00164   eigRatioForApply_ (STS::nan ()),
00165   userLambdaMax_ (STS::nan ()),
00166   userLambdaMin_ (STS::nan ()),
00167   userEigRatio_ (Teuchos::as<ST> (30)),
00168   minDiagVal_ (STS::eps ()),
00169   numIters_ (1),
00170   eigMaxIters_ (10),
00171   zeroStartingSolution_ (true),
00172   assumeMatrixUnchanged_ (false),
00173   textbookAlgorithm_ (false),
00174   computeMaxResNorm_ (false)
00175 {
00176   checkConstructorInput ();
00177 }
00178 
00179 template<class ScalarType, class MV>
00180 Chebyshev<ScalarType, MV>::
00181 Chebyshev (Teuchos::RCP<const row_matrix_type> A, Teuchos::ParameterList& params) :
00182   A_ (A),
00183   savedDiagOffsets_ (false),
00184   computedLambdaMax_ (STS::nan ()),
00185   computedLambdaMin_ (STS::nan ()),
00186   lambdaMaxForApply_ (STS::nan ()),
00187   lambdaMinForApply_ (STS::nan ()),
00188   eigRatioForApply_ (STS::nan ()),
00189   userLambdaMax_ (STS::nan ()),
00190   userLambdaMin_ (STS::nan ()),
00191   userEigRatio_ (Teuchos::as<ST> (30)),
00192   minDiagVal_ (STS::eps ()),
00193   numIters_ (1),
00194   eigMaxIters_ (10),
00195   zeroStartingSolution_ (true),
00196   assumeMatrixUnchanged_ (false),
00197   textbookAlgorithm_ (false),
00198   computeMaxResNorm_ (false)
00199 {
00200   checkConstructorInput ();
00201   setParameters (params);
00202 }
00203 
00204 template<class ScalarType, class MV>
00205 void
00206 Chebyshev<ScalarType, MV>::
00207 setParameters (Teuchos::ParameterList& plist)
00208 {
00209   using Teuchos::RCP;
00210   using Teuchos::rcp;
00211   using Teuchos::rcp_const_cast;
00212   using Teuchos::rcpFromRef;
00213 
00214   // Note to developers: The logic for this method is complicated,
00215   // because we want to accept Ifpack and ML parameters whenever
00216   // possible, but we don't want to add their default values to the
00217   // user's ParameterList.  That's why we do all the isParameter()
00218   // checks, instead of using the two-argument version of get()
00219   // everywhere.  The min and max eigenvalue parameters are also a
00220   // special case, because we decide whether or not to do eigenvalue
00221   // analysis based on whether the user supplied the max eigenvalue.
00222 
00223   // Default values of all the parameters.
00224   const ST defaultLambdaMax = STS::nan ();
00225   const ST defaultLambdaMin = STS::nan ();
00226   // 30 is Ifpack::Chebyshev's default.  ML has a different default
00227   // eigRatio for smoothers and the coarse-grid solve (if using
00228   // Chebyshev for that).  The former uses 20; the latter uses 30.
00229   // We're testing smoothers here, so use 20.  (However, if you give
00230   // ML an Epetra matrix, it will use Ifpack for Chebyshev, in which
00231   // case it would defer to Ifpack's default settings.)
00232   const ST defaultEigRatio = Teuchos::as<ST> (30);
00233   const ST defaultMinDiagVal = STS::eps ();
00234   const int defaultNumIters = 1;
00235   const int defaultEigMaxIters = 10;
00236   const bool defaultZeroStartingSolution = true; // Ifpack::Chebyshev default
00237   const bool defaultAssumeMatrixUnchanged = false;
00238   const bool defaultTextbookAlgorithm = false;
00239   const bool defaultComputeMaxResNorm = false;
00240 
00241   // We'll set the instance data transactionally, after all reads
00242   // from the ParameterList.  That way, if any of the ParameterList
00243   // reads fail (e.g., due to the wrong parameter type), we will not
00244   // have left the instance data in a half-changed state.
00245   RCP<const V> userInvDiag;
00246   ST lambdaMax = defaultLambdaMax;
00247   ST lambdaMin = defaultLambdaMin;
00248   ST eigRatio = defaultEigRatio;
00249   ST minDiagVal = defaultMinDiagVal;
00250   int numIters = defaultNumIters;
00251   int eigMaxIters = defaultEigMaxIters;
00252   bool zeroStartingSolution = defaultZeroStartingSolution;
00253   bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
00254   bool textbookAlgorithm = defaultTextbookAlgorithm;
00255   bool computeMaxResNorm = defaultComputeMaxResNorm;
00256 
00257   // Fetch the parameters from the ParameterList.  Defer all
00258   // externally visible side effects until we have finished all
00259   // ParameterList interaction.  This makes the method satisfy the
00260   // strong exception guarantee.
00261 
00262   // Get the user-supplied inverse diagonal.
00263   //
00264   // Check for a raw pointer (const V* or V*), for Ifpack
00265   // compatibility, as well as for RCP<const V>, RCP<V>, const V, or
00266   // V.  We'll copy the vector anyway, so it doesn't matter whether
00267   // it's const or nonconst.
00268   if (plist.isParameter ("chebyshev: operator inv diagonal")) {
00269     try { // Could the type be const V*?
00270       const V* rawUserInvDiag =
00271         plist.get<const V*> ("chebyshev: operator inv diagonal");
00272       // If it's a raw pointer, then we have to copy it, since we
00273       // can't otherwise ensure that the user won't deallocate the
00274       // Vector before we use it.
00275       userInvDiag = rcp (new V (*rawUserInvDiag));
00276     } catch (Teuchos::Exceptions::InvalidParameterType&) {
00277     }
00278     if (userInvDiag.is_null ()) {
00279       try { // Could the type be V*?
00280         V* rawUserInvDiag = plist.get<V*> ("chebyshev: operator inv diagonal");
00281         // If it's a raw pointer, then we have to copy it, since we
00282         // can't otherwise ensure that the user won't deallocate the
00283         // Vector before we use it.
00284         userInvDiag = rcp (new V (*rawUserInvDiag));
00285       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00286       }
00287     }
00288     if (userInvDiag.is_null ()) {
00289       try { // Could the type be RCP<const V>?
00290         // If the type is RCP<const V>, then the user has promised
00291         // that they won't change the Vector.  Thus, it's safe just to
00292         // keep the RCP; we don't have to make a copy.
00293         userInvDiag =
00294           plist.get<RCP<const V> > ("chebyshev: operator inv diagonal");
00295       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00296       }
00297     }
00298     if (userInvDiag.is_null ()) {
00299       try { // Could the type be RCP<V>?
00300         RCP<V> userInvDiagNonconst =
00301           plist.get<RCP<V> > ("chebyshev: operator inv diagonal");
00302         // If the type is RCP<V>, the user could still change the
00303         // Vector.  That means we have to make a deep copy.
00304         userInvDiag = rcp (new V (*userInvDiagNonconst));
00305       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00306       }
00307     }
00308     if (userInvDiag.is_null ()) {
00309       try { // Could the type be const V?
00310         // The line below does a deep copy (V::operator=).
00311         userInvDiag =
00312           rcp (new V (plist.get<const V> ("chebyshev: operator inv diagonal")));
00313       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00314       }
00315     }
00316     if (userInvDiag.is_null ()) {
00317       try { // Could the type be V?
00318         // The line below does a deep copy (V::operator=).
00319         userInvDiag =
00320           rcp (new V (plist.get<V> ("chebyshev: operator inv diagonal")));
00321       } catch (Teuchos::Exceptions::InvalidParameterType&) {
00322       }
00323     }
00324     // We don't necessarily have a range Map yet.  compute() is the
00325     // proper place to compute the range Map version of userInvDiag.
00326   }
00327 
00328   // Don't fill in defaults for the max or min eigenvalue, because
00329   // this class uses the existence of those parameters to determine
00330   // whether it should do eigenanalysis.
00331   if (plist.isParameter ("chebyshev: max eigenvalue")) {
00332     if (plist.isType<double>("chebyshev: max eigenvalue"))
00333       lambdaMax = plist.get<double> ("chebyshev: max eigenvalue");
00334     else
00335       lambdaMax = plist.get<ST> ("chebyshev: max eigenvalue");
00336     TEUCHOS_TEST_FOR_EXCEPTION(
00337       STS::isnaninf (lambdaMax), std::invalid_argument,
00338       "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
00339       "parameter is NaN or Inf.  This parameter is optional, but if you "
00340       "choose to supply it, it must have a finite value.");
00341   }
00342   if (plist.isParameter ("chebyshev: min eigenvalue")) {
00343     if (plist.isType<double>("chebyshev: min eigenvalue"))
00344       lambdaMin = plist.get<double> ("chebyshev: min eigenvalue");
00345     else
00346       lambdaMin = plist.get<ST> ("chebyshev: min eigenvalue");
00347     TEUCHOS_TEST_FOR_EXCEPTION(
00348       STS::isnaninf (lambdaMin), std::invalid_argument,
00349       "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
00350       "parameter is NaN or Inf.  This parameter is optional, but if you "
00351       "choose to supply it, it must have a finite value.");
00352   }
00353 
00354   // Only fill in Ifpack2's name for the default parameter, not ML's.
00355   if (plist.isParameter ("smoother: Chebyshev alpha")) { // ML compatibility
00356     if (plist.isType<double>("smoother: Chebyshev alpha"))
00357       eigRatio = plist.get<double> ("smoother: Chebyshev alpha");
00358     else
00359       eigRatio = plist.get<ST> ("smoother: Chebyshev alpha");
00360   }
00361   // Ifpack2's name overrides ML's name.
00362   eigRatio = plist.get ("chebyshev: ratio eigenvalue", eigRatio);
00363   TEUCHOS_TEST_FOR_EXCEPTION(
00364     STS::isnaninf (eigRatio), std::invalid_argument,
00365     "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
00366     "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf.  "
00367     "This parameter is optional, but if you choose to supply it, it must have "
00368     "a finite value.");
00369   // mfh 11 Feb 2013: This class is currently only correct for real
00370   // Scalar types, but we still want it to build for complex Scalar
00371   // type so that users of Ifpack2::Factory can build their
00372   // executables for real or complex Scalar type.  Thus, we take the
00373   // real parts here, which are always less-than comparable.
00374   TEUCHOS_TEST_FOR_EXCEPTION(
00375     STS::real (eigRatio) < STS::real (STS::one ()),
00376     std::invalid_argument,
00377     "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
00378     "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
00379     "but you supplied the value " << eigRatio << ".");
00380 
00381   // Same name in Ifpack2 and Ifpack.
00382   minDiagVal = plist.get ("chebyshev: min diagonal value", minDiagVal);
00383   TEUCHOS_TEST_FOR_EXCEPTION(
00384     STS::isnaninf (minDiagVal), std::invalid_argument,
00385     "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
00386     "parameter is NaN or Inf.  This parameter is optional, but if you choose "
00387     "to supply it, it must have a finite value.");
00388 
00389   // Only fill in Ifpack2's name, not ML's or Ifpack's.
00390   if (plist.isParameter ("smoother: sweeps")) { // ML compatibility
00391     numIters = plist.get<int> ("smoother: sweeps");
00392   } // Ifpack's name overrides ML's name.
00393   if (plist.isParameter ("relaxation: sweeps")) { // Ifpack compatibility
00394     numIters = plist.get<int> ("relaxation: sweeps");
00395   } // Ifpack2's name overrides Ifpack's name.
00396   numIters = plist.get ("chebyshev: degree", numIters);
00397   TEUCHOS_TEST_FOR_EXCEPTION(
00398     numIters < 0, std::invalid_argument,
00399     "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
00400     "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
00401     "nonnegative integer.  You gave a value of " << numIters << ".");
00402 
00403   // The last parameter name overrides the first.
00404   if (plist.isParameter ("eigen-analysis: iterations")) { // ML compatibility
00405     eigMaxIters = plist.get<int> ("eigen-analysis: iterations");
00406   } // Ifpack2's name overrides ML's name.
00407   eigMaxIters = plist.get ("chebyshev: eigenvalue max iterations", eigMaxIters);
00408   TEUCHOS_TEST_FOR_EXCEPTION(
00409     eigMaxIters < 0, std::invalid_argument,
00410     "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
00411     "\" parameter (also called \"eigen-analysis: iterations\") must be a "
00412     "nonnegative integer.  You gave a value of " << eigMaxIters << ".");
00413 
00414   zeroStartingSolution = plist.get ("chebyshev: zero starting solution",
00415                                     zeroStartingSolution);
00416   assumeMatrixUnchanged = plist.get ("chebyshev: assume matrix does not change",
00417                                      assumeMatrixUnchanged);
00418 
00419   // We don't want to fill these parameters in, because they shouldn't
00420   // be visible to Ifpack2::Chebyshev users.
00421   if (plist.isParameter ("chebyshev: textbook algorithm")) {
00422     textbookAlgorithm = plist.get<bool> ("chebyshev: textbook algorithm");
00423   }
00424   if (plist.isParameter ("chebyshev: compute max residual norm")) {
00425     computeMaxResNorm = plist.get<bool> ("chebyshev: compute max residual norm");
00426   }
00427 
00428   // Test for Ifpack parameters that we won't ever implement here.
00429   // Be careful to use the one-argument version of get(), since the
00430   // two-argment version adds the parameter if it's not there.
00431   TEUCHOS_TEST_FOR_EXCEPTION(
00432     plist.isParameter ("chebyshev: use block mode") &&
00433     ! plist.get<bool> ("chebyshev: use block mode"),
00434     std::invalid_argument,
00435     "Ifpack2::Chebyshev requires that if you supply the Ifpack parameter "
00436     "\"chebyshev: use block mode\", it must be set to false.  Ifpack2's "
00437     "Chebyshev does not implement Ifpack's block mode.");
00438   TEUCHOS_TEST_FOR_EXCEPTION(
00439     plist.isParameter ("chebyshev: solve normal equations") &&
00440     ! plist.get<bool> ("chebyshev: solve normal equations"),
00441     std::invalid_argument,
00442     "Ifpack2::Chebyshev does not and will never implement the Ifpack "
00443     "parameter \"chebyshev: solve normal equations\".  If you want to solve "
00444     "the normal equations, construct a Tpetra::Operator that implements "
00445     "A^* A, and use Chebyshev to solve A^* A x = A^* b.");
00446 
00447   // Test for Ifpack parameters that we haven't implemented yet.
00448   //
00449   // For now, we only check that this ML parameter, if provided, has
00450   // the one value that we support.  We consider other values "invalid
00451   // arguments" rather than "logic errors," because Ifpack does not
00452   // implement eigenanalyses other than the power method.
00453   std::string eigenAnalysisType ("power-method");
00454   if (plist.isParameter ("eigen-analysis: type")) {
00455     eigenAnalysisType = plist.get<std::string> ("eigen-analysis: type");
00456     TEUCHOS_TEST_FOR_EXCEPTION(
00457       eigenAnalysisType != "power-method" &&
00458       eigenAnalysisType != "power method",
00459       std::invalid_argument,
00460       "Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-"
00461       "analysis: type\" for backwards compatibility.  However, Ifpack2 "
00462       "currently only supports the \"power-method\" option for this "
00463       "parameter.  This imitates Ifpack, which only implements the power "
00464       "method for eigenanalysis.");
00465   }
00466 
00467   // We've validated all the parameters, so it's safe now to "commit" them.
00468   userInvDiag_ = userInvDiag;
00469   userLambdaMax_ = lambdaMax;
00470   userLambdaMin_ = lambdaMin;
00471   userEigRatio_ = eigRatio;
00472   minDiagVal_ = minDiagVal;
00473   numIters_ = numIters;
00474   eigMaxIters_ = eigMaxIters;
00475   zeroStartingSolution_ = zeroStartingSolution;
00476   assumeMatrixUnchanged_ = assumeMatrixUnchanged;
00477   textbookAlgorithm_ = textbookAlgorithm;
00478   computeMaxResNorm_ = computeMaxResNorm;
00479 }
00480 
00481 
00482 template<class ScalarType, class MV>
00483 void
00484 Chebyshev<ScalarType, MV>::reset ()
00485 {
00486   D_ = Teuchos::null;
00487   diagOffsets_ = Teuchos::null;
00488   savedDiagOffsets_ = false;
00489   V_ = Teuchos::null;
00490   W_ = Teuchos::null;
00491   computedLambdaMax_ = STS::nan ();
00492   computedLambdaMin_ = STS::nan ();
00493 }
00494 
00495 
00496 template<class ScalarType, class MV>
00497 void
00498 Chebyshev<ScalarType, MV>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
00499 {
00500   if (A.getRawPtr () != A_.getRawPtr ()) {
00501     if (! assumeMatrixUnchanged_) {
00502       reset ();
00503     }
00504     A_ = A;
00505   }
00506 }
00507 
00508 
00509 template<class ScalarType, class MV>
00510 void
00511 Chebyshev<ScalarType, MV>::compute ()
00512 {
00513   using std::endl;
00514   // Some of the optimizations below only work if A_ is a
00515   // Tpetra::CrsMatrix.  We'll make our best guess about its type
00516   // here, since we have no way to get back the original fifth
00517   // template parameter.
00518   typedef Tpetra::CrsMatrix<typename MV::scalar_type,
00519     typename MV::local_ordinal_type,
00520     typename MV::global_ordinal_type,
00521     typename MV::node_type> crs_matrix_type;
00522 
00523   TEUCHOS_TEST_FOR_EXCEPTION(
00524     A_.is_null (), std::runtime_error, "Ifpack2::Chebyshev::compute: The input "
00525     "matrix A is null.  Please call setMatrix() with a nonnull input matrix "
00526     "before calling this method.");
00527 
00528   // If A_ is a CrsMatrix and its graph is constant, we presume that
00529   // the user plans to reuse the structure of A_, but possibly change
00530   // A_'s values before each compute() call.  This is the intended use
00531   // case for caching the offsets of the diagonal entries of A_, to
00532   // speed up extraction of diagonal entries on subsequent compute()
00533   // calls.
00534 
00535   // FIXME (mfh 22 Jan 2013, 10 Feb 2013) In all cases when we use
00536   // isnaninf() in this method, we really only want to check if the
00537   // number is NaN.  Inf means something different.  However,
00538   // Teuchos::ScalarTraits doesn't distinguish the two cases.
00539 
00540   // makeInverseDiagonal() returns a range Map Vector.
00541   if (userInvDiag_.is_null ()) {
00542     Teuchos::RCP<const crs_matrix_type> A_crsMat =
00543       Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
00544 
00545     if (D_.is_null ()) { // We haven't computed D_ before
00546       if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
00547         // It's a CrsMatrix with a const graph; cache diagonal offsets.
00548         A_crsMat->getLocalDiagOffsets (diagOffsets_);
00549         savedDiagOffsets_ = true;
00550         D_ = makeInverseDiagonal (*A_, true);
00551       }
00552       else { // either A_ is not a CrsMatrix, or its graph is nonconst
00553         D_ = makeInverseDiagonal (*A_);
00554       }
00555     }
00556     else if (! assumeMatrixUnchanged_) { // D_ exists but A_ may have changed
00557       if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
00558         // It's a CrsMatrix with a const graph; cache diagonal offsets
00559         // if we haven't already.
00560         if (! savedDiagOffsets_) {
00561           A_crsMat->getLocalDiagOffsets (diagOffsets_);
00562           savedDiagOffsets_ = true;
00563         }
00564         // Now we're guaranteed to have cached diagonal offsets.
00565         D_ = makeInverseDiagonal (*A_, true);
00566       }
00567       else { // either A_ is not a CrsMatrix, or its graph is nonconst
00568         D_ = makeInverseDiagonal (*A_);
00569       }
00570     }
00571   }
00572   else { // the user provided an inverse diagonal
00573     D_ = makeRangeMapVectorConst (userInvDiag_);
00574   }
00575 
00576   // Have we estimated eigenvalues before?
00577   const bool computedEigenvalueEstimates =
00578     STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
00579 
00580   // Only recompute the eigenvalue estimates if
00581   // - we are supposed to assume that the matrix may have changed, or
00582   // - they haven't been computed before, and the user hasn't given
00583   //   us at least an estimate of the max eigenvalue.
00584   //
00585   // We at least need an estimate of the max eigenvalue.  This is the
00586   // most important one if using Chebyshev as a smoother.
00587   if (! assumeMatrixUnchanged_ ||
00588       (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
00589     const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
00590     TEUCHOS_TEST_FOR_EXCEPTION(
00591       STS::isnaninf (computedLambdaMax),
00592       std::runtime_error,
00593       "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
00594       "of D^{-1} A failed, by producing Inf or NaN.  This probably means that "
00595       "the matrix contains Inf or NaN values, or that it is badly scaled.");
00596     TEUCHOS_TEST_FOR_EXCEPTION(
00597       STS::isnaninf (userEigRatio_),
00598       std::logic_error,
00599       "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
00600       << endl << "This should be impossible." << endl <<
00601       "Please report this bug to the Ifpack2 developers.");
00602 
00603     // The power method doesn't estimate the min eigenvalue, so we
00604     // do our best to provide an estimate.  userEigRatio_ has a
00605     // reasonable default value, and if the user provided it, we
00606     // have already checked that its value is finite and >= 1.
00607     const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
00608 
00609     // Defer "committing" results until all computations succeeded.
00610     computedLambdaMax_ = computedLambdaMax;
00611     computedLambdaMin_ = computedLambdaMin;
00612   } else {
00613     TEUCHOS_TEST_FOR_EXCEPTION(
00614       STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
00615       std::logic_error,
00616       "Ifpack2::Chebyshev::compute: " << endl <<
00617       "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
00618       << endl << "This should be impossible." << endl <<
00619       "Please report this bug to the Ifpack2 developers.");
00620   }
00621 
00623   // Figure out the eigenvalue estimates that apply() will use.
00625 
00626   // Always favor the user's max eigenvalue estimate, if provided.
00627   if (STS::isnaninf (userLambdaMax_)) {
00628     lambdaMaxForApply_ = computedLambdaMax_;
00629   } else {
00630     lambdaMaxForApply_ = userLambdaMax_;
00631   }
00632   // mfh 11 Feb 2013: For now, we imitate Ifpack by ignoring the
00633   // user's min eigenvalue estimate, and using the given eigenvalue
00634   // ratio to estimate the min eigenvalue.  We could instead do this:
00635   // favor the user's eigenvalue ratio estimate, but if it's not
00636   // provided, use lambdaMax / lambdaMin.  However, we want Chebyshev
00637   // to have sensible smoother behavior if the user did not provide
00638   // eigenvalue estimates.  Ifpack's behavior attempts to push down
00639   // the error terms associated with the largest eigenvalues, while
00640   // expecting that users will only want a small number of iterations,
00641   // so that error terms associated with the smallest eigenvalues
00642   // won't grow too much.  This is sensible behavior for a smoother.
00643   lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
00644   eigRatioForApply_ = userEigRatio_;
00645 
00646   if (! textbookAlgorithm_) {
00647     // Ifpack has a special-case modification of the eigenvalue bounds
00648     // for the case where the max eigenvalue estimate is close to one.
00649     const ST one = Teuchos::as<ST> (1);
00650     // FIXME (mfh 20 Nov 2013) Should scale this 1.0e-6 term
00651     // appropriately for MT's machine precision.
00652     if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
00653       lambdaMinForApply_ = one;
00654       lambdaMaxForApply_ = lambdaMinForApply_;
00655       eigRatioForApply_ = one; // Ifpack doesn't include this line.
00656     }
00657   }
00658 }
00659 
00660 
00661 template<class ScalarType, class MV>
00662 ScalarType
00663 Chebyshev<ScalarType, MV>::
00664 getLambdaMaxForApply() const {
00665   return lambdaMaxForApply_;
00666 }
00667 
00668 
00669 template<class ScalarType, class MV>
00670 typename Chebyshev<ScalarType, MV>::MT
00671 Chebyshev<ScalarType, MV>::apply (const MV& B, MV& X)
00672 {
00673   TEUCHOS_TEST_FOR_EXCEPTION(
00674     A_.is_null (), std::runtime_error, "Ifpack2::Chebyshev::apply: The input "
00675     "matrix A is null.  Please call setMatrix() with a nonnull input matrix, "
00676     "and then call compute(), before calling this method.");
00677   TEUCHOS_TEST_FOR_EXCEPTION(
00678     STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
00679     "Ifpack2::Chebyshev::apply: There is no estimate for the max eigenvalue."
00680     << std::endl << computeBeforeApplyReminder);
00681   TEUCHOS_TEST_FOR_EXCEPTION(
00682     STS::isnaninf (lambdaMinForApply_), std::runtime_error,
00683     "Ifpack2::Chebyshev::apply: There is no estimate for the min eigenvalue."
00684     << std::endl << computeBeforeApplyReminder);
00685   TEUCHOS_TEST_FOR_EXCEPTION(
00686     STS::isnaninf (eigRatioForApply_), std::runtime_error,
00687     "Ifpack2::Chebyshev::apply: There is no estimate for the ratio of the max "
00688     "eigenvalue to the min eigenvalue."
00689     << std::endl << computeBeforeApplyReminder);
00690   TEUCHOS_TEST_FOR_EXCEPTION(
00691     D_.is_null (), std::runtime_error,
00692     "Ifpack2::Chebyshev::apply: "
00693     "The vector of inverse diagonal entries of the matrix has not yet been "
00694     "computed." << std::endl << computeBeforeApplyReminder);
00695 
00696   if (textbookAlgorithm_) {
00697     textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
00698                        lambdaMinForApply_, eigRatioForApply_, *D_);
00699   } else {
00700     ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
00701                      lambdaMinForApply_, eigRatioForApply_, *D_);
00702   }
00703 
00704   if (computeMaxResNorm_ && B.getNumVectors () > 0) {
00705     MV R (B.getMap (), B.getNumVectors ());
00706     computeResidual (R, B, *A_, X);
00707     Teuchos::Array<MT> norms (B.getNumVectors ());
00708     R.norm2 (norms ());
00709     return *std::max_element (norms.begin (), norms.end ());
00710   } else {
00711     return Teuchos::ScalarTraits<MT>::zero ();
00712   }
00713 }
00714 
00715 template<class ScalarType, class MV>
00716 void
00717 Chebyshev<ScalarType, MV>::
00718 print (std::ostream& out) {
00719   this->describe (* (Teuchos::getFancyOStream (Teuchos::rcpFromRef (out))),
00720                   Teuchos::VERB_MEDIUM);
00721 }
00722 
00723 template<class ScalarType, class MV>
00724 void
00725 Chebyshev<ScalarType, MV>::
00726 computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
00727                  const Teuchos::ETransp mode)
00728 {
00729   Tpetra::deep_copy(R, B);
00730   A.apply (X, R, mode, -STS::one(), STS::one());
00731 }
00732 
00733 template<class ScalarType, class MV>
00734 void
00735 Chebyshev<ScalarType, MV>::
00736 solve (MV& Z, const V& D_inv, const MV& R) {
00737   Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
00738 }
00739 
00740 template<class ScalarType, class MV>
00741 void
00742 Chebyshev<ScalarType, MV>::
00743 solve (MV& Z, const ST alpha, const V& D_inv, const MV& R) {
00744   Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
00745 }
00746 
00747 template<class ScalarType, class MV>
00748 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
00749 Chebyshev<ScalarType, MV>::
00750 makeInverseDiagonal (const row_matrix_type& A, const bool useDiagOffsets) const
00751 {
00752   using Teuchos::RCP;
00753   using Teuchos::rcpFromRef;
00754   using Teuchos::rcp_dynamic_cast;
00755 
00756   RCP<V> D_rowMap (new V (A.getGraph ()->getRowMap ()));
00757   if (useDiagOffsets) {
00758     // The optimizations below only work if A_ is a Tpetra::CrsMatrix.
00759     // We'll make our best guess about its type here, since we have no
00760     // way to get back the original fifth template parameter.
00761     typedef Tpetra::CrsMatrix<typename MV::scalar_type,
00762       typename MV::local_ordinal_type,
00763       typename MV::global_ordinal_type,
00764       typename MV::node_type> crs_matrix_type;
00765     RCP<const crs_matrix_type> A_crsMat =
00766       rcp_dynamic_cast<const crs_matrix_type> (rcpFromRef (A));
00767     if (! A_crsMat.is_null ()) {
00768       TEUCHOS_TEST_FOR_EXCEPTION(
00769         ! savedDiagOffsets_, std::logic_error,
00770         "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
00771         "It is not allowed to call this method with useDiagOffsets=true, "
00772         "if you have not previously saved offsets of diagonal entries.  "
00773         "This situation should never arise if this class is used properly.  "
00774         "Please report this bug to the Ifpack2 developers.");
00775       A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_ ());
00776     }
00777   }
00778   else {
00779     // This always works for a Tpetra::RowMatrix, even if it is not a
00780     // Tpetra::CrsMatrix.  We just don't have offsets in this case.
00781     A.getLocalDiagCopy (*D_rowMap);
00782   }
00783   RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
00784 
00785   // Invert the diagonal entries, replacing entries less (in
00786   // magnitude) than the user-specified value with that value.
00787   typedef KokkosClassic::MultiVector<ST, typename MV::node_type> KMV;
00788   KMV localDiag = D_rangeMap->getLocalMV ();
00789   typedef KokkosClassic::DefaultArithmetic<KMV> KMVT;
00790   KMVT::ReciprocalThreshold (localDiag, minDiagVal_);
00791   return Teuchos::rcp_const_cast<const V> (D_rangeMap);
00792 }
00793 
00794 
00795 template<class ScalarType, class MV>
00796 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
00797 Chebyshev<ScalarType, MV>::
00798 makeRangeMapVectorConst (const Teuchos::RCP<const V>& D) const
00799 {
00800   using Teuchos::RCP;
00801   using Teuchos::rcp;
00802   typedef Tpetra::Export<typename MV::local_ordinal_type,
00803                          typename MV::global_ordinal_type,
00804                          typename MV::node_type> export_type;
00805   // This throws logic_error instead of runtime_error, because the
00806   // methods that call makeRangeMapVector should all have checked
00807   // whether A_ is null before calling this method.
00808   TEUCHOS_TEST_FOR_EXCEPTION(
00809     A_.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
00810     "makeRangeMapVector: The input matrix A is null.  Please call setMatrix() "
00811     "with a nonnull input matrix before calling this method.  This is probably "
00812     "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
00813   TEUCHOS_TEST_FOR_EXCEPTION(
00814     D.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
00815     "makeRangeMapVector: The input Vector D is null.  This is probably "
00816     "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
00817 
00818   RCP<const map_type> sourceMap = D->getMap ();
00819   RCP<const map_type> rangeMap = A_->getRangeMap ();
00820   RCP<const map_type> rowMap = A_->getRowMap ();
00821 
00822   if (rangeMap->isSameAs (*sourceMap)) {
00823     // The given vector's Map is the same as the matrix's range Map.
00824     // That means we don't need to Export.  This is the fast path.
00825     return D;
00826   }
00827   else { // We need to Export.
00828     RCP<const export_type> exporter;
00829     // Making an Export object from scratch is expensive enough that
00830     // it's worth the O(1) global reductions to call isSameAs(), to
00831     // see if we can avoid that cost.
00832     if (sourceMap->isSameAs (*rowMap)) {
00833       // We can reuse the matrix's Export object, if there is one.
00834       exporter = A_->getGraph ()->getExporter ();
00835     }
00836     else { // We have to make a new Export object.
00837       exporter = rcp (new export_type (sourceMap, rangeMap));
00838     }
00839 
00840     if (exporter.is_null ()) {
00841       return D; // Row Map and range Map are the same; no need to Export.
00842     }
00843     else { // Row Map and range Map are _not_ the same; must Export.
00844       RCP<V> D_out = rcp (new V (*D));
00845       D_out->doExport (*D, *exporter, Tpetra::ADD);
00846       return Teuchos::rcp_const_cast<const V> (D_out);
00847     }
00848   }
00849 }
00850 
00851 
00852 template<class ScalarType, class MV>
00853 Teuchos::RCP<typename Chebyshev<ScalarType, MV>::V>
00854 Chebyshev<ScalarType, MV>::
00855 makeRangeMapVector (const Teuchos::RCP<V>& D) const
00856 {
00857   using Teuchos::rcp_const_cast;
00858   return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
00859 }
00860 
00861 
00862 template<class ScalarType, class MV>
00863 void
00864 Chebyshev<ScalarType, MV>::
00865 textbookApplyImpl (const op_type& A,
00866                    const MV& B,
00867                    MV& X,
00868                    const int numIters,
00869                    const ST lambdaMax,
00870                    const ST lambdaMin,
00871                    const ST eigRatio,
00872                    const V& D_inv) const
00873 {
00874   (void) lambdaMin; // Forestall compiler warning.
00875   const ST myLambdaMin = lambdaMax / eigRatio;
00876 
00877   const ST zero = Teuchos::as<ST> (0);
00878   const ST one = Teuchos::as<ST> (1);
00879   const ST two = Teuchos::as<ST> (2);
00880   const ST d = (lambdaMax + myLambdaMin) / two; // Ifpack2 calls this theta
00881   const ST c = (lambdaMax - myLambdaMin) / two; // Ifpack2 calls this 1/delta
00882 
00883   if (zeroStartingSolution_ && numIters > 0) {
00884     // If zero iterations, then input X is output X.
00885     X.putScalar (zero);
00886   }
00887   MV R (B.getMap (), B.getNumVectors (), false);
00888   MV P (B.getMap (), B.getNumVectors (), false);
00889   MV Z (B.getMap (), B.getNumVectors (), false);
00890   ST alpha, beta;
00891   for (int i = 0; i < numIters; ++i) {
00892     computeResidual (R, B, A, X); // R = B - A*X
00893     solve (Z, D_inv, R); // z = D_inv * R, that is, D \ R.
00894     if (i == 0) {
00895       P = Z;
00896       alpha = two / d;
00897     } else {
00898       //beta = (c * alpha / two)^2;
00899       //const ST sqrtBeta = c * alpha / two;
00900       //beta = sqrtBeta * sqrtBeta;
00901       beta = alpha * (c/two) * (c/two);
00902       alpha = one / (d - beta);
00903       P.update (one, Z, beta); // P = Z + beta*P
00904     }
00905     X.update (alpha, P, one); // X = X + alpha*P
00906     // If we compute the residual here, we could either do R = B -
00907     // A*X, or R = R - alpha*A*P.  Since we choose the former, we
00908     // can move the computeResidual call to the top of the loop.
00909   }
00910 }
00911 
00912 template<class ScalarType, class MV>
00913 typename Chebyshev<ScalarType, MV>::MT
00914 Chebyshev<ScalarType, MV>::maxNormInf (const MV& X) {
00915   Teuchos::Array<MT> norms (X.getNumVectors ());
00916   X.normInf (norms());
00917   return *std::max_element (norms.begin (), norms.end ());
00918 }
00919 
00920 template<class ScalarType, class MV>
00921 void
00922 Chebyshev<ScalarType, MV>::
00923 ifpackApplyImpl (const op_type& A,
00924                  const MV& B,
00925                  MV& X,
00926                  const int numIters,
00927                  const ST lambdaMax,
00928                  const ST lambdaMin,
00929                  const ST eigRatio,
00930                  const V& D_inv)
00931 {
00932 #ifdef HAVE_TEUCHOS_DEBUG
00933 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00934   using std::cerr;
00935   using std::endl;
00936   cerr << "\\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
00937   cerr << "\\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
00938 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00939 #endif // HAVE_TEUCHOS_DEBUG
00940 
00941   if (numIters <= 0) {
00942     return;
00943   }
00944   const ST one = Teuchos::as<ST> (1);
00945   const ST two = Teuchos::as<ST> (2);
00946 
00947   // Quick solve when the matrix A is the identity.
00948   if (lambdaMin == one && lambdaMax == lambdaMin) {
00949     solve (X, D_inv, B);
00950     return;
00951   }
00952 
00953   // Initialize coefficients
00954   const ST alpha = lambdaMax / eigRatio;
00955   const ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
00956   const ST delta = two / (beta - alpha);
00957   const ST theta = (beta + alpha) / two;
00958   const ST s1 = theta * delta;
00959 
00960 #ifdef HAVE_TEUCHOS_DEBUG
00961 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00962   cerr << "alpha = " << alpha << endl
00963        << "beta = " << beta << endl
00964        << "delta = " << delta << endl
00965        << "theta = " << theta << endl
00966        << "s1 = " << s1 << endl;
00967 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00968 #endif // HAVE_TEUCHOS_DEBUG
00969 
00970   // Fetch cached temporary vectors.
00971   Teuchos::RCP<MV> V_ptr, W_ptr;
00972   makeTempMultiVectors (V_ptr, W_ptr, B);
00973 
00974   // mfh 28 Jan 2013: We write V1 instead of V, so as not to confuse
00975   // the multivector V with the typedef V (for Tpetra::Vector).
00976   //MV V1 (B.getMap (), B.getNumVectors (), false);
00977   //MV W (B.getMap (), B.getNumVectors (), false);
00978   MV& V1 = *V_ptr;
00979   MV& W = *W_ptr;
00980 
00981 #ifdef HAVE_TEUCHOS_DEBUG
00982 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00983   cerr << "Iteration " << 1 << ":" << endl
00984        << "- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
00985 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00986 #endif // HAVE_TEUCHOS_DEBUG
00987 
00988   // Special case for the first iteration.
00989   if (! zeroStartingSolution_) {
00990     computeResidual (V1, B, A, X); // V1 = B - A*X
00991 
00992 #ifdef HAVE_TEUCHOS_DEBUG
00993 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
00994     cerr << "- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
00995 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
00996 #endif // HAVE_TEUCHOS_DEBUG
00997 
00998     solve (W, one/theta, D_inv, V1); // W = (1/theta)*D_inv*(B-A*X)
00999 
01000 #ifdef HAVE_TEUCHOS_DEBUG
01001 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
01002     cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
01003 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
01004 #endif // HAVE_TEUCHOS_DEBUG
01005 
01006     X.update (one, W, one); // X = X + W
01007   } else {
01008     solve (W, one/theta, D_inv, B); // W = (1/theta)*D_inv*B
01009 
01010 #ifdef HAVE_TEUCHOS_DEBUG
01011 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
01012     cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
01013 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
01014 #endif // HAVE_TEUCHOS_DEBUG
01015 
01016     Tpetra::deep_copy(X, W); // X = 0 + W
01017   }
01018 #ifdef HAVE_TEUCHOS_DEBUG
01019 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
01020   cerr << "- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
01021 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
01022 #endif // HAVE_TEUCHOS_DEBUG
01023 
01024   // The rest of the iterations.
01025   ST rhok = one / s1;
01026   ST rhokp1, dtemp1, dtemp2;
01027   for (int deg = 1; deg < numIters; ++deg) {
01028 
01029 #ifdef HAVE_TEUCHOS_DEBUG
01030 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
01031     cerr << "Iteration " << deg+1 << ":" << endl;
01032     cerr << "- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
01033     cerr << "- \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
01034     cerr << "- \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm () << endl;
01035     cerr << "- rhok = " << rhok << endl;
01036     V1.putScalar (STS::zero ());
01037 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
01038 #endif // HAVE_TEUCHOS_DEBUG
01039 
01040     computeResidual (V1, B, A, X); // V1 = B - A*X
01041 
01042 #ifdef HAVE_TEUCHOS_DEBUG
01043 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
01044     cerr << "- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
01045 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
01046 #endif // HAVE_TEUCHOS_DEBUG
01047 
01048     rhokp1 = one / (two * s1 - rhok);
01049     dtemp1 = rhokp1 * rhok;
01050     dtemp2 = two * rhokp1 * delta;
01051     rhok = rhokp1;
01052 
01053 #ifdef HAVE_TEUCHOS_DEBUG
01054 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
01055     cerr << "- dtemp1 = " << dtemp1 << endl
01056          << "- dtemp2 = " << dtemp2 << endl;
01057 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
01058 #endif // HAVE_TEUCHOS_DEBUG
01059 
01060     W.scale (dtemp1);
01061     W.elementWiseMultiply (dtemp2, D_inv, V1, one);
01062     X.update (one, W, one);
01063 
01064 #ifdef HAVE_TEUCHOS_DEBUG
01065 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
01066     cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
01067     cerr << "- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
01068 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
01069 #endif // HAVE_TEUCHOS_DEBUG
01070   }
01071 }
01072 
01073 template<class ScalarType, class MV>
01074 typename Chebyshev<ScalarType, MV>::ST
01075 Chebyshev<ScalarType, MV>::
01076 powerMethod (const op_type& A, const V& D_inv, const int numIters)
01077 {
01078   const ST zero = Teuchos::as<ST> (0);
01079   const ST one = Teuchos::as<ST> (1);
01080   ST lambdaMax = zero;
01081   ST RQ_top, RQ_bottom, norm;
01082 
01083   V x (A.getDomainMap ());
01084   V y (A.getRangeMap ());
01085   x.randomize ();
01086   norm = x.norm2 ();
01087   TEUCHOS_TEST_FOR_EXCEPTION(norm == zero, std::runtime_error,
01088     "Ifpack2::Chebyshev::powerMethod: "
01089     "Tpetra::Vector's randomize() method filled the vector "
01090     "with zeros.  This is not impossible, but is unlikely.  "
01091     "It's far more likely that there is a bug in Tpetra.");
01092 
01093   x.scale (one / norm);
01094   for (int iter = 0; iter < numIters; ++iter) {
01095     A.apply (x, y);
01096     solve (y, D_inv, y);
01097     RQ_top = y.dot (x);
01098     RQ_bottom = x.dot (x);
01099     lambdaMax = RQ_top / RQ_bottom;
01100     norm = y.norm2 ();
01101     if (norm == zero) { // Return something reasonable.
01102       return zero;
01103     }
01104     x.update (one / norm, y, zero);
01105   }
01106   return lambdaMax;
01107 }
01108 
01109 template<class ScalarType, class MV>
01110 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
01111 Chebyshev<ScalarType, MV>::getMatrix () const {
01112   return A_;
01113 }
01114 
01115 template<class ScalarType, class MV>
01116 bool
01117 Chebyshev<ScalarType, MV>::
01118 hasTransposeApply () const {
01119   // Technically, this is true, because the matrix must be symmetric.
01120   return true;
01121 }
01122 
01123 template<class ScalarType, class MV>
01124 void
01125 Chebyshev<ScalarType, MV>::
01126 makeTempMultiVectors (Teuchos::RCP<MV>& V1,
01127                       Teuchos::RCP<MV>& W,
01128                       const MV& X)
01129 {
01130   if (V_.is_null ()) {
01131     V_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), false));
01132   }
01133   //W must be initialized to zero when it is used as a multigrid smoother.
01134   if (W_.is_null ()) {
01135     W_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), true));
01136   }
01137   V1 = V_;
01138   W = W_;
01139 }
01140 
01141 template<class ScalarType, class MV>
01142 std::string
01143 Chebyshev<ScalarType, MV>::
01144 description () const {
01145   std::ostringstream oss;
01146   // YAML requires quoting the key in this case, to distinguish
01147   // key's colons from the colon that separates key from value.
01148   oss << "\"Ifpack2::Details::Chebyshev\":"
01149       << "{"
01150       << "degree: " << numIters_
01151       << ", lambdaMax: " << lambdaMaxForApply_
01152       << ", alpha: " << eigRatioForApply_
01153       << ", lambdaMin: " << lambdaMinForApply_
01154       << "}";
01155   return oss.str();
01156 }
01157 
01158 template<class ScalarType, class MV>
01159 void
01160 Chebyshev<ScalarType, MV>::
01161 describe (Teuchos::FancyOStream& out,
01162           const Teuchos::EVerbosityLevel verbLevel) const
01163 {
01164   using Teuchos::TypeNameTraits;
01165   using std::endl;
01166 
01167   const Teuchos::EVerbosityLevel vl =
01168     (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
01169   if (vl > Teuchos::VERB_NONE) {
01170     if (vl == Teuchos::VERB_LOW) {
01171       out << description () << endl;
01172     } else { // vl > Teuchos::VERB_LOW
01173       // YAML requires quoting the key in this case, to distinguish
01174       // key's colons from the colon that separates key from value.
01175       out << "\"Ifpack2::Details::Chebyshev\":" << endl;
01176       Teuchos::OSTab tab1 (Teuchos::rcpFromRef (out));
01177       out << "Template parameters:" << endl;
01178       {
01179         Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
01180         out << "ScalarType: \"" << TypeNameTraits<ScalarType>::name () << "\"" << endl
01181             << "MV: \"" << TypeNameTraits<MV>::name () << "\"" << endl;
01182       }
01183       // "Computed parameters" literally means "parameters whose
01184       // values were computed by compute()."
01185       out << endl << "Computed parameters:" << endl;
01186       {
01187         Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
01188         // Users might want to see the values in the computed inverse
01189         // diagonal, so we print them out at the highest verbosity.
01190         out << "D_: ";
01191         if (D_.is_null ()) {
01192           out << "unset" << endl;
01193         } else if (vl <= Teuchos::VERB_HIGH) {
01194           out << "set" << endl;
01195         } else { // D_ not null and vl > Teuchos::VERB_HIGH
01196           out << endl;
01197           {
01198             Teuchos::OSTab tab3 (Teuchos::rcpFromRef (out));
01199             D_->describe (out, vl);
01200           }
01201           out << endl;
01202         }
01203         // V_ and W_ are scratch space; their values are irrelevant.
01204         // All that matters is whether or not they have been set.
01205         out << "V_: " << (V_.is_null () ? "unset" : "set") << endl
01206             << "W_: " << (W_.is_null () ? "unset" : "set") << endl
01207             << "computedLambdaMax_: " << computedLambdaMax_ << endl
01208             << "computedLambdaMin_: " << computedLambdaMin_ << endl
01209             << "lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
01210             << "lambdaMinForApply_: " << lambdaMinForApply_ << endl
01211             << "eigRatioForApply_: " << eigRatioForApply_ << endl;
01212       }
01213       out << "User parameters:" << endl;
01214       {
01215         Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
01216         out << "userInvDiag_: ";
01217         if (userInvDiag_.is_null ()) {
01218           out << "unset" << endl;
01219         } else if (vl <= Teuchos::VERB_HIGH) {
01220           out << "set" << endl;
01221         } else { // userInvDiag_ not null and vl > Teuchos::VERB_HIGH
01222           out << endl;
01223           {
01224             Teuchos::OSTab tab3 (Teuchos::rcpFromRef (out));
01225             userInvDiag_->describe (out, vl);
01226           }
01227           out << endl;
01228         }
01229         out << "userLambdaMax_: " << userLambdaMax_ << endl
01230             << "userLambdaMin_: " << userLambdaMin_ << endl
01231             << "userEigRatio_: " << userEigRatio_ << endl
01232             << "numIters_: " << numIters_ << endl
01233             << "eigMaxIters_: " << eigMaxIters_ << endl
01234             << "zeroStartingSolution_: " << zeroStartingSolution_ << endl
01235             << "assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
01236             << "textbookAlgorithm_: " << textbookAlgorithm_ << endl
01237             << "computeMaxResNorm_: " << computeMaxResNorm_ << endl;
01238       }
01239       out << endl;
01240     }
01241   }
01242 }
01243 
01244 } // namespace Details
01245 } // namespace Ifpack2
01246 
01247 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
01248   template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
01249 
01250 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends