IFPACK Development
Ifpack_Chebyshev.cpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) 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 #include "Ifpack_ConfigDefs.h"
00031 #include <iomanip>
00032 #include "Epetra_Operator.h"
00033 #include "Epetra_RowMatrix.h"
00034 #include "Epetra_Comm.h"
00035 #include "Epetra_Map.h"
00036 #include "Epetra_MultiVector.h"
00037 #include "Epetra_Vector.h"
00038 #include "Epetra_Time.h"
00039 #include "Ifpack_Chebyshev.h"
00040 #include "Ifpack_Utils.h"
00041 #include "Ifpack_Condest.h"
00042 #ifdef HAVE_IFPACK_AZTECOO
00043 #include "Ifpack_DiagPreconditioner.h"
00044 #include "AztecOO.h"
00045 #endif
00046 
00047 #ifdef HAVE_IFPACK_EPETRAEXT
00048 #include "Epetra_CrsMatrix.h"
00049 #include "EpetraExt_PointToBlockDiagPermute.h"
00050 #endif
00051 
00052 
00053 #define ABS(x) ((x)>0?(x):-(x))
00054 
00055 // Helper function for normal equations
00056 inline void Apply_Transpose(Teuchos::RCP<const Epetra_Operator> Operator_,const Epetra_MultiVector &X,Epetra_MultiVector &Y){
00057   Epetra_Operator * Operator=const_cast<Epetra_Operator*>(&*Operator_);
00058   Operator->SetUseTranspose(true);
00059   Operator->Apply(X,Y);
00060   Operator->SetUseTranspose(false);
00061 }
00062 
00063 
00064 
00065 
00066 //==============================================================================
00067 // NOTE: any change to the default values should be committed to the other
00068 //       constructor as well.
00069 Ifpack_Chebyshev::
00070 Ifpack_Chebyshev(const Epetra_Operator* Operator) :
00071   IsInitialized_(false),
00072   IsComputed_(false),
00073   NumInitialize_(0),
00074   NumCompute_(0),
00075   NumApplyInverse_(0),
00076   InitializeTime_(0.0),
00077   ComputeTime_(0.0),
00078   ApplyInverseTime_(0.0),
00079   ComputeFlops_(0.0),
00080   ApplyInverseFlops_(0.0),
00081   PolyDegree_(1),
00082   UseTranspose_(false),
00083   Condest_(-1.0),
00084   ComputeCondest_(false),
00085   EigRatio_(30.0),
00086   Label_(),
00087   LambdaMin_(0.0),
00088   LambdaMax_(100.0),
00089   MinDiagonalValue_(0.0),
00090   NumMyRows_(0),
00091   NumMyNonzeros_(0),
00092   NumGlobalRows_(0),
00093   NumGlobalNonzeros_(0),
00094   Operator_(Teuchos::rcp(Operator,false)),
00095   UseBlockMode_(false),  
00096   SolveNormalEquations_(false),
00097   IsRowMatrix_(false), 
00098   ZeroStartingSolution_(true)
00099 {
00100 }
00101 
00102 //==============================================================================
00103 // NOTE: This constructor has been introduced because SWIG does not appear
00104 //       to appreciate dynamic_cast. An instruction of type
00105 //       Matrix_ = dynamic_cast<const Epetra_RowMatrix*> in the
00106 //       other construction does not work in PyTrilinos -- of course
00107 //       it does in any C++ code (for an Epetra_Operator that is also
00108 //       an Epetra_RowMatrix).
00109 //
00110 // FIXME: move declarations into a separate method?
00111 Ifpack_Chebyshev::
00112 Ifpack_Chebyshev(const Epetra_RowMatrix* Operator) :
00113   IsInitialized_(false),
00114   IsComputed_(false),
00115   NumInitialize_(0),
00116   NumCompute_(0),
00117   NumApplyInverse_(0),
00118   InitializeTime_(0.0),
00119   ComputeTime_(0.0),
00120   ApplyInverseTime_(0.0),
00121   ComputeFlops_(0.0),
00122   ApplyInverseFlops_(0.0),
00123   PolyDegree_(1),
00124   UseTranspose_(false),
00125   Condest_(-1.0),
00126   ComputeCondest_(false),
00127   EigRatio_(30.0),
00128   EigMaxIters_(10),
00129   Label_(),
00130   LambdaMin_(0.0),
00131   LambdaMax_(100.0),
00132   MinDiagonalValue_(0.0),
00133   NumMyRows_(0),
00134   NumMyNonzeros_(0),
00135   NumGlobalRows_(0),
00136   NumGlobalNonzeros_(0),
00137   Operator_(Teuchos::rcp(Operator,false)),
00138   Matrix_(Teuchos::rcp(Operator,false)),
00139   UseBlockMode_(false),
00140   SolveNormalEquations_(false),
00141   IsRowMatrix_(true), 
00142   ZeroStartingSolution_(true)
00143 {
00144 }
00145 
00146 //==============================================================================
00147 int Ifpack_Chebyshev::SetParameters(Teuchos::ParameterList& List)
00148 {
00149 
00150   EigRatio_             = List.get("chebyshev: ratio eigenvalue", EigRatio_);
00151   LambdaMin_            = List.get("chebyshev: min eigenvalue", LambdaMin_);
00152   LambdaMax_            = List.get("chebyshev: max eigenvalue", LambdaMax_);
00153   PolyDegree_           = List.get("chebyshev: degree",PolyDegree_);
00154   MinDiagonalValue_     = List.get("chebyshev: min diagonal value", 
00155                                    MinDiagonalValue_);
00156   ZeroStartingSolution_ = List.get("chebyshev: zero starting solution", 
00157                                    ZeroStartingSolution_);
00158 
00159   Epetra_Vector* ID     = List.get("chebyshev: operator inv diagonal", 
00160                                    (Epetra_Vector*)0);
00161   EigMaxIters_          = List.get("chebyshev: eigenvalue max iterations",EigMaxIters_);
00162 
00163 #ifdef HAVE_IFPACK_EPETRAEXT
00164   // This is *always* false if EpetraExt isn't enabled
00165   UseBlockMode_         = List.get("chebyshev: use block mode",UseBlockMode_);
00166   if(!List.isParameter("chebyshev: block list")) UseBlockMode_=false;
00167   else{
00168     BlockList_          = List.get("chebyshev: block list",BlockList_);
00169 
00170     // Since we know we're doing a matrix inverse, clobber the block list
00171     // w/"invert" if it's set to multiply      
00172     Teuchos::ParameterList Blist;
00173     Blist=BlockList_.get("blockdiagmatrix: list",Blist);
00174     string dummy("invert");
00175     string ApplyMode=Blist.get("apply mode",dummy);
00176     if(ApplyMode==string("multiply")){
00177       Blist.set("apply mode","invert");
00178       BlockList_.set("blockdiagmatrix: list",Blist);
00179     }    
00180   }  
00181 #endif
00182 
00183   SolveNormalEquations_ = List.get("chebyshev: solve normal equations",SolveNormalEquations_);
00184 
00185   if (ID != 0) 
00186   {
00187     InvDiagonal_ = Teuchos::rcp( new Epetra_Vector(*ID) );
00188   }
00189 
00190   SetLabel();
00191 
00192   return(0);
00193 }
00194 
00195 //==============================================================================
00196 const Epetra_Comm& Ifpack_Chebyshev::Comm() const
00197 {
00198   return(Operator_->Comm());
00199 }
00200 
00201 //==============================================================================
00202 const Epetra_Map& Ifpack_Chebyshev::OperatorDomainMap() const
00203 {
00204   return(Operator_->OperatorDomainMap());
00205 }
00206 
00207 //==============================================================================
00208 const Epetra_Map& Ifpack_Chebyshev::OperatorRangeMap() const
00209 {
00210   return(Operator_->OperatorRangeMap());
00211 }
00212 
00213 //==============================================================================
00214 int Ifpack_Chebyshev::
00215 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00216 {
00217   if (IsComputed() == false)
00218     IFPACK_CHK_ERR(-3);
00219 
00220   if (X.NumVectors() != Y.NumVectors())
00221     IFPACK_CHK_ERR(-2);
00222 
00223   if (IsRowMatrix_)
00224   {
00225     IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
00226   }
00227   else
00228   {
00229     IFPACK_CHK_ERR(Operator_->Apply(X,Y));
00230   }
00231 
00232   return(0);
00233 }
00234 
00235 //==============================================================================
00236 int Ifpack_Chebyshev::Initialize()
00237 {
00238   IsInitialized_ = false;
00239 
00240   if (Operator_ == Teuchos::null)
00241     IFPACK_CHK_ERR(-2);
00242 
00243   if (Time_ == Teuchos::null)
00244     Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
00245 
00246   if (IsRowMatrix_)
00247   {
00248     if (Matrix().NumGlobalRows() != Matrix().NumGlobalCols())
00249       IFPACK_CHK_ERR(-2); // only square matrices
00250 
00251     NumMyRows_ = Matrix_->NumMyRows();
00252     NumMyNonzeros_ = Matrix_->NumMyNonzeros();
00253     NumGlobalRows_ = Matrix_->NumGlobalRows();
00254     NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros();
00255   }
00256   else
00257   {
00258     if (Operator_->OperatorDomainMap().NumGlobalElements() !=       
00259         Operator_->OperatorRangeMap().NumGlobalElements())
00260       IFPACK_CHK_ERR(-2); // only square operators
00261   }
00262 
00263   ++NumInitialize_;
00264   InitializeTime_ += Time_->ElapsedTime();
00265   IsInitialized_ = true;
00266   return(0);
00267 }
00268 
00269 //==============================================================================
00270 int Ifpack_Chebyshev::Compute()
00271 {
00272   if (!IsInitialized())
00273     IFPACK_CHK_ERR(Initialize());
00274 
00275   Time_->ResetStartTime();
00276 
00277   // reset values
00278   IsComputed_ = false;
00279   Condest_ = -1.0;
00280 
00281   if (PolyDegree_ <= 0)
00282     IFPACK_CHK_ERR(-2); // at least one application
00283 
00284 #ifdef HAVE_IFPACK_EPETRAEXT
00285   // Check to see if we can run in block mode
00286   if(IsRowMatrix_ && InvDiagonal_ == Teuchos::null && UseBlockMode_){
00287     const Epetra_CrsMatrix *CrsMatrix=dynamic_cast<const Epetra_CrsMatrix*>(&*Matrix_);
00288     
00289     // If we don't have CrsMatrix, we can't use the block preconditioner
00290     if(!CrsMatrix) UseBlockMode_=false;    
00291     else{
00292       int ierr;
00293       InvBlockDiagonal_=Teuchos::rcp(new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
00294       if(InvBlockDiagonal_==Teuchos::null) IFPACK_CHK_ERR(-6);
00295 
00296       ierr=InvBlockDiagonal_->SetParameters(BlockList_);
00297       if(ierr) IFPACK_CHK_ERR(-7);
00298 
00299       ierr=InvBlockDiagonal_->Compute();
00300       if(ierr) IFPACK_CHK_ERR(-8);
00301     }
00302 
00303     // Automatically Compute Eigenvalues
00304     double lambda_max=0;
00305     PowerMethod(EigMaxIters_,lambda_max);
00306     LambdaMax_=lambda_max;
00307     // Test for Exact Preconditioned case
00308     if(ABS(LambdaMax_-1) < 1e-6) LambdaMax_=LambdaMin_=1.0;
00309     else LambdaMin_=LambdaMax_/EigRatio_;
00310   }    
00311 #endif
00312   
00313   if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_)
00314   {
00315     InvDiagonal_ = Teuchos::rcp( new Epetra_Vector(Matrix().Map()) );
00316 
00317     if (InvDiagonal_ == Teuchos::null)
00318       IFPACK_CHK_ERR(-5);
00319 
00320     IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*InvDiagonal_));
00321 
00322     // Inverse diagonal elements
00323     // Replace zeros with 1.0
00324     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00325       double diag = (*InvDiagonal_)[i];
00326       if (IFPACK_ABS(diag) < MinDiagonalValue_)
00327         (*InvDiagonal_)[i] = MinDiagonalValue_;
00328       else
00329         (*InvDiagonal_)[i] = 1.0 / diag;
00330     }
00331   }
00332   // otherwise the inverse of the diagonal has been given by the user
00333 #ifdef IFPACK_FLOPCOUNTERS
00334   ComputeFlops_ += NumMyRows_;
00335 #endif
00336 
00337   ++NumCompute_;
00338   ComputeTime_ += Time_->ElapsedTime();
00339   IsComputed_ = true;
00340 
00341   return(0);
00342 }
00343 
00344 //==============================================================================
00345 ostream& Ifpack_Chebyshev::Print(ostream & os) const
00346 {
00347 
00348   double MyMinVal, MyMaxVal;
00349   double MinVal, MaxVal;
00350 
00351   if (IsComputed_) {
00352     InvDiagonal_->MinValue(&MyMinVal);
00353     InvDiagonal_->MaxValue(&MyMaxVal);
00354     Comm().MinAll(&MyMinVal,&MinVal,1);
00355     Comm().MinAll(&MyMaxVal,&MaxVal,1);
00356   }
00357 
00358   if (!Comm().MyPID()) {
00359     os << endl;
00360     os << "================================================================================" << endl;
00361     os << "Ifpack_Chebyshev" << endl;
00362     os << "Degree of polynomial      = " << PolyDegree_ << endl;
00363     os << "Condition number estimate = " << Condest() << endl;
00364     os << "Global number of rows     = " << Operator_->OperatorRangeMap().NumGlobalElements() << endl;
00365     if (IsComputed_) {
00366       os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
00367       os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
00368     }
00369     if (ZeroStartingSolution_) 
00370       os << "Using zero starting solution" << endl;
00371     else
00372       os << "Using input starting solution" << endl;
00373     os << endl;
00374     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00375     os << "-----           -------   --------------       ------------     --------" << endl;
00376     os << "Initialize()    "   << std::setw(5) << NumInitialize_ 
00377        << "  " << std::setw(15) << InitializeTime_ 
00378        << "              0.0              0.0" << endl;
00379     os << "Compute()       "   << std::setw(5) << NumCompute_ 
00380        << "  " << std::setw(15) << ComputeTime_
00381        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_;
00382     if (ComputeTime_ != 0.0)
00383       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
00384     else
00385       os << "  " << std::setw(15) << 0.0 << endl;
00386     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse_ 
00387        << "  " << std::setw(15) << ApplyInverseTime_
00388        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
00389     if (ApplyInverseTime_ != 0.0)
00390       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
00391     else
00392       os << "  " << std::setw(15) << 0.0 << endl;
00393     os << "================================================================================" << endl;
00394     os << endl;
00395   }
00396 
00397   return(os);
00398 }
00399 
00400 //==============================================================================
00401 double Ifpack_Chebyshev::
00402 Condest(const Ifpack_CondestType CT, 
00403         const int MaxIters, const double Tol,
00404     Epetra_RowMatrix* Matrix_in)
00405 {
00406   if (!IsComputed()) // cannot compute right now
00407     return(-1.0);
00408 
00409   // always computes it. Call Condest() with no parameters to get
00410   // the previous estimate.
00411   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00412 
00413   return(Condest_);
00414 }
00415 
00416 //==============================================================================
00417 void Ifpack_Chebyshev::SetLabel()
00418 {
00419   Label_ = "IFPACK (Chebyshev polynomial), degree=" + Ifpack_toString(PolyDegree_);
00420 }
00421 
00422 //==============================================================================
00423 int Ifpack_Chebyshev::
00424 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00425 {
00426   
00427   if (!IsComputed())
00428     IFPACK_CHK_ERR(-3);
00429 
00430   if (PolyDegree_ == 0)
00431     return 0;
00432 
00433   int nVec = X.NumVectors();
00434   int len = X.MyLength();
00435   if (nVec != Y.NumVectors())
00436     IFPACK_CHK_ERR(-2);
00437 
00438   Time_->ResetStartTime();
00439 
00440   // AztecOO gives X and Y pointing to the same memory location,
00441   // need to create an auxiliary vector, Xcopy
00442   Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
00443   if (X.Pointers()[0] == Y.Pointers()[0])
00444     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00445   else
00446     Xcopy = Teuchos::rcp( &X, false );
00447 
00448   double **xPtr = 0, **yPtr = 0;
00449   Xcopy->ExtractView(&xPtr);
00450   Y.ExtractView(&yPtr);
00451 
00452 #ifdef HAVE_IFPACK_EPETRAEXT
00453   EpetraExt_PointToBlockDiagPermute* IBD=0;
00454   if (UseBlockMode_) IBD=&*InvBlockDiagonal_;
00455 #endif
00456   
00457 
00458   //--- Do a quick solve when the matrix is identity
00459   double *invDiag=0;
00460   if(!UseBlockMode_) invDiag=InvDiagonal_->Values();
00461   if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
00462 #ifdef HAVE_IFPACK_EPETRAEXT
00463     if(UseBlockMode_) IBD->ApplyInverse(*Xcopy,Y);
00464     else
00465 #endif
00466     if (nVec == 1) {
00467       double *yPointer = yPtr[0], *xPointer = xPtr[0];
00468       for (int i = 0; i < len; ++i)
00469         yPointer[i] = xPointer[i]*invDiag[i];
00470     }
00471     else {
00472       int i, k;
00473       for (i = 0; i < len; ++i) {
00474         double coeff = invDiag[i];
00475         for (k = 0; k < nVec; ++k)
00476           yPtr[k][i] = xPtr[k][i] * coeff;
00477       }
00478     } // if (nVec == 1)
00479     return 0;
00480   } // if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_))
00481 
00482   //--- Initialize coefficients
00483   // Note that delta stores the inverse of ML_Cheby::delta
00484   double alpha = LambdaMax_ / EigRatio_;
00485   double beta = 1.1 * LambdaMax_;
00486   double delta = 2.0 / (beta - alpha);
00487   double theta = 0.5 * (beta + alpha);
00488   double s1 = theta * delta;
00489 
00490   //--- Define vectors
00491   // In ML_Cheby, V corresponds to pAux and W to dk
00492   Epetra_MultiVector V(X);
00493   Epetra_MultiVector W(X);
00494 #ifdef HAVE_IFPACK_EPETRAEXT
00495   Epetra_MultiVector Temp(X);
00496 #endif
00497   
00498   double *vPointer = V.Values(), *wPointer = W.Values();
00499 
00500   double oneOverTheta = 1.0/theta;
00501   int i, j, k;
00502 
00503 
00504   //--- If solving normal equations, multiply RHS by A^T
00505   if(SolveNormalEquations_){
00506     Apply_Transpose(Operator_,Y,V);
00507     Y=V;
00508   }
00509 
00510   // Do the smoothing when block scaling is turned OFF
00511   // --- Treat the initial guess
00512   if (ZeroStartingSolution_ == false) {
00513     Operator_->Apply(Y, V);
00514     // Compute W = invDiag * ( X - V )/ Theta
00515 #ifdef HAVE_IFPACK_EPETRAEXT    
00516     if(UseBlockMode_) {
00517       Temp.Update(oneOverTheta,X,-oneOverTheta,V,0.0);
00518       IBD->ApplyInverse(Temp,W);
00519 
00520       // Perform additional matvecs for normal equations
00521       // CMS: Testing this only in block mode FOR NOW
00522       if(SolveNormalEquations_){
00523     IBD->ApplyInverse(W,Temp);
00524     Apply_Transpose(Operator_,Temp,W);
00525       }
00526     }
00527     else
00528 #endif
00529     if (nVec == 1) {
00530       double *xPointer = xPtr[0];
00531       for (i = 0; i < len; ++i)
00532         wPointer[i] = invDiag[i] * (xPointer[i] - vPointer[i]) * oneOverTheta;
00533     }
00534     else {
00535       for (i = 0; i < len; ++i) {
00536         double coeff = invDiag[i]*oneOverTheta;
00537         double *wi = wPointer + i, *vi = vPointer + i;
00538         for (k = 0; k < nVec; ++k) {
00539           *wi = (xPtr[k][i] - (*vi)) * coeff;
00540           wi = wi + len; vi = vi + len;
00541         }
00542       }
00543     } // if (nVec == 1)
00544     // Update the vector Y
00545     Y.Update(1.0, W, 1.0);
00546   }
00547   else {
00548     // Compute W = invDiag * X / Theta
00549 #ifdef HAVE_IFPACK_EPETRAEXT    
00550     if(UseBlockMode_) {
00551       IBD->ApplyInverse(X,W);
00552 
00553       // Perform additional matvecs for normal equations
00554       // CMS: Testing this only in block mode FOR NOW
00555       if(SolveNormalEquations_){
00556     IBD->ApplyInverse(W,Temp);
00557     Apply_Transpose(Operator_,Temp,W);
00558       }
00559 
00560       W.Scale(oneOverTheta);
00561       Y.Update(1.0, W, 0.0);      
00562     }
00563     else
00564 #endif
00565     if (nVec == 1) {
00566       double *xPointer = xPtr[0];
00567       for (i = 0; i < len; ++i){
00568         wPointer[i] = invDiag[i] * xPointer[i] * oneOverTheta;
00569       }
00570       memcpy(yPtr[0], wPointer, len*sizeof(double));
00571     }
00572     else {
00573       for (i = 0; i < len; ++i) {
00574         double coeff = invDiag[i]*oneOverTheta;
00575         double *wi = wPointer + i;
00576         for (k = 0; k < nVec; ++k) {
00577           *wi = xPtr[k][i] * coeff;
00578           wi = wi + len;
00579         }
00580       }
00581       for (k = 0; k < nVec; ++k)
00582         memcpy(yPtr[k], wPointer + k*len, len*sizeof(double));
00583     } // if (nVec == 1)
00584   } // if (ZeroStartingSolution_ == false)
00585   
00586   //--- Apply the polynomial
00587   double rhok = 1.0/s1, rhokp1;
00588   double dtemp1, dtemp2;
00589   int degreeMinusOne = PolyDegree_ - 1;
00590   if (nVec == 1) {
00591     double *xPointer = xPtr[0];
00592     for (k = 0; k < degreeMinusOne; ++k) {
00593       Operator_->Apply(Y, V);
00594       rhokp1 = 1.0 / (2.0*s1 - rhok);
00595       dtemp1 = rhokp1 * rhok;
00596       dtemp2 = 2.0 * rhokp1 * delta;
00597       rhok = rhokp1;
00598       // Compute W = dtemp1 * W
00599       W.Scale(dtemp1);
00600       // Compute W = W + dtemp2 * invDiag * ( X - V )
00601 #ifdef HAVE_IFPACK_EPETRAEXT    
00602     if(UseBlockMode_) {
00603       //NTS: We can clobber V since it will be reset in the Apply
00604       V.Update(dtemp2,X,-dtemp2);
00605       IBD->ApplyInverse(V,Temp);
00606 
00607       // Perform additional matvecs for normal equations
00608       // CMS: Testing this only in block mode FOR NOW
00609       if(SolveNormalEquations_){
00610     IBD->ApplyInverse(V,Temp);
00611     Apply_Transpose(Operator_,Temp,V);
00612       }
00613 
00614       W.Update(1.0,Temp,1.0);
00615     }
00616     else{
00617 #endif
00618       for (i = 0; i < len; ++i)
00619         wPointer[i] += dtemp2* invDiag[i] * (xPointer[i] - vPointer[i]);
00620 #ifdef HAVE_IFPACK_EPETRAEXT
00621     }
00622 #endif
00623 
00624       // Update the vector Y
00625       Y.Update(1.0, W, 1.0);
00626     } // for (k = 0; k < degreeMinusOne; ++k)
00627   }
00628   else {
00629     for (k = 0; k < degreeMinusOne; ++k) {
00630       Operator_->Apply(Y, V);
00631       rhokp1 = 1.0 / (2.0*s1 - rhok);
00632       dtemp1 = rhokp1 * rhok;
00633       dtemp2 = 2.0 * rhokp1 * delta;
00634       rhok = rhokp1;
00635       // Compute W = dtemp1 * W
00636       W.Scale(dtemp1);
00637       // Compute W = W + dtemp2 * invDiag * ( X - V )
00638 #ifdef HAVE_IFPACK_EPETRAEXT    
00639     if(UseBlockMode_) {
00640       //We can clobber V since it will be reset in the Apply
00641       V.Update(dtemp2,X,-dtemp2);
00642       IBD->ApplyInverse(V,Temp);
00643 
00644       // Perform additional matvecs for normal equations
00645       // CMS: Testing this only in block mode FOR NOW
00646       if(SolveNormalEquations_){
00647     IBD->ApplyInverse(V,Temp);
00648     Apply_Transpose(Operator_,Temp,V);
00649       }
00650 
00651 
00652       W.Update(1.0,Temp,1.0);
00653     }
00654     else{
00655 #endif
00656       for (i = 0; i < len; ++i) {
00657         double coeff = invDiag[i]*dtemp2;
00658         double *wi = wPointer + i, *vi = vPointer + i;
00659         for (j = 0; j < nVec; ++j) {
00660           *wi += (xPtr[j][i] - (*vi)) * coeff;
00661           wi = wi + len; vi = vi + len;
00662         }
00663       }
00664 #ifdef HAVE_IFPACK_EPETRAEXT
00665     }
00666 #endif      
00667       // Update the vector Y
00668       Y.Update(1.0, W, 1.0);
00669     } // for (k = 0; k < degreeMinusOne; ++k)
00670   } // if (nVec == 1)
00671 
00672   
00673   // Flops are updated in each of the following. 
00674   ++NumApplyInverse_;
00675   ApplyInverseTime_ += Time_->ElapsedTime();
00676   return(0);
00677 }
00678 
00679 //==============================================================================
00680 int Ifpack_Chebyshev::
00681 PowerMethod(const Epetra_Operator& Operator, 
00682             const Epetra_Vector& InvPointDiagonal, 
00683             const int MaximumIterations, 
00684             double& lambda_max)
00685 {
00686   // this is a simple power method
00687   lambda_max = 0.0;
00688   double RQ_top, RQ_bottom, norm;
00689   Epetra_Vector x(Operator.OperatorDomainMap());
00690   Epetra_Vector y(Operator.OperatorRangeMap());
00691   x.Random();
00692   x.Norm2(&norm);
00693   if (norm == 0.0) IFPACK_CHK_ERR(-1);
00694 
00695   x.Scale(1.0 / norm);
00696 
00697   for (int iter = 0; iter < MaximumIterations; ++iter)
00698   {
00699     Operator.Apply(x, y);
00700     IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
00701     IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
00702     IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
00703     lambda_max = RQ_top / RQ_bottom;
00704     IFPACK_CHK_ERR(y.Norm2(&norm));
00705     if (norm == 0.0) IFPACK_CHK_ERR(-1);
00706     IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
00707   }
00708 
00709   return(0);
00710 }
00711 
00712 //==============================================================================
00713 int Ifpack_Chebyshev::
00714 CG(const Epetra_Operator& Operator, 
00715    const Epetra_Vector& InvPointDiagonal, 
00716    const int MaximumIterations, 
00717    double& lambda_min, double& lambda_max)
00718 {
00719 #ifdef HAVE_IFPACK_AZTECOO
00720   Epetra_Vector x(Operator.OperatorDomainMap());
00721   Epetra_Vector y(Operator.OperatorRangeMap());
00722   x.Random();
00723   y.PutScalar(0.0);
00724 
00725   Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
00726   AztecOO solver(LP);
00727   solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
00728   solver.SetAztecOption(AZ_output, AZ_none);
00729 
00730   Ifpack_DiagPreconditioner diag(Operator.OperatorDomainMap(),
00731                                  Operator.OperatorRangeMap(),
00732                                  InvPointDiagonal);
00733   solver.SetPrecOperator(&diag);
00734   solver.Iterate(MaximumIterations, 1e-10);
00735 
00736   const double* status = solver.GetAztecStatus();
00737 
00738   lambda_min = status[AZ_lambda_min];
00739   lambda_max = status[AZ_lambda_max];
00740 
00741   return(0);
00742 #else
00743   cout << "You need to configure IFPACK with support for AztecOO" << endl;
00744   cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
00745   cout << "in your configure script." << endl;
00746   IFPACK_CHK_ERR(-1);
00747 #endif
00748 }
00749 
00750 //==============================================================================
00751 #ifdef HAVE_IFPACK_EPETRAEXT
00752 int Ifpack_Chebyshev::
00753 PowerMethod(const int MaximumIterations,  double& lambda_max)
00754 {
00755 
00756   if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
00757   // this is a simple power method
00758   lambda_max = 0.0;
00759   double RQ_top, RQ_bottom, norm;
00760   Epetra_Vector x(Operator_->OperatorDomainMap());
00761   Epetra_Vector y(Operator_->OperatorRangeMap());
00762   Epetra_Vector z(Operator_->OperatorRangeMap());
00763   x.Random();
00764   x.Norm2(&norm);
00765   if (norm == 0.0) IFPACK_CHK_ERR(-1);
00766 
00767   x.Scale(1.0 / norm);
00768 
00769   for (int iter = 0; iter < MaximumIterations; ++iter)
00770   {
00771     Operator_->Apply(x, z);
00772     InvBlockDiagonal_->ApplyInverse(z,y);
00773     if(SolveNormalEquations_){
00774       InvBlockDiagonal_->ApplyInverse(y,z);
00775       Apply_Transpose(Operator_,z, y);
00776     }
00777 
00778     IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
00779     IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
00780     lambda_max = RQ_top / RQ_bottom;
00781     IFPACK_CHK_ERR(y.Norm2(&norm));
00782     if (norm == 0.0) IFPACK_CHK_ERR(-1);
00783     IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
00784   }
00785 
00786   return(0);
00787 }
00788 #endif
00789 
00790 //==============================================================================
00791 #ifdef HAVE_IFPACK_EPETRAEXT
00792 int Ifpack_Chebyshev::
00793 CG(const int MaximumIterations, 
00794    double& lambda_min, double& lambda_max)
00795 {
00796   IFPACK_CHK_ERR(-1);// NTS: This always seems to yield errors in AztecOO, ergo,
00797                      // I turned it off.
00798 
00799   if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
00800   
00801 #ifdef HAVE_IFPACK_AZTECOO
00802   Epetra_Vector x(Operator_->OperatorDomainMap());
00803   Epetra_Vector y(Operator_->OperatorRangeMap());
00804   x.Random();
00805   y.PutScalar(0.0);
00806   Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
00807   
00808   AztecOO solver(LP);
00809   solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
00810   solver.SetAztecOption(AZ_output, AZ_none);
00811 
00812   solver.SetPrecOperator(&*InvBlockDiagonal_);
00813   solver.Iterate(MaximumIterations, 1e-10);
00814 
00815   const double* status = solver.GetAztecStatus();
00816 
00817   lambda_min = status[AZ_lambda_min];
00818   lambda_max = status[AZ_lambda_max];
00819   
00820   return(0);
00821 #else
00822   cout << "You need to configure IFPACK with support for AztecOO" << endl;
00823   cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
00824   cout << "in your configure script." << endl;
00825   IFPACK_CHK_ERR(-1);
00826 #endif
00827 }
00828 #endif
 All Classes Files Functions Variables Enumerations Friends