Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Chebyshev_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Tempated 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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 */
00042 
00043 #ifndef IFPACK2_CHEBYSHEV_DEF_HPP
00044 #define IFPACK2_CHEBYSHEV_DEF_HPP
00045 
00046 #include <Ifpack2_Condest.hpp>
00047 #include <Ifpack2_Parameters.hpp>
00048 #include <Teuchos_TimeMonitor.hpp>
00049 
00050 namespace Ifpack2 {
00051 
00052 template<class MatrixType>
00053 Chebyshev<MatrixType>::
00054 Chebyshev (const Teuchos::RCP<const row_matrix_type>& A)
00055   : impl_ (A),
00056     Condest_ ( -Teuchos::ScalarTraits<magnitude_type>::one() ),
00057     IsInitialized_ (false),
00058     IsComputed_ (false),
00059     NumInitialize_ (0),
00060     NumCompute_ (0),
00061     NumApply_ (0),
00062     InitializeTime_ (0.0),
00063     ComputeTime_ (0.0),
00064     ApplyTime_ (0.0),
00065     ComputeFlops_ (0.0),
00066     ApplyFlops_ (0.0)
00067 {
00068   this->setObjectLabel ("Ifpack2::Chebyshev");
00069 }
00070 
00071 
00072 template<class MatrixType>
00073 Chebyshev<MatrixType>::~Chebyshev() {
00074 }
00075 
00076 
00077 template<class MatrixType>
00078 void Chebyshev<MatrixType>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
00079 {
00080   if (A.getRawPtr () != impl_.getMatrix ().getRawPtr ()) {
00081     IsInitialized_ = false;
00082     IsComputed_ = false;
00083     impl_.setMatrix (A);
00084   }
00085 }
00086 
00087 
00088 template<class MatrixType>
00089 void
00090 Chebyshev<MatrixType>::setParameters (const Teuchos::ParameterList& List)
00091 {
00092   // FIXME (mfh 25 Jan 2013) Casting away const is bad here.
00093   impl_.setParameters (const_cast<Teuchos::ParameterList&> (List));
00094 }
00095 
00096 
00097 template<class MatrixType>
00098 Teuchos::RCP<const Teuchos::Comm<int> >
00099 Chebyshev<MatrixType>::getComm () const
00100 {
00101   Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
00102   TEUCHOS_TEST_FOR_EXCEPTION(
00103     A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getComm: The input "
00104     "matrix A is null.  Please call setMatrix() with a nonnull input matrix "
00105     "before calling this method.");
00106   return A->getRowMap ()->getComm ();
00107 }
00108 
00109 
00110 template<class MatrixType>
00111 Teuchos::RCP<const typename Chebyshev<MatrixType>::row_matrix_type>
00112 Chebyshev<MatrixType>::
00113 getMatrix() const {
00114   return impl_.getMatrix ();
00115 }
00116 
00117 
00118 template<class MatrixType>
00119 Teuchos::RCP<const MatrixType>
00120 Chebyshev<MatrixType>::
00121 getCrsMatrix() const {
00122   return Teuchos::rcp_dynamic_cast<const MatrixType> (impl_.getMatrix ());
00123 }
00124 
00125 
00126 template<class MatrixType>
00127 Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
00128 Chebyshev<MatrixType>::
00129 getDomainMap () const
00130 {
00131   Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
00132   TEUCHOS_TEST_FOR_EXCEPTION(
00133     A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getDomainMap: The "
00134     "input matrix A is null.  Please call setMatrix() with a nonnull input "
00135     "matrix before calling this method.");
00136   return A->getDomainMap ();
00137 }
00138 
00139 
00140 template<class MatrixType>
00141 Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
00142 Chebyshev<MatrixType>::
00143 getRangeMap () const
00144 {
00145   Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
00146   TEUCHOS_TEST_FOR_EXCEPTION(
00147     A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getRangeMap: The "
00148     "input matrix A is null.  Please call setMatrix() with a nonnull input "
00149     "matrix before calling this method.");
00150   return A->getRangeMap ();
00151 }
00152 
00153 
00154 template<class MatrixType>
00155 bool Chebyshev<MatrixType>::hasTransposeApply() const {
00156   return impl_.hasTransposeApply ();
00157 }
00158 
00159 
00160 template<class MatrixType>
00161 int Chebyshev<MatrixType>::getNumInitialize() const {
00162   return NumInitialize_;
00163 }
00164 
00165 
00166 template<class MatrixType>
00167 int Chebyshev<MatrixType>::getNumCompute() const {
00168   return NumCompute_;
00169 }
00170 
00171 
00172 template<class MatrixType>
00173 int Chebyshev<MatrixType>::getNumApply() const {
00174   return NumApply_;
00175 }
00176 
00177 
00178 template<class MatrixType>
00179 double Chebyshev<MatrixType>::getInitializeTime() const {
00180   return InitializeTime_;
00181 }
00182 
00183 
00184 template<class MatrixType>
00185 double Chebyshev<MatrixType>::getComputeTime() const {
00186   return ComputeTime_;
00187 }
00188 
00189 
00190 template<class MatrixType>
00191 double Chebyshev<MatrixType>::getApplyTime() const {
00192   return ApplyTime_;
00193 }
00194 
00195 
00196 template<class MatrixType>
00197 double Chebyshev<MatrixType>::getComputeFlops () const {
00198   return ComputeFlops_;
00199 }
00200 
00201 
00202 template<class MatrixType>
00203 double Chebyshev<MatrixType>::getApplyFlops () const {
00204   return ApplyFlops_;
00205 }
00206 
00207 
00208 template<class MatrixType>
00209 typename Chebyshev<MatrixType>::magnitude_type
00210 Chebyshev<MatrixType>::getCondEst () const {
00211   return Condest_;
00212 }
00213 
00214 
00215 template<class MatrixType>
00216 typename Chebyshev<MatrixType>::magnitude_type
00217 Chebyshev<MatrixType>::
00218 computeCondEst (CondestType CT,
00219                 local_ordinal_type MaxIters,
00220                 magnitude_type Tol,
00221                 const Teuchos::Ptr<const row_matrix_type>& matrix)
00222 {
00223   if (! isComputed ()) {
00224     return -Teuchos::ScalarTraits<magnitude_type>::one ();
00225   }
00226   else {
00227     // Always compute it. Call Condest() with no parameters to get
00228     // the previous estimate.
00229     Condest_ = Ifpack2::Condest (*this, CT, MaxIters, Tol, matrix);
00230     return Condest_;
00231   }
00232 }
00233 
00234 
00235 template<class MatrixType>
00236 void
00237 Chebyshev<MatrixType>::
00238 apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
00239        Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
00240        Teuchos::ETransp mode,
00241        scalar_type alpha,
00242        scalar_type beta) const
00243 {
00244   const std::string timerName ("Ifpack2::Chebyshev::apply");
00245   Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
00246   if (timer.is_null ()) {
00247     timer = Teuchos::TimeMonitor::getNewCounter (timerName);
00248   }
00249 
00250   // Start timing here.
00251   {
00252     Teuchos::TimeMonitor timeMon (*timer);
00253 
00254     // compute() calls initialize() if it hasn't already been called.
00255     // Thus, we only need to check isComputed().
00256     TEUCHOS_TEST_FOR_EXCEPTION(
00257       ! isComputed (), std::runtime_error,
00258       "Ifpack2::Chebyshev::apply(): You must call the compute() method before "
00259       "you may call apply().");
00260     TEUCHOS_TEST_FOR_EXCEPTION(
00261       X.getNumVectors () != Y.getNumVectors (), std::runtime_error,
00262       "Ifpack2::Chebyshev::apply(): X and Y must have the same number of "
00263       "columns.  X.getNumVectors() = " << X.getNumVectors() << " != "
00264       << "Y.getNumVectors() = " << Y.getNumVectors() << ".");
00265     applyImpl (X, Y, mode, alpha, beta);
00266   }
00267   ++NumApply_;
00268 
00269   // timer->totalElapsedTime() returns the total time over all timer
00270   // calls.  Thus, we use = instead of +=.
00271   ApplyTime_ = timer->totalElapsedTime ();
00272 }
00273 
00274 
00275 template<class MatrixType>
00276 void
00277 Chebyshev<MatrixType>::
00278 applyMat (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
00279           Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
00280           Teuchos::ETransp mode) const
00281 {
00282   TEUCHOS_TEST_FOR_EXCEPTION(
00283     X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
00284     "Ifpack2::Chebyshev::applyMat: X.getNumVectors() != Y.getNumVectors().");
00285 
00286   Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
00287   TEUCHOS_TEST_FOR_EXCEPTION(
00288     A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::applyMat: The input "
00289     "matrix A is null.  Please call setMatrix() with a nonnull input matrix "
00290     "before calling this method.");
00291 
00292   A->apply (X, Y, mode);
00293 }
00294 
00295 
00296 template<class MatrixType>
00297 void Chebyshev<MatrixType>::initialize () {
00298   // We create the timer, but this method doesn't do anything, so
00299   // there is no need to start the timer.  The resulting total time
00300   // will always be zero.
00301   const std::string timerName ("Ifpack2::Chebyshev::initialize");
00302   Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
00303   if (timer.is_null ()) {
00304     timer = Teuchos::TimeMonitor::getNewCounter (timerName);
00305   }
00306   IsInitialized_ = true;
00307   ++NumInitialize_;
00308 }
00309 
00310 
00311 template<class MatrixType>
00312 void Chebyshev<MatrixType>::compute ()
00313 {
00314   const std::string timerName ("Ifpack2::Chebyshev::compute");
00315   Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
00316   if (timer.is_null ()) {
00317     timer = Teuchos::TimeMonitor::getNewCounter (timerName);
00318   }
00319 
00320   // Start timing here.
00321   {
00322     Teuchos::TimeMonitor timeMon (*timer);
00323     if (! isInitialized ()) {
00324       initialize ();
00325     }
00326     IsComputed_ = false;
00327     Condest_ = -Teuchos::ScalarTraits<magnitude_type>::one();
00328     impl_.compute ();
00329   }
00330   IsComputed_ = true;
00331   ++NumCompute_;
00332 
00333   // timer->totalElapsedTime() returns the total time over all timer
00334   // calls.  Thus, we use = instead of +=.
00335   ComputeTime_ = timer->totalElapsedTime ();
00336 }
00337 
00338 
00339 template<class MatrixType>
00340 void
00341 Chebyshev<MatrixType>::
00342 PowerMethod (const Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Operator,
00343              const Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& InvPointDiagonal,
00344              const int MaximumIterations,
00345              scalar_type& lambda_max)
00346 {
00347   TEUCHOS_TEST_FOR_EXCEPTION(
00348     true, std::logic_error, "Ifpack2::Chebyshev::PowerMethod: "
00349     "This method is deprecated.  Please do not call it.");
00350 
00351   // const scalar_type one = STS::one();
00352   // const scalar_type zero = STS::zero();
00353 
00354   // lambda_max = zero;
00355   // Teuchos::Array<scalar_type> RQ_top(1), RQ_bottom(1);
00356   // vector_type x (Operator.getDomainMap ());
00357   // vector_type y (Operator.getRangeMap ());
00358   // x.randomize ();
00359   // Teuchos::Array<magnitude_type> norms (x.getNumVectors ());
00360   // x.norm2 (norms ());
00361   // x.scale (one / norms[0]);
00362 
00363   // for (int iter = 0; iter < MaximumIterations; ++iter) {
00364   //   Operator.apply (x, y);
00365   //   y.elementWiseMultiply (one, InvPointDiagonal, y, zero);
00366   //   y.dot (x, RQ_top ());
00367   //   x.dot (x, RQ_bottom ());
00368   //   lambda_max = RQ_top[0] / RQ_bottom[0];
00369   //   y.norm2 (norms ());
00370   //   TEUCHOS_TEST_FOR_EXCEPTION(
00371   //     norms[0] == zero,
00372   //     std::runtime_error,
00373   //     "Ifpack2::Chebyshev::PowerMethod: norm == 0 at iteration " << (iter+1)
00374   //     << " of " << MaximumIterations);
00375   //   x.update (one / norms[0], y, zero);
00376   // }
00377 }
00378 
00379 
00380 template<class MatrixType>
00381 void Chebyshev<MatrixType>::
00382 CG (const Tpetra::Operator<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Operator,
00383     const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& InvPointDiagonal,
00384    const int MaximumIterations,
00385    scalar_type& lambda_min, scalar_type& lambda_max)
00386 {
00387   TEUCHOS_TEST_FOR_EXCEPTION(
00388     true, std::logic_error,
00389     "Ifpack2::Chebyshev::CG: Not implemented.  "
00390     "Please use Belos' implementation of CG with Tpetra objects.");
00391 }
00392 
00393 
00394 template <class MatrixType>
00395 std::string Chebyshev<MatrixType>::description () const {
00396   std::ostringstream out;
00397 
00398   // Output is a valid YAML dictionary in flow style.  If you don't
00399   // like everything on a single line, you should call describe()
00400   // instead.
00401   out << "\"Ifpack2::Chebyshev\": {";
00402   out << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
00403       << "Computed: " << (isComputed () ? "true" : "false") << ", ";
00404 
00405   out << impl_.description() << ", ";
00406 
00407   if (impl_.getMatrix ().is_null ()) {
00408     out << "Matrix: null";
00409   }
00410   else {
00411     out << "Global matrix dimensions: ["
00412         << impl_.getMatrix ()->getGlobalNumRows () << ", "
00413         << impl_.getMatrix ()->getGlobalNumCols () << "]"
00414         << ", Global nnz: " << impl_.getMatrix ()->getGlobalNumEntries();
00415   }
00416 
00417   out << "}";
00418   return out.str ();
00419 }
00420 
00421 
00422 template <class MatrixType>
00423 void Chebyshev<MatrixType>::
00424 describe (Teuchos::FancyOStream &out,
00425           const Teuchos::EVerbosityLevel verbLevel) const
00426 {
00427   const Teuchos::EVerbosityLevel vl =
00428     (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
00429   const int myRank = this->getComm ()->getRank ();
00430 
00431   if (vl != Teuchos::VERB_NONE && myRank == 0) {
00432     // By convention, describe() starts with a tab.
00433     Teuchos::OSTab tab0 (out);
00434     out << description ();
00435   }
00436 
00437 #if 0
00438   using Teuchos::Comm;
00439   using Teuchos::RCP;
00440   using Teuchos::VERB_DEFAULT;
00441   using Teuchos::VERB_NONE;
00442   using Teuchos::VERB_LOW;
00443   using Teuchos::VERB_MEDIUM;
00444   using Teuchos::VERB_HIGH;
00445   using Teuchos::VERB_EXTREME;
00446   using std::endl;
00447   using std::setw;
00448 
00449   Teuchos::EVerbosityLevel vl = verbLevel;
00450   if (vl == VERB_DEFAULT) {
00451     vl = VERB_LOW;
00452   }
00453   RCP<const Comm<int> > comm = A_->getRowMap ()->getComm ();
00454 
00455   const int myImageID = comm->getRank();
00456   Teuchos::OSTab tab(out);
00457 
00458   scalar_type MinVal, MaxVal;
00459   if (IsComputed_) {
00460     Teuchos::ArrayRCP<const scalar_type> DiagView = InvDiagonal_->get1dView();
00461     scalar_type myMinVal = DiagView[0];
00462     scalar_type myMaxVal = DiagView[0];
00463     for(typename Teuchos::ArrayRCP<scalar_type>::size_type i=1; i<DiagView.size(); ++i) {
00464       if (STS::magnitude(myMinVal) > STS::magnitude(DiagView[i])) myMinVal = DiagView[i];
00465       if (STS::magnitude(myMaxVal) < STS::magnitude(DiagView[i])) myMaxVal = DiagView[i];
00466     }
00467     Teuchos::reduceAll(*comm, Teuchos::REDUCE_MIN, 1, &myMinVal, &MinVal);
00468     Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, 1, &myMaxVal, &MaxVal);
00469   }
00470 
00471   //    none: print nothing
00472   //     low: print O(1) info from node 0
00473   //  medium:
00474   //    high:
00475   // extreme:
00476   if (vl != VERB_NONE && myImageID == 0) {
00477     out << this->description() << endl;
00478     out << endl;
00479     out << "===============================================================================" << std::endl;
00480     out << "Degree of polynomial      = " << PolyDegree_ << std::endl;
00481     if   (ZeroStartingSolution_) { out << "Using zero starting solution" << endl; }
00482     else                         { out << "Using input starting solution" << endl; }
00483     if   (Condest_ == - Teuchos::ScalarTraits<magnitude_type>::one()) {
00484       out << "Condition number estimate       = N/A" << endl;
00485     }
00486     else                    { out << "Condition number estimate       = " << Condest_ << endl; }
00487     if (IsComputed_) {
00488       out << "Minimum value on stored inverse diagonal = " << MinVal << std::endl;
00489       out << "Maximum value on stored inverse diagonal = " << MaxVal << std::endl;
00490     }
00491     out << std::endl;
00492     out << "Phase           # calls    Total Time (s)     Total MFlops      MFlops/s       " << endl;
00493     out << "------------    -------    ---------------    ---------------   ---------------" << endl;
00494     out << setw(12) << "initialize()" << setw(5) << getNumInitialize() << "    " << setw(15) << getInitializeTime() << endl;
00495     out << setw(12) << "compute()" << setw(5) << getNumCompute()    << "    " << setw(15) << getComputeTime() << "    "
00496         << setw(15) << getComputeFlops() << "    "
00497         << setw(15) << (getComputeTime() != 0.0 ? getComputeFlops() / getComputeTime() * 1.0e-6 : 0.0) << endl;
00498     out << setw(12) << "apply()" << setw(5) << getNumApply()    << "    " << setw(15) << getApplyTime() << "    "
00499         << setw(15) << getApplyFlops() << "    "
00500         << setw(15) << (getApplyTime() != 0.0 ? getApplyFlops() / getApplyTime() * 1.0e-6 : 0.0) << endl;
00501     out << "===============================================================================" << std::endl;
00502     out << endl;
00503   }
00504 #endif // 0
00505 }
00506 
00507 template<class MatrixType>
00508 void
00509 Chebyshev<MatrixType>::
00510 applyImpl (const MV& X,
00511            MV& Y,
00512            Teuchos::ETransp mode,
00513            scalar_type alpha,
00514            scalar_type beta) const
00515 {
00516   using Teuchos::ArrayRCP;
00517   using Teuchos::as;
00518   using Teuchos::RCP;
00519   using Teuchos::rcp;
00520   using Teuchos::rcp_const_cast;
00521   using Teuchos::rcpFromRef;
00522 
00523   const scalar_type zero = STS::zero();
00524   const scalar_type one = STS::one();
00525 
00526   // Y = beta*Y + alpha*M*X.
00527 
00528   // If alpha == 0, then we don't need to do Chebyshev at all.
00529   if (alpha == zero) {
00530     if (beta == zero) { // Obey Sparse BLAS rules; avoid 0*NaN.
00531       Y.putScalar (zero);
00532     }
00533     else {
00534       Y.scale (beta);
00535     }
00536     return;
00537   }
00538 
00539   // If beta != 0, then we need to keep a copy of the initial value of
00540   // Y, so that we can add beta*it to the Chebyshev result at the end.
00541   // Usually this method is called with beta == 0, so we don't have to
00542   // worry about caching Y_org.
00543   RCP<MV> Y_orig;
00544   if (beta != zero) {
00545     Y_orig = rcp (new MV (createCopy(Y)));
00546   }
00547 
00548   // If X and Y point to the same memory location, we need to use a
00549   // copy of X (X_copy) as the input MV.  Otherwise, just let X_copy
00550   // point to X.
00551   //
00552   // This is hopefully an uncommon use case, so we don't bother to
00553   // optimize for it by caching X_copy.
00554   RCP<const MV> X_copy;
00555   bool copiedInput = false;
00556   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues()) {
00557     X_copy = rcp (new MV (createCopy(X)));
00558     copiedInput = true;
00559   }
00560   else {
00561     X_copy = rcpFromRef (X);
00562   }
00563 
00564   // If alpha != 1, fold alpha into (a copy of) X.
00565   //
00566   // This is an uncommon use case, so we don't bother to optimize for
00567   // it by caching X_copy.  However, we do check whether we've already
00568   // copied X above, to avoid a second copy.
00569   if (alpha != one) {
00570     RCP<MV> X_copy_nonConst = rcp_const_cast<MV> (X_copy);
00571     if (! copiedInput) {
00572       X_copy_nonConst = rcp (new MV (createCopy(X)));
00573       copiedInput = true;
00574     }
00575     X_copy_nonConst->scale (alpha);
00576     X_copy = rcp_const_cast<const MV> (X_copy_nonConst);
00577   }
00578 
00579   impl_.apply (*X_copy, Y);
00580 
00581   if (beta != zero) {
00582     Y.update (beta, *Y_orig, one); // Y = beta * Y_orig + 1 * Y
00583   }
00584 }
00585 
00586 
00587 template<class MatrixType>
00588 typename MatrixType::scalar_type Chebyshev<MatrixType>::getLambdaMaxForApply () const {
00589   return impl_.getLambdaMaxForApply ();
00590 }
00591 
00592 
00593 
00594 }//namespace Ifpack2
00595 
00596 #define IFPACK2_CHEBYSHEV_INSTANT(S,LO,GO,N)                            \
00597   template class Ifpack2::Chebyshev< Tpetra::CrsMatrix<S, LO, GO, N> >; \
00598   template class Ifpack2::Chebyshev< Tpetra::RowMatrix<S, LO, GO, N> >;
00599 
00600 #endif // IFPACK2_CHEBYSHEV_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends