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

Generated on Wed May 12 21:30:18 2010 for IFPACK by  doxygen 1.4.7