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 // 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 //Definitions for the Chebyshev methods:
00037 
00038 //==========================================================================
00039 template<class MatrixType>
00040 Chebyshev<MatrixType>::Chebyshev(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A)
00041 : A_(A),
00042   Comm_(A->getRowMap()->getComm()),
00043   Time_( Teuchos::rcp( new Teuchos::Time("Ifpack2::Chebyshev") ) ),
00044   PolyDegree_(1),
00045   EigRatio_(30.0),
00046   LambdaMin_(0.0),
00047   LambdaMax_(100.0),
00048   MinDiagonalValue_(0.0),
00049   ZeroStartingSolution_(true),
00050   Condest_(-1.0),
00051   IsInitialized_(false),
00052   IsComputed_(false),
00053   NumInitialize_(0),
00054   NumCompute_(0),
00055   NumApply_(0),
00056   InitializeTime_(0.0),
00057   ComputeTime_(0.0),
00058   ApplyTime_(0.0),
00059   ComputeFlops_(0.0),
00060   ApplyFlops_(0.0),
00061   NumMyRows_(0),
00062   NumGlobalRows_(0),
00063   NumGlobalNonzeros_(0)
00064 {
00065   TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error, 
00066       Teuchos::typeName(*this) << "::Chebyshev(): input matrix reference was null.");
00067 }
00068 
00069 //==========================================================================
00070 template<class MatrixType>
00071 Chebyshev<MatrixType>::~Chebyshev() {
00072 }
00073 
00074 //==========================================================================
00075 template<class MatrixType>
00076 void Chebyshev<MatrixType>::setParameters(const Teuchos::ParameterList& List) {
00077 
00078   Ifpack2::getParameter(List, "chebyshev: ratio eigenvalue", EigRatio_);
00079   Ifpack2::getParameter(List, "chebyshev: min eigenvalue", LambdaMin_);
00080   Ifpack2::getParameter(List, "chebyshev: max eigenvalue", LambdaMax_);
00081   Ifpack2::getParameter(List, "chebyshev: degree",PolyDegree_);
00082   Ifpack2::getParameter(List, "chebyshev: min diagonal value", MinDiagonalValue_);
00083   Ifpack2::getParameter(List, "chebyshev: zero starting solution", ZeroStartingSolution_);
00084 
00085   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>* ID = 0;
00086   Ifpack2::getParameter(List, "chebyshev: operator inv diagonal", ID);
00087 
00088   if (ID != 0) {
00089     InvDiagonal_ = Teuchos::rcp( new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(*ID) );
00090   }
00091 }
00092 
00093 //==========================================================================
00094 template<class MatrixType>
00095 const Teuchos::RCP<const Teuchos::Comm<int> > & 
00096 Chebyshev<MatrixType>::getComm() const{
00097   return(Comm_);
00098 }
00099 
00100 //==========================================================================
00101 template<class MatrixType>
00102 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
00103 Chebyshev<MatrixType>::getMatrix() const {
00104   return(A_);
00105 }
00106 
00107 //==========================================================================
00108 template<class MatrixType>
00109 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >&
00110 Chebyshev<MatrixType>::getDomainMap() const {
00111   return A_->getDomainMap();
00112 }
00113 
00114 //==========================================================================
00115 template<class MatrixType>
00116 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >&
00117 Chebyshev<MatrixType>::getRangeMap() const {
00118   return A_->getRangeMap();
00119 }
00120 
00121 //==========================================================================
00122 template<class MatrixType>
00123 bool Chebyshev<MatrixType>::hasTransposeApply() const {
00124   return true;
00125 }
00126 
00127 //==========================================================================
00128 template<class MatrixType>
00129 int Chebyshev<MatrixType>::getNumInitialize() const {
00130   return(NumInitialize_);
00131 }
00132 
00133 //==========================================================================
00134 template<class MatrixType>
00135 int Chebyshev<MatrixType>::getNumCompute() const {
00136   return(NumCompute_);
00137 }
00138 
00139 //==========================================================================
00140 template<class MatrixType>
00141 int Chebyshev<MatrixType>::getNumApply() const {
00142   return(NumApply_);
00143 }
00144 
00145 //==========================================================================
00146 template<class MatrixType>
00147 double Chebyshev<MatrixType>::getInitializeTime() const {
00148   return(InitializeTime_);
00149 }
00150 
00151 //==========================================================================
00152 template<class MatrixType>
00153 double Chebyshev<MatrixType>::getComputeTime() const {
00154   return(ComputeTime_);
00155 }
00156 
00157 //==========================================================================
00158 template<class MatrixType>
00159 double Chebyshev<MatrixType>::getApplyTime() const {
00160   return(ApplyTime_);
00161 }
00162 
00163 //==========================================================================
00164 template<class MatrixType>
00165 double Chebyshev<MatrixType>::getComputeFlops() const {
00166   return(ComputeFlops_);
00167 }
00168 
00169 //==========================================================================
00170 template<class MatrixType>
00171 double Chebyshev<MatrixType>::getApplyFlops() const {
00172   return(ApplyFlops_);
00173 }
00174 
00175 //==========================================================================
00176 template<class MatrixType>
00177 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00178 Chebyshev<MatrixType>::getCondEst() const {
00179   return(Condest_);
00180 }
00181 
00182 //==========================================================================
00183 template<class MatrixType>
00184 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00185 Chebyshev<MatrixType>::computeCondEst(
00186                      CondestType CT,
00187                      LocalOrdinal MaxIters, 
00188                      magnitudeType Tol,
00189                      const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &matrix) {
00190   if (!isComputed()) // cannot compute right now
00191     return(-1.0);
00192 
00193   // always compute it. Call Condest() with no parameters to get
00194   // the previous estimate.
00195   Condest_ = Ifpack2::Condest(*this, CT, MaxIters, Tol, matrix);
00196 
00197   return(Condest_);
00198 }
00199 
00200 //==========================================================================
00201 template<class MatrixType>
00202 void Chebyshev<MatrixType>::apply(
00203            const Tpetra::MultiVector<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& X,
00204                  Tpetra::MultiVector<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& Y,
00205                 Teuchos::ETransp mode,
00206                  Scalar alpha,
00207                  Scalar beta) const {
00208   TEUCHOS_TEST_FOR_EXCEPTION(!isComputed(), std::runtime_error, 
00209       "Ifpack2::Chebyshev::apply() ERROR, not yet computed.");
00210 
00211   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00212      "Ifpack2::Chebyshev::apply() ERROR: X.getNumVectors() != Y.getNumVectors().");
00213 
00214   Time_->start();
00215 
00216   // If X and Y are pointing to the same memory location,
00217   // we need to create an auxiliary vector, Xcopy
00218   Teuchos::RCP< const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xcopy;
00219   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues())
00220     Xcopy = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X) );
00221   else
00222     Xcopy = Teuchos::rcp( &X, false );
00223 
00224   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > xView = Xcopy->get2dView();
00225   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > yView = Y.get2dViewNonConst();
00226   Teuchos::ArrayRCP<const Scalar> invdiag = InvDiagonal_->get1dView();
00227 
00228   size_t nVecs = Y.getNumVectors();
00229 
00230   //--- Do a quick solve when the matrix is identity
00231   if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
00232     if (nVecs == 1) {
00233       Teuchos::ArrayRCP<Scalar> y = yView[0];
00234       Teuchos::ArrayRCP<const Scalar> x = xView[0];
00235       for (size_t i = 0; i < NumMyRows_; ++i)
00236         y[i] = x[i]*invdiag[i];
00237     }
00238     else {
00239       for (size_t i = 0; i < NumMyRows_; ++i) {
00240         const Scalar& coeff = invdiag[i];
00241         for (size_t k = 0; k < nVecs; ++k)
00242           yView[k][i] = xView[k][i] * coeff;
00243       }
00244     } // if (nVec == 1)
00245     return;
00246   } // if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_))
00247 
00248   //--- initialize coefficients
00249   // Note that delta stores the inverse of ML_Cheby::delta
00250   Scalar alpha_cheby = LambdaMax_ / EigRatio_;
00251   Scalar beta_cheby = 1.1 * LambdaMax_;
00252   Scalar delta = 2.0 / (beta_cheby - alpha_cheby);
00253   Scalar theta = 0.5 * (beta_cheby + alpha_cheby);
00254   Scalar s1 = theta * delta;
00255 
00256   //--- Define vectors
00257   // In ML_Cheby, V corresponds to pAux and W to dk
00258   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> V(X);
00259   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> W(X);
00260 
00261   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > vView = V.get2dView();
00262   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > wView = W.get2dViewNonConst();
00263 
00264   Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00265 
00266   Scalar oneOverTheta = one/theta;
00267 
00268   // Do the smoothing when block scaling is turned OFF
00269   // --- Treat the initial guess
00270   if (ZeroStartingSolution_ == false) {
00271     A_->apply(Y, V);
00272     // compute W = invDiag * ( X - V )/ Theta
00273     if (nVecs == 1) {
00274       Teuchos::ArrayRCP<const Scalar> x = xView[0];
00275       Teuchos::ArrayRCP<Scalar> w = wView[0];
00276       Teuchos::ArrayRCP<const Scalar> v = vView[0];
00277       for (size_t i = 0; i < NumMyRows_; ++i)
00278         w[i] = invdiag[i] * (x[i] - v[i]) * oneOverTheta;
00279     }
00280     else {
00281       for (size_t k = 0; k < nVecs; ++k) {
00282         Teuchos::ArrayRCP<Scalar> wk = wView[k];
00283         Teuchos::ArrayRCP<const Scalar> vk = vView[k];
00284         for (size_t i = 0; i < NumMyRows_; ++i) {
00285           Scalar coeff = invdiag[i]*oneOverTheta;
00286           wk[i] = (xView[k][i] - (vk[i])) * coeff;
00287         }
00288       }
00289     } // if (nVec == 1)
00290     // Update the vector Y
00291     Y.update(one, W, one);
00292   }
00293   else {
00294     // compute W = invDiag * X / Theta
00295     if (nVecs == 1) {
00296       Teuchos::ArrayRCP<const Scalar> x= xView[0];
00297       Teuchos::ArrayRCP<Scalar> w = wView[0];
00298       Teuchos::ArrayRCP<Scalar> y = yView[0];
00299       for (size_t i = 0; i < NumMyRows_; ++i) {
00300         w[i] = invdiag[i] * x[i] * oneOverTheta;
00301         y[i] = w[i];
00302       }
00303     }
00304     else {
00305       for (size_t k = 0; k < nVecs; ++k) {
00306         for (size_t i = 0; i < NumMyRows_; ++i) {
00307           Scalar coeff = invdiag[i]*oneOverTheta;
00308           wView[k][i] = xView[k][i] * coeff;
00309           yView[k][i] = wView[k][i];
00310         }
00311       }
00312     } // if (nVec == 1)
00313   } // if (ZeroStartingSolution_ == false)
00314 
00315   //--- Apply the polynomial
00316   Scalar rhok = 1.0/s1, rhokp1;
00317   Scalar dtemp1, dtemp2;
00318   int degreeMinusOne = PolyDegree_ - 1;
00319   Scalar two = 2.0;
00320   for (int deg = 0; deg < degreeMinusOne; ++deg) {
00321     A_->apply(Y, V);
00322     rhokp1 = one / (two *s1 - rhok);
00323     dtemp1 = rhokp1 * rhok;
00324     dtemp2 = two * rhokp1 * delta;
00325     rhok = rhokp1;
00326     // compute W = dtemp1 * W
00327     W.scale(dtemp1);
00328     // compute W = W + dtemp2 * invDiag * ( X - V )
00329     for (size_t k = 0; k < nVecs; ++k) {
00330       for (size_t i = 0; i < NumMyRows_; ++i) {
00331         Scalar coeff = invdiag[i]*dtemp2;
00332         wView[k][i] += (xView[k][i] - (vView[k][i])) * coeff;
00333       }
00334     }
00335     // Update the vector Y
00336     Y.update(one, W, one);
00337   } // for (deg = 0; deg < degreeMinusOne; ++deg)
00338 
00339   // Flops are updated in each of the following. 
00340 
00341   ++NumApply_;
00342   Time_->stop();
00343   ApplyTime_ += Time_->totalElapsedTime();
00344 }
00345 
00346 //==========================================================================
00347 template<class MatrixType>
00348 void Chebyshev<MatrixType>::applyMat(
00349            const Tpetra::MultiVector<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& X,
00350                  Tpetra::MultiVector<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& Y,
00351              Teuchos::ETransp mode) const
00352 {
00353   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == false, std::runtime_error,
00354      "Ifpack2::Chebyshev::applyMat() ERROR: isComputed() must be true prior to calling applyMat().");
00355   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00356      "Ifpack2::Chebyshev::applyMat() ERROR: X.getNumVectors() != Y.getNumVectors().");
00357   A_->apply(X, Y, mode);
00358 }
00359 
00360 //==========================================================================
00361 template<class MatrixType>
00362 void Chebyshev<MatrixType>::initialize() {
00363   IsInitialized_ = false;
00364 
00365   TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error, 
00366       "Ifpack2::Chebyshev::initialize ERROR: A_ == Teuchos::null");
00367 
00368   TEUCHOS_TEST_FOR_EXCEPTION(A_->getGlobalNumRows() != A_->getGlobalNumCols(), std::runtime_error,
00369      "Ifpack2::Chebyshev::initialize ERROR: only square matrices are supported");
00370 
00371   NumMyRows_ = A_->getNodeNumRows();
00372   NumGlobalRows_ = A_->getGlobalNumRows();
00373   NumGlobalNonzeros_ = A_->getGlobalNumEntries();
00374 
00375   ++NumInitialize_;
00376   Time_->stop();
00377   InitializeTime_ += Time_->totalElapsedTime();
00378   IsInitialized_ = true;
00379 }
00380 
00381 //==========================================================================
00382 template<class MatrixType>
00383 void Chebyshev<MatrixType>::compute()
00384 {
00385   if (!isInitialized()) {
00386     initialize();
00387   }
00388 
00389   Time_->start(true);
00390 
00391   // reset values
00392   IsComputed_ = false;
00393   Condest_ = -1.0;
00394 
00395   TEUCHOS_TEST_FOR_EXCEPTION(PolyDegree_ <= 0, std::runtime_error,
00396     "Ifpack2::Chebyshev::compute() ERROR: PolyDegree_ must be at least one");
00397   
00398   if (InvDiagonal_ == Teuchos::null)
00399   {
00400     InvDiagonal_ = Teuchos::rcp( new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A_->getRowMap()) );
00401 
00402     A_->getLocalDiagCopy(*InvDiagonal_);
00403 
00404     // Inverse diagonal elements
00405     // Replace zeros with 1.0
00406     Teuchos::ArrayRCP<Scalar> diagvals = InvDiagonal_->get1dViewNonConst();
00407     for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00408       Scalar& diag = diagvals[i];
00409       if (Teuchos::ScalarTraits<Scalar>::magnitude(diag) < Teuchos::ScalarTraits<Scalar>::magnitude(MinDiagonalValue_))
00410         diag = MinDiagonalValue_;
00411       else
00412         diag = 1.0 / diag;
00413     }
00414   }
00415   // otherwise the inverse of the diagonal has been given by the user
00416 
00417   ComputeFlops_ += NumMyRows_;
00418 
00419   ++NumCompute_;
00420   Time_->stop();
00421   ComputeTime_ += Time_->totalElapsedTime();
00422   IsComputed_ = true;
00423 }
00424 
00425 //==========================================================================
00426 template<class MatrixType>
00427 void Chebyshev<MatrixType>::
00428 PowerMethod(const Tpetra::Operator<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& Operator, 
00429             const Tpetra::Vector<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& InvPointDiagonal, 
00430             const int MaximumIterations, 
00431             typename MatrixType::scalar_type& lambda_max)
00432 {
00433   // this is a simple power method
00434   lambda_max = 0.0;
00435   Teuchos::Array<Scalar> RQ_top(1), RQ_bottom(1);
00436   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> x(Operator.getDomainMap());
00437   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> y(Operator.getRangeMap());
00438   x.randomize();
00439   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
00440   Teuchos::Array<magnitudeType> norms(x.getNumVectors());
00441   x.norm2(norms());
00442 
00443   x.scale(1.0 / norms[0]);
00444 
00445   Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00446   Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00447 
00448   for (int iter = 0; iter < MaximumIterations; ++iter)
00449   {
00450     Operator.apply(x, y);
00451     y.elementWiseMultiply(1.0, InvPointDiagonal, y, 0.0);
00452     y.dot(x, RQ_top());
00453     x.dot(x, RQ_bottom());
00454     lambda_max = RQ_top[0] / RQ_bottom[0];
00455     y.norm2(norms());
00456     TEUCHOS_TEST_FOR_EXCEPTION(norms[0] == zero, std::runtime_error, "Ifpack2::Chebyshev::PowerMethod ERROR, norm == 0");
00457     x.update( one / norms[0], y, zero);
00458   }
00459 }
00460 
00461 //==========================================================================
00462 template<class MatrixType>
00463 void Chebyshev<MatrixType>::
00464 CG(const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Operator, 
00465             const Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& InvPointDiagonal, 
00466    const int MaximumIterations, 
00467    Scalar& lambda_min, Scalar& lambda_max)
00468 {
00469 #ifdef HAVE_IFPACK2_AZTECOO
00470   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> x(Operator.getDomainMap());
00471   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> y(Operator.getRangeMap());
00472   x.Random();
00473   y.putScalar(0.0);
00474 
00475   Tpetra::LinearProblem LP(const_cast<Tpetra::Operator*>(&Operator), &x, &y);
00476   AztecOO solver(LP);
00477   solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
00478   solver.SetAztecOption(AZ_output, AZ_none);
00479 
00480   Ifpack2_DiagPreconditioner diag(Operator.OperatorDomainMap(),
00481                                  Operator.OperatorRangeMap(),
00482                                  InvPointDiagonal);
00483   solver.SetPrecOperator(&diag);
00484   solver.Iterate(MaximumIterations, 1e-10);
00485 
00486   const double* status = solver.GetAztecStatus();
00487 
00488   lambda_min = status[AZ_lambda_min];
00489   lambda_max = status[AZ_lambda_max];
00490 
00491   return(0);
00492 #else
00493   throw std::runtime_error("Ifpack2::Chebyshev::CG: support for AztecOO not currently implemented.");
00494 #endif
00495 }
00496 
00497 //==========================================================================
00498 template <class MatrixType>
00499 std::string Chebyshev<MatrixType>::description() const {
00500   std::ostringstream oss;
00501   oss << Teuchos::Describable::description();
00502   if (isInitialized()) {
00503     if (isComputed()) {
00504       oss << "{status = initialized, computed";
00505     }
00506     else {
00507       oss << "{status = initialized, not computed";
00508     }
00509   }
00510   else {
00511     oss << "{status = not initialized, not computed";
00512   }
00513   //
00514   oss << ", global rows = " << A_->getGlobalNumRows()
00515       << ", global cols = " << A_->getGlobalNumCols()
00516       << "}";
00517   return oss.str();
00518 }
00519 
00520 //==========================================================================
00521 template <class MatrixType>
00522 void Chebyshev<MatrixType>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00523   using std::endl;
00524   using std::setw;
00525   using Teuchos::VERB_DEFAULT;
00526   using Teuchos::VERB_NONE;
00527   using Teuchos::VERB_LOW;
00528   using Teuchos::VERB_MEDIUM;
00529   using Teuchos::VERB_HIGH;
00530   using Teuchos::VERB_EXTREME;
00531   Teuchos::EVerbosityLevel vl = verbLevel;
00532   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00533   const int myImageID = Comm_->getRank();
00534   Teuchos::OSTab tab(out);
00535 
00536   Scalar MinVal, MaxVal;
00537   if (IsComputed_) {
00538     Teuchos::ArrayRCP<const Scalar> DiagView = InvDiagonal_->get1dView();
00539     Scalar myMinVal = DiagView[0];
00540     Scalar myMaxVal = DiagView[0];
00541     for(typename Teuchos::ArrayRCP<Scalar>::size_type i=1; i<DiagView.size(); ++i) {
00542       if (Teuchos::ScalarTraits<Scalar>::magnitude(myMinVal) > Teuchos::ScalarTraits<Scalar>::magnitude(DiagView[i])) myMinVal = DiagView[i];
00543       if (Teuchos::ScalarTraits<Scalar>::magnitude(myMaxVal) < Teuchos::ScalarTraits<Scalar>::magnitude(DiagView[i])) myMaxVal = DiagView[i];
00544     }
00545     Teuchos::reduceAll(*Comm_, Teuchos::REDUCE_MIN, 1, &myMinVal, &MinVal);
00546     Teuchos::reduceAll(*Comm_, Teuchos::REDUCE_MAX, 1, &myMaxVal, &MaxVal);
00547   }
00548 
00549   //    none: print nothing
00550   //     low: print O(1) info from node 0
00551   //  medium: 
00552   //    high: 
00553   // extreme: 
00554   if (vl != VERB_NONE && myImageID == 0) {
00555     out << this->description() << endl;
00556     out << endl;
00557     out << "===============================================================================" << std::endl;
00558     out << "Degree of polynomial      = " << PolyDegree_ << std::endl;
00559     if   (ZeroStartingSolution_) { out << "Using zero starting solution" << endl; }
00560     else                         { out << "Using input starting solution" << endl; }
00561     if   (Condest_ == -1.0) { out << "Condition number estimate       = N/A" << endl; }
00562     else                    { out << "Condition number estimate       = " << Condest_ << endl; }
00563     if (IsComputed_) {
00564       out << "Minimum value on stored inverse diagonal = " << MinVal << std::endl;
00565       out << "Maximum value on stored inverse diagonal = " << MaxVal << std::endl;
00566     }
00567     out << std::endl;
00568     out << "Phase           # calls    Total Time (s)     Total MFlops      MFlops/s       " << endl;
00569     out << "------------    -------    ---------------    ---------------   ---------------" << endl;
00570     out << setw(12) << "initialize()" << setw(5) << getNumInitialize() << "    " << setw(15) << getInitializeTime() << endl;
00571     out << setw(12) << "compute()" << setw(5) << getNumCompute()    << "    " << setw(15) << getComputeTime() << "    " 
00572         << setw(15) << getComputeFlops() << "    " 
00573         << setw(15) << (getComputeTime() != 0.0 ? getComputeFlops() / getComputeTime() * 1.0e-6 : 0.0) << endl;
00574     out << setw(12) << "apply()" << setw(5) << getNumApply()    << "    " << setw(15) << getApplyTime() << "    " 
00575         << setw(15) << getApplyFlops() << "    " 
00576         << setw(15) << (getApplyTime() != 0.0 ? getApplyFlops() / getApplyTime() * 1.0e-6 : 0.0) << endl;
00577     out << "===============================================================================" << std::endl;
00578     out << endl;
00579   }
00580 }
00581 
00582 }//namespace Ifpack2
00583 
00584 #endif // IFPACK2_CHEBYSHEV_DEF_HPP
00585 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends