Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Chebyshev_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #ifndef IFPACK2_CHEBYSHEV_DEF_HPP
00031 #define IFPACK2_CHEBYSHEV_DEF_HPP
00032 
00033 
00034 namespace Ifpack2 {
00035 
00036 template<class MatrixType>
00037 Chebyshev<MatrixType>::
00038 Chebyshev (const Teuchos::RCP<const row_matrix_type>& A)
00039   : impl_ (A),
00040     Time_ (Teuchos::rcp (new Teuchos::Time ("Ifpack2::Chebyshev"))),
00041     Condest_ (-1.0),
00042     IsInitialized_ (false),
00043     IsComputed_ (false),
00044     NumInitialize_ (0),
00045     NumCompute_ (0),
00046     NumApply_ (0),
00047     InitializeTime_ (0.0),
00048     ComputeTime_ (0.0),
00049     ApplyTime_ (0.0),
00050     ComputeFlops_ (0.0),
00051     ApplyFlops_ (0.0)
00052 {}
00053 
00054 //==========================================================================
00055 template<class MatrixType>
00056 Chebyshev<MatrixType>::~Chebyshev() {
00057 }
00058 
00059 //==========================================================================
00060 template<class MatrixType>
00061 void 
00062 Chebyshev<MatrixType>::setParameters (const Teuchos::ParameterList& List) 
00063 {
00064   // FIXME (mfh 25 Jan 2013) Casting away const is bad here.
00065   impl_.setParameters (const_cast<Teuchos::ParameterList&> (List));
00066 }
00067 
00068 
00069 //==========================================================================
00070 template<class MatrixType>
00071 const Teuchos::RCP<const Teuchos::Comm<int> > & 
00072 Chebyshev<MatrixType>::getComm() const {
00073   return impl_.getMatrix ()->getRowMap ()->getComm ();
00074 }
00075 
00076 //==========================================================================
00077 template<class MatrixType>
00078 Teuchos::RCP<const typename Chebyshev<MatrixType>::row_matrix_type>
00079 Chebyshev<MatrixType>::
00080 getMatrix() const {
00081   return impl_.getMatrix ();
00082 }
00083 
00084 //==========================================================================
00085 template<class MatrixType>
00086 const Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>&
00087 Chebyshev<MatrixType>::
00088 getDomainMap() const {
00089   return impl_.getMatrix ()->getDomainMap ();
00090 }
00091 
00092 //==========================================================================
00093 template<class MatrixType>
00094 const Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>&
00095 Chebyshev<MatrixType>::
00096 getRangeMap() const {
00097   return impl_.getMatrix ()->getRangeMap ();
00098 }
00099 
00100 //==========================================================================
00101 template<class MatrixType>
00102 bool Chebyshev<MatrixType>::hasTransposeApply() const {
00103   return impl_.hasTransposeApply ();
00104 }
00105 
00106 //==========================================================================
00107 template<class MatrixType>
00108 int Chebyshev<MatrixType>::getNumInitialize() const {
00109   return(NumInitialize_);
00110 }
00111 
00112 //==========================================================================
00113 template<class MatrixType>
00114 int Chebyshev<MatrixType>::getNumCompute() const {
00115   return(NumCompute_);
00116 }
00117 
00118 //==========================================================================
00119 template<class MatrixType>
00120 int Chebyshev<MatrixType>::getNumApply() const {
00121   return(NumApply_);
00122 }
00123 
00124 //==========================================================================
00125 template<class MatrixType>
00126 double Chebyshev<MatrixType>::getInitializeTime() const {
00127   return(InitializeTime_);
00128 }
00129 
00130 //==========================================================================
00131 template<class MatrixType>
00132 double Chebyshev<MatrixType>::getComputeTime() const {
00133   return(ComputeTime_);
00134 }
00135 
00136 //==========================================================================
00137 template<class MatrixType>
00138 double Chebyshev<MatrixType>::getApplyTime() const {
00139   return(ApplyTime_);
00140 }
00141 
00142 //==========================================================================
00143 template<class MatrixType>
00144 double Chebyshev<MatrixType>::getComputeFlops() const {
00145   return(ComputeFlops_);
00146 }
00147 
00148 //==========================================================================
00149 template<class MatrixType>
00150 double Chebyshev<MatrixType>::getApplyFlops() const {
00151   return(ApplyFlops_);
00152 }
00153 
00154 //==========================================================================
00155 template<class MatrixType>
00156 typename Chebyshev<MatrixType>::magnitude_type
00157 Chebyshev<MatrixType>::
00158 getCondEst() const {
00159   return(Condest_);
00160 }
00161 
00162 //==========================================================================
00163 template<class MatrixType>
00164 typename Chebyshev<MatrixType>::magnitude_type
00165 Chebyshev<MatrixType>::
00166 computeCondEst (CondestType CT,
00167     local_ordinal_type MaxIters, 
00168     magnitude_type Tol,
00169     const Teuchos::Ptr<const row_matrix_type>& matrix) 
00170 {
00171   if (! isComputed ()) {
00172     return -Teuchos::ScalarTraits<magnitude_type>::one ();
00173   }
00174   else {
00175     // Always compute it. Call Condest() with no parameters to get
00176     // the previous estimate.
00177     Condest_ = Ifpack2::Condest (*this, CT, MaxIters, Tol, matrix);
00178     return Condest_;
00179   }
00180 }
00181 
00182 
00183 template<class MatrixType>
00184 void 
00185 Chebyshev<MatrixType>::
00186 apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
00187        Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
00188        Teuchos::ETransp mode,
00189        scalar_type alpha,
00190        scalar_type beta) const 
00191 {
00192   {
00193     Teuchos::TimeMonitor timeMon (*Time_);
00194 
00195     // compute() calls initialize() if it hasn't already been called.
00196     // Thus, we only need to check isComputed().
00197     TEUCHOS_TEST_FOR_EXCEPTION(! isComputed(), std::runtime_error, 
00198       "Ifpack2::Chebyshev::apply(): You must call the compute() method before "
00199       "you may call apply().");
00200     TEUCHOS_TEST_FOR_EXCEPTION(
00201        X.getNumVectors() != Y.getNumVectors(), 
00202        std::runtime_error,
00203        "Ifpack2::Chebyshev::apply(): X and Y must have the same number of "
00204        "columns.  X.getNumVectors() = " << X.getNumVectors() << " != "
00205        << "Y.getNumVectors() = " << Y.getNumVectors() << ".");
00206 #ifdef HAVE_TEUCHOS_DEBUG
00207     {
00208       // The relation 'isSameAs' is transitive.  It's also a collective,
00209       // so we don't have to do a "shared" test for exception (i.e., a
00210       // global reduction on the test value).
00211       TEUCHOS_TEST_FOR_EXCEPTION(
00212          ! X.getMap ()->isSameAs (*getDomainMap ()),
00213          std::runtime_error,
00214          "Ifpack2::Chebyshev: The domain Map of the matrix must be the same as "
00215    "the Map of the input vector(s) X.");
00216       TEUCHOS_TEST_FOR_EXCEPTION(
00217          ! Y.getMap ()->isSameAs (*getRangeMap ()),
00218          std::runtime_error,
00219          "Ifpack2::Chebyshev: The range Map of the matrix must be the same as "
00220    "the Map of the output vector(s) Y.");
00221     }
00222 #endif // HAVE_TEUCHOS_DEBUG
00223     applyImpl (X, Y, mode, alpha, beta);
00224   }
00225   ++NumApply_;
00226   ApplyTime_ += Time_->totalElapsedTime ();
00227 }
00228 
00229 
00230 template<class MatrixType>
00231 void 
00232 Chebyshev<MatrixType>::
00233 applyMat (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
00234     Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
00235     Teuchos::ETransp mode) const
00236 {
00237   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00238    "Ifpack2::Chebyshev::applyMat(): X.getNumVectors() != Y.getNumVectors().");
00239   impl_.getMatrix ()->apply (X, Y, mode);
00240 }
00241 
00242 
00243 template<class MatrixType>
00244 void Chebyshev<MatrixType>::initialize() {
00245   // This method doesn't do anything anymore, so there's no need to time it.
00246   IsInitialized_ = true;
00247   ++NumInitialize_;
00248   // Note to developers: Defer fetching any data that relate to the
00249   // structure of the matrix until compute().  That way, it will
00250   // always be correct to omit calling initialize(), even if the
00251   // number of entries in the sparse matrix has changed since the last
00252   // call to compute().
00253 }
00254 
00255 
00256 template<class MatrixType>
00257 void Chebyshev<MatrixType>::compute()
00258 {
00259   {
00260     Teuchos::TimeMonitor timeMon (*Time_);
00261     if (! isInitialized ()) {
00262       initialize ();
00263     }
00264     IsComputed_ = false;
00265     Condest_ = -1.0;  
00266     impl_.compute ();
00267   }
00268   IsComputed_ = true;
00269   ++NumCompute_;
00270   ComputeTime_ += Time_->totalElapsedTime();
00271 }
00272 
00273 
00274 template<class MatrixType>
00275 void 
00276 Chebyshev<MatrixType>::
00277 PowerMethod (const Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Operator, 
00278        const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& InvPointDiagonal, 
00279        const int MaximumIterations, 
00280        scalar_type& lambda_max)
00281 {
00282   const scalar_type one = STS::one();
00283   const scalar_type zero = STS::zero();
00284 
00285   lambda_max = zero;
00286   Teuchos::Array<scalar_type> RQ_top(1), RQ_bottom(1);
00287   vector_type x (Operator.getDomainMap ());
00288   vector_type y (Operator.getRangeMap ());
00289   x.randomize ();
00290   Teuchos::Array<magnitude_type> norms (x.getNumVectors ());
00291   x.norm2 (norms ());
00292   x.scale (one / norms[0]);
00293 
00294   for (int iter = 0; iter < MaximumIterations; ++iter) {
00295     Operator.apply (x, y);
00296     y.elementWiseMultiply (one, InvPointDiagonal, y, zero);
00297     y.dot (x, RQ_top ());
00298     x.dot (x, RQ_bottom ());
00299     lambda_max = RQ_top[0] / RQ_bottom[0];
00300     y.norm2 (norms ());
00301     TEUCHOS_TEST_FOR_EXCEPTION(
00302       norms[0] == zero, 
00303       std::runtime_error, 
00304       "Ifpack2::Chebyshev::PowerMethod: norm == 0 at iteration " << (iter+1) 
00305       << " of " << MaximumIterations);
00306     x.update (one / norms[0], y, zero);
00307   }
00308 }
00309 
00310 //==========================================================================
00311 template<class MatrixType>
00312 void Chebyshev<MatrixType>::
00313 CG(const Tpetra::Operator<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Operator, 
00314             const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& InvPointDiagonal, 
00315    const int MaximumIterations, 
00316    scalar_type& lambda_min, scalar_type& lambda_max)
00317 {
00318   TEUCHOS_TEST_FOR_EXCEPTION(
00319     true, std::logic_error,
00320     "Ifpack2::Chebyshev::CG: Not implemented.  "
00321     "Please use Belos' implementation of CG with Tpetra objects.");
00322 }
00323 
00324 //==========================================================================
00325 template <class MatrixType>
00326 std::string Chebyshev<MatrixType>::description() const {
00327   std::ostringstream oss;
00328   oss << Teuchos::Describable::description();
00329   if (isInitialized()) {
00330     if (isComputed()) {
00331       oss << "{status = initialized, computed";
00332     }
00333     else {
00334       oss << "{status = initialized, not computed";
00335     }
00336   }
00337   else {
00338     oss << "{status = not initialized, not computed";
00339   }
00340   //
00341   oss << ", global rows = " << impl_.getMatrix ()->getGlobalNumRows()
00342       << ", global cols = " << impl_.getMatrix ()->getGlobalNumCols()
00343       << "}";
00344   return oss.str();
00345 }
00346 
00347 //==========================================================================
00348 template <class MatrixType>
00349 void Chebyshev<MatrixType>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00350 #if 0
00351   using Teuchos::Comm;
00352   using Teuchos::RCP;
00353   using Teuchos::VERB_DEFAULT;
00354   using Teuchos::VERB_NONE;
00355   using Teuchos::VERB_LOW;
00356   using Teuchos::VERB_MEDIUM;
00357   using Teuchos::VERB_HIGH;
00358   using Teuchos::VERB_EXTREME;
00359   using std::endl;
00360   using std::setw;
00361 
00362   Teuchos::EVerbosityLevel vl = verbLevel;
00363   if (vl == VERB_DEFAULT) {
00364     vl = VERB_LOW;
00365   }
00366   RCP<const Comm<int> > comm = A_->getRowMap ()->getComm ();
00367 
00368   const int myImageID = comm->getRank();
00369   Teuchos::OSTab tab(out);
00370 
00371   scalar_type MinVal, MaxVal;
00372   if (IsComputed_) {
00373     Teuchos::ArrayRCP<const scalar_type> DiagView = InvDiagonal_->get1dView();
00374     scalar_type myMinVal = DiagView[0];
00375     scalar_type myMaxVal = DiagView[0];
00376     for(typename Teuchos::ArrayRCP<scalar_type>::size_type i=1; i<DiagView.size(); ++i) {
00377       if (STS::magnitude(myMinVal) > STS::magnitude(DiagView[i])) myMinVal = DiagView[i];
00378       if (STS::magnitude(myMaxVal) < STS::magnitude(DiagView[i])) myMaxVal = DiagView[i];
00379     }
00380     Teuchos::reduceAll(*comm, Teuchos::REDUCE_MIN, 1, &myMinVal, &MinVal);
00381     Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, 1, &myMaxVal, &MaxVal);
00382   }
00383 
00384   //    none: print nothing
00385   //     low: print O(1) info from node 0
00386   //  medium: 
00387   //    high: 
00388   // extreme: 
00389   if (vl != VERB_NONE && myImageID == 0) {
00390     out << this->description() << endl;
00391     out << endl;
00392     out << "===============================================================================" << std::endl;
00393     out << "Degree of polynomial      = " << PolyDegree_ << std::endl;
00394     if   (ZeroStartingSolution_) { out << "Using zero starting solution" << endl; }
00395     else                         { out << "Using input starting solution" << endl; }
00396     if   (Condest_ == -1.0) { out << "Condition number estimate       = N/A" << endl; }
00397     else                    { out << "Condition number estimate       = " << Condest_ << endl; }
00398     if (IsComputed_) {
00399       out << "Minimum value on stored inverse diagonal = " << MinVal << std::endl;
00400       out << "Maximum value on stored inverse diagonal = " << MaxVal << std::endl;
00401     }
00402     out << std::endl;
00403     out << "Phase           # calls    Total Time (s)     Total MFlops      MFlops/s       " << endl;
00404     out << "------------    -------    ---------------    ---------------   ---------------" << endl;
00405     out << setw(12) << "initialize()" << setw(5) << getNumInitialize() << "    " << setw(15) << getInitializeTime() << endl;
00406     out << setw(12) << "compute()" << setw(5) << getNumCompute()    << "    " << setw(15) << getComputeTime() << "    " 
00407         << setw(15) << getComputeFlops() << "    " 
00408         << setw(15) << (getComputeTime() != 0.0 ? getComputeFlops() / getComputeTime() * 1.0e-6 : 0.0) << endl;
00409     out << setw(12) << "apply()" << setw(5) << getNumApply()    << "    " << setw(15) << getApplyTime() << "    " 
00410         << setw(15) << getApplyFlops() << "    " 
00411         << setw(15) << (getApplyTime() != 0.0 ? getApplyFlops() / getApplyTime() * 1.0e-6 : 0.0) << endl;
00412     out << "===============================================================================" << std::endl;
00413     out << endl;
00414   }
00415 #endif // 0
00416 }
00417 
00418 template<class MatrixType>
00419 void 
00420 Chebyshev<MatrixType>::
00421 applyImpl (const MV& X,
00422      MV& Y,
00423      Teuchos::ETransp mode,
00424      scalar_type alpha,
00425      scalar_type beta) const 
00426 {
00427   using Teuchos::ArrayRCP;
00428   using Teuchos::as;
00429   using Teuchos::RCP;
00430   using Teuchos::rcp;
00431   using Teuchos::rcp_const_cast;
00432   using Teuchos::rcpFromRef;
00433 
00434   const scalar_type zero = STS::zero();
00435   const scalar_type one = STS::one();
00436 
00437   // Y = beta*Y + alpha*M*X.
00438 
00439   // If alpha == 0, then we don't need to do Chebyshev at all.
00440   if (alpha == zero) {
00441     if (beta == zero) { // Obey Sparse BLAS rules; avoid 0*NaN.
00442       Y.putScalar (zero);
00443     }
00444     else {
00445       Y.scale (beta);
00446     }
00447     return;
00448   }
00449 
00450   // If beta != 0, then we need to keep a copy of the initial value of
00451   // Y, so that we can add beta*it to the Chebyshev result at the end.
00452   // Usually this method is called with beta == 0, so we don't have to 
00453   // worry about caching Y_org.
00454   RCP<MV> Y_orig;
00455   if (beta != zero) {
00456     Y_orig = rcp (new MV (Y));
00457   }
00458 
00459   // If X and Y point to the same memory location, we need to use a
00460   // copy of X (X_copy) as the input MV.  Otherwise, just let X_copy
00461   // point to X.
00462   //
00463   // This is hopefully an uncommon use case, so we don't bother to
00464   // optimize for it by caching X_copy.
00465   RCP<const MV> X_copy;
00466   bool copiedInput = false;
00467   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues()) {
00468     X_copy = rcp (new MV (X));
00469     copiedInput = true;
00470   }
00471   else {
00472     X_copy = rcpFromRef (X);
00473   }
00474   
00475   // If alpha != 1, fold alpha into (a copy of) X.
00476   //
00477   // This is an uncommon use case, so we don't bother to optimize for
00478   // it by caching X_copy.  However, we do check whether we've already
00479   // copied X above, to avoid a second copy.
00480   if (alpha != one) {
00481     RCP<MV> X_copy_nonConst = rcp_const_cast<MV> (X_copy);
00482     if (! copiedInput) {
00483       X_copy_nonConst = rcp (new MV (X));
00484       copiedInput = true;
00485     }
00486     X_copy_nonConst->scale (alpha);
00487     X_copy = rcp_const_cast<const MV> (X_copy_nonConst);
00488   }
00489 
00490   impl_.apply (*X_copy, Y);
00491 
00492   if (beta != zero) {
00493     Y.update (beta, *Y_orig, one); // Y = beta * Y_orig + 1 * Y
00494   }
00495 }
00496 
00497 
00498 }//namespace Ifpack2
00499 
00500 #endif // IFPACK2_CHEBYSHEV_DEF_HPP
00501 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends