Ifpack_Chebyshev.cpp

Go to the documentation of this file.
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 
00334   ComputeFlops_ += NumMyRows_;
00335 
00336   ++NumCompute_;
00337   ComputeTime_ += Time_->ElapsedTime();
00338   IsComputed_ = true;
00339 
00340   return(0);
00341 }
00342 
00343 //==============================================================================
00344 ostream& Ifpack_Chebyshev::Print(ostream & os) const
00345 {
00346 
00347   double MyMinVal, MyMaxVal;
00348   double MinVal, MaxVal;
00349 
00350   if (IsComputed_) {
00351     InvDiagonal_->MinValue(&MyMinVal);
00352     InvDiagonal_->MaxValue(&MyMaxVal);
00353     Comm().MinAll(&MyMinVal,&MinVal,1);
00354     Comm().MinAll(&MyMaxVal,&MaxVal,1);
00355   }
00356 
00357   if (!Comm().MyPID()) {
00358     os << endl;
00359     os << "================================================================================" << endl;
00360     os << "Ifpack_Chebyshev" << endl;
00361     os << "Degree of polynomial      = " << PolyDegree_ << endl;
00362     os << "Condition number estimate = " << Condest() << endl;
00363     os << "Global number of rows     = " << Operator_->OperatorRangeMap().NumGlobalElements() << endl;
00364     if (IsComputed_) {
00365       os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
00366       os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
00367     }
00368     if (ZeroStartingSolution_) 
00369       os << "Using zero starting solution" << endl;
00370     else
00371       os << "Using input starting solution" << endl;
00372     os << endl;
00373     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00374     os << "-----           -------   --------------       ------------     --------" << endl;
00375     os << "Initialize()    "   << std::setw(5) << NumInitialize_ 
00376        << "  " << std::setw(15) << InitializeTime_ 
00377        << "              0.0              0.0" << endl;
00378     os << "Compute()       "   << std::setw(5) << NumCompute_ 
00379        << "  " << std::setw(15) << ComputeTime_
00380        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_;
00381     if (ComputeTime_ != 0.0)
00382       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
00383     else
00384       os << "  " << std::setw(15) << 0.0 << endl;
00385     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse_ 
00386        << "  " << std::setw(15) << ApplyInverseTime_
00387        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
00388     if (ApplyInverseTime_ != 0.0)
00389       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
00390     else
00391       os << "  " << std::setw(15) << 0.0 << endl;
00392     os << "================================================================================" << endl;
00393     os << endl;
00394   }
00395 
00396   return(os);
00397 }
00398 
00399 //==============================================================================
00400 double Ifpack_Chebyshev::
00401 Condest(const Ifpack_CondestType CT, 
00402         const int MaxIters, const double Tol,
00403   Epetra_RowMatrix* Matrix_in)
00404 {
00405   if (!IsComputed()) // cannot compute right now
00406     return(-1.0);
00407 
00408   // always computes it. Call Condest() with no parameters to get
00409   // the previous estimate.
00410   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00411 
00412   return(Condest_);
00413 }
00414 
00415 //==============================================================================
00416 void Ifpack_Chebyshev::SetLabel()
00417 {
00418   Label_ = "IFPACK (Chebyshev polynomial), degree=" + Ifpack_toString(PolyDegree_);
00419 }
00420 
00421 //==============================================================================
00422 int Ifpack_Chebyshev::
00423 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00424 {
00425   
00426   if (!IsComputed())
00427     IFPACK_CHK_ERR(-3);
00428 
00429   if (PolyDegree_ == 0)
00430     return 0;
00431 
00432   int nVec = X.NumVectors();
00433   int len = X.MyLength();
00434   if (nVec != Y.NumVectors())
00435     IFPACK_CHK_ERR(-2);
00436 
00437   Time_->ResetStartTime();
00438 
00439   // AztecOO gives X and Y pointing to the same memory location,
00440   // need to create an auxiliary vector, Xcopy
00441   Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
00442   if (X.Pointers()[0] == Y.Pointers()[0])
00443     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00444   else
00445     Xcopy = Teuchos::rcp( &X, false );
00446 
00447   double **xPtr = 0, **yPtr = 0;
00448   Xcopy->ExtractView(&xPtr);
00449   Y.ExtractView(&yPtr);
00450 
00451 #ifdef HAVE_IFPACK_EPETRAEXT
00452   EpetraExt_PointToBlockDiagPermute* IBD=0;
00453   if (UseBlockMode_) IBD=&*InvBlockDiagonal_;
00454 #endif
00455   
00456 
00457   //--- Do a quick solve when the matrix is identity
00458   double *invDiag=0;
00459   if(!UseBlockMode_) invDiag=InvDiagonal_->Values();
00460   if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
00461 #ifdef HAVE_IFPACK_EPETRAEXT
00462     if(UseBlockMode_) IBD->ApplyInverse(*Xcopy,Y);
00463     else
00464 #endif
00465     if (nVec == 1) {
00466       double *yPointer = yPtr[0], *xPointer = xPtr[0];
00467       for (int i = 0; i < len; ++i)
00468         yPointer[i] = xPointer[i]*invDiag[i];
00469     }
00470     else {
00471       int i, k;
00472       for (i = 0; i < len; ++i) {
00473         double coeff = invDiag[i];
00474         for (k = 0; k < nVec; ++k)
00475           yPtr[k][i] = xPtr[k][i] * coeff;
00476       }
00477     } // if (nVec == 1)
00478     return 0;
00479   } // if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_))
00480 
00481   //--- Initialize coefficients
00482   // Note that delta stores the inverse of ML_Cheby::delta
00483   double alpha = LambdaMax_ / EigRatio_;
00484   double beta = 1.1 * LambdaMax_;
00485   double delta = 2.0 / (beta - alpha);
00486   double theta = 0.5 * (beta + alpha);
00487   double s1 = theta * delta;
00488 
00489   //--- Define vectors
00490   // In ML_Cheby, V corresponds to pAux and W to dk
00491   Epetra_MultiVector V(X);
00492   Epetra_MultiVector W(X);
00493 #ifdef HAVE_IFPACK_EPETRAEXT
00494   Epetra_MultiVector Temp(X);
00495 #endif
00496   
00497   double *vPointer = V.Values(), *wPointer = W.Values();
00498 
00499   double oneOverTheta = 1.0/theta;
00500   int i, j, k;
00501 
00502 
00503   //--- If solving normal equations, multiply RHS by A^T
00504   if(SolveNormalEquations_){
00505     Apply_Transpose(Operator_,Y,V);
00506     Y=V;
00507   }
00508 
00509   // Do the smoothing when block scaling is turned OFF
00510   // --- Treat the initial guess
00511   if (ZeroStartingSolution_ == false) {
00512     Operator_->Apply(Y, V);
00513     // Compute W = invDiag * ( X - V )/ Theta
00514 #ifdef HAVE_IFPACK_EPETRAEXT    
00515     if(UseBlockMode_) {
00516       Temp.Update(oneOverTheta,X,-oneOverTheta,V,0.0);
00517       IBD->ApplyInverse(Temp,W);
00518 
00519       // Perform additional matvecs for normal equations
00520       // CMS: Testing this only in block mode FOR NOW
00521       if(SolveNormalEquations_){
00522   IBD->ApplyInverse(W,Temp);
00523   Apply_Transpose(Operator_,Temp,W);
00524       }
00525     }
00526     else
00527 #endif
00528     if (nVec == 1) {
00529       double *xPointer = xPtr[0];
00530       for (i = 0; i < len; ++i)
00531         wPointer[i] = invDiag[i] * (xPointer[i] - vPointer[i]) * oneOverTheta;
00532     }
00533     else {
00534       for (i = 0; i < len; ++i) {
00535         double coeff = invDiag[i]*oneOverTheta;
00536         double *wi = wPointer + i, *vi = vPointer + i;
00537         for (k = 0; k < nVec; ++k) {
00538           *wi = (xPtr[k][i] - (*vi)) * coeff;
00539           wi = wi + len; vi = vi + len;
00540         }
00541       }
00542     } // if (nVec == 1)
00543     // Update the vector Y
00544     Y.Update(1.0, W, 1.0);
00545   }
00546   else {
00547     // Compute W = invDiag * X / Theta
00548 #ifdef HAVE_IFPACK_EPETRAEXT    
00549     if(UseBlockMode_) {
00550       IBD->ApplyInverse(X,W);
00551 
00552       // Perform additional matvecs for normal equations
00553       // CMS: Testing this only in block mode FOR NOW
00554       if(SolveNormalEquations_){
00555   IBD->ApplyInverse(W,Temp);
00556   Apply_Transpose(Operator_,Temp,W);
00557       }
00558 
00559       W.Scale(oneOverTheta);
00560       Y.Update(1.0, W, 0.0);      
00561     }
00562     else
00563 #endif
00564     if (nVec == 1) {
00565       double *xPointer = xPtr[0];
00566       for (i = 0; i < len; ++i){
00567         wPointer[i] = invDiag[i] * xPointer[i] * oneOverTheta;
00568       }
00569       memcpy(yPtr[0], wPointer, len*sizeof(double));
00570     }
00571     else {
00572       for (i = 0; i < len; ++i) {
00573         double coeff = invDiag[i]*oneOverTheta;
00574         double *wi = wPointer + i;
00575         for (k = 0; k < nVec; ++k) {
00576           *wi = xPtr[k][i] * coeff;
00577           wi = wi + len;
00578         }
00579       }
00580       for (k = 0; k < nVec; ++k)
00581         memcpy(yPtr[k], wPointer + k*len, len*sizeof(double));
00582     } // if (nVec == 1)
00583   } // if (ZeroStartingSolution_ == false)
00584   
00585   //--- Apply the polynomial
00586   double rhok = 1.0/s1, rhokp1;
00587   double dtemp1, dtemp2;
00588   int degreeMinusOne = PolyDegree_ - 1;
00589   if (nVec == 1) {
00590     double *xPointer = xPtr[0];
00591     for (k = 0; k < degreeMinusOne; ++k) {
00592       Operator_->Apply(Y, V);
00593       rhokp1 = 1.0 / (2.0*s1 - rhok);
00594       dtemp1 = rhokp1 * rhok;
00595       dtemp2 = 2.0 * rhokp1 * delta;
00596       rhok = rhokp1;
00597       // Compute W = dtemp1 * W
00598       W.Scale(dtemp1);
00599       // Compute W = W + dtemp2 * invDiag * ( X - V )
00600 #ifdef HAVE_IFPACK_EPETRAEXT    
00601     if(UseBlockMode_) {
00602       //NTS: We can clobber V since it will be reset in the Apply
00603       V.Update(dtemp2,X,-dtemp2);
00604       IBD->ApplyInverse(V,Temp);
00605 
00606       // Perform additional matvecs for normal equations
00607       // CMS: Testing this only in block mode FOR NOW
00608       if(SolveNormalEquations_){
00609   IBD->ApplyInverse(V,Temp);
00610   Apply_Transpose(Operator_,Temp,V);
00611       }
00612 
00613       W.Update(1.0,Temp,1.0);
00614     }
00615     else{
00616 #endif
00617       for (i = 0; i < len; ++i)
00618         wPointer[i] += dtemp2* invDiag[i] * (xPointer[i] - vPointer[i]);
00619 #ifdef HAVE_IFPACK_EPETRAEXT
00620     }
00621 #endif
00622 
00623       // Update the vector Y
00624       Y.Update(1.0, W, 1.0);
00625     } // for (k = 0; k < degreeMinusOne; ++k)
00626   }
00627   else {
00628     for (k = 0; k < degreeMinusOne; ++k) {
00629       Operator_->Apply(Y, V);
00630       rhokp1 = 1.0 / (2.0*s1 - rhok);
00631       dtemp1 = rhokp1 * rhok;
00632       dtemp2 = 2.0 * rhokp1 * delta;
00633       rhok = rhokp1;
00634       // Compute W = dtemp1 * W
00635       W.Scale(dtemp1);
00636       // Compute W = W + dtemp2 * invDiag * ( X - V )
00637 #ifdef HAVE_IFPACK_EPETRAEXT    
00638     if(UseBlockMode_) {
00639       //We can clobber V since it will be reset in the Apply
00640       V.Update(dtemp2,X,-dtemp2);
00641       IBD->ApplyInverse(V,Temp);
00642 
00643       // Perform additional matvecs for normal equations
00644       // CMS: Testing this only in block mode FOR NOW
00645       if(SolveNormalEquations_){
00646   IBD->ApplyInverse(V,Temp);
00647   Apply_Transpose(Operator_,Temp,V);
00648       }
00649 
00650 
00651       W.Update(1.0,Temp,1.0);
00652     }
00653     else{
00654 #endif
00655       for (i = 0; i < len; ++i) {
00656         double coeff = invDiag[i]*dtemp2;
00657         double *wi = wPointer + i, *vi = vPointer + i;
00658         for (j = 0; j < nVec; ++j) {
00659           *wi += (xPtr[j][i] - (*vi)) * coeff;
00660           wi = wi + len; vi = vi + len;
00661         }
00662       }
00663 #ifdef HAVE_IFPACK_EPETRAEXT
00664     }
00665 #endif      
00666       // Update the vector Y
00667       Y.Update(1.0, W, 1.0);
00668     } // for (k = 0; k < degreeMinusOne; ++k)
00669   } // if (nVec == 1)
00670 
00671   
00672   // Flops are updated in each of the following. 
00673   ++NumApplyInverse_;
00674   ApplyInverseTime_ += Time_->ElapsedTime();
00675   return(0);
00676 }
00677 
00678 //==============================================================================
00679 int Ifpack_Chebyshev::
00680 PowerMethod(const Epetra_Operator& Operator, 
00681             const Epetra_Vector& InvPointDiagonal, 
00682             const int MaximumIterations, 
00683             double& lambda_max)
00684 {
00685   // this is a simple power method
00686   lambda_max = 0.0;
00687   double RQ_top, RQ_bottom, norm;
00688   Epetra_Vector x(Operator.OperatorDomainMap());
00689   Epetra_Vector y(Operator.OperatorRangeMap());
00690   x.Random();
00691   x.Norm2(&norm);
00692   if (norm == 0.0) IFPACK_CHK_ERR(-1);
00693 
00694   x.Scale(1.0 / norm);
00695 
00696   for (int iter = 0; iter < MaximumIterations; ++iter)
00697   {
00698     Operator.Apply(x, y);
00699     IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
00700     IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
00701     IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
00702     lambda_max = RQ_top / RQ_bottom;
00703     IFPACK_CHK_ERR(y.Norm2(&norm));
00704     if (norm == 0.0) IFPACK_CHK_ERR(-1);
00705     IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
00706   }
00707 
00708   return(0);
00709 }
00710 
00711 //==============================================================================
00712 int Ifpack_Chebyshev::
00713 CG(const Epetra_Operator& Operator, 
00714    const Epetra_Vector& InvPointDiagonal, 
00715    const int MaximumIterations, 
00716    double& lambda_min, double& lambda_max)
00717 {
00718 #ifdef HAVE_IFPACK_AZTECOO
00719   Epetra_Vector x(Operator.OperatorDomainMap());
00720   Epetra_Vector y(Operator.OperatorRangeMap());
00721   x.Random();
00722   y.PutScalar(0.0);
00723 
00724   Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
00725   AztecOO solver(LP);
00726   solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
00727   solver.SetAztecOption(AZ_output, AZ_none);
00728 
00729   Ifpack_DiagPreconditioner diag(Operator.OperatorDomainMap(),
00730                                  Operator.OperatorRangeMap(),
00731                                  InvPointDiagonal);
00732   solver.SetPrecOperator(&diag);
00733   solver.Iterate(MaximumIterations, 1e-10);
00734 
00735   const double* status = solver.GetAztecStatus();
00736 
00737   lambda_min = status[AZ_lambda_min];
00738   lambda_max = status[AZ_lambda_max];
00739 
00740   return(0);
00741 #else
00742   cout << "You need to configure IFPACK with support for AztecOO" << endl;
00743   cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
00744   cout << "in your configure script." << endl;
00745   IFPACK_CHK_ERR(-1);
00746 #endif
00747 }
00748 
00749 //==============================================================================
00750 #ifdef HAVE_IFPACK_EPETRAEXT
00751 int Ifpack_Chebyshev::
00752 PowerMethod(const int MaximumIterations,  double& lambda_max)
00753 {
00754 
00755   if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
00756   // this is a simple power method
00757   lambda_max = 0.0;
00758   double RQ_top, RQ_bottom, norm;
00759   Epetra_Vector x(Operator_->OperatorDomainMap());
00760   Epetra_Vector y(Operator_->OperatorRangeMap());
00761   Epetra_Vector z(Operator_->OperatorRangeMap());
00762   x.Random();
00763   x.Norm2(&norm);
00764   if (norm == 0.0) IFPACK_CHK_ERR(-1);
00765 
00766   x.Scale(1.0 / norm);
00767 
00768   for (int iter = 0; iter < MaximumIterations; ++iter)
00769   {
00770     Operator_->Apply(x, z);
00771     InvBlockDiagonal_->ApplyInverse(z,y);
00772     if(SolveNormalEquations_){
00773       InvBlockDiagonal_->ApplyInverse(y,z);
00774       Apply_Transpose(Operator_,z, y);
00775     }
00776 
00777     IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
00778     IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
00779     lambda_max = RQ_top / RQ_bottom;
00780     IFPACK_CHK_ERR(y.Norm2(&norm));
00781     if (norm == 0.0) IFPACK_CHK_ERR(-1);
00782     IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
00783   }
00784 
00785   return(0);
00786 }
00787 #endif
00788 
00789 //==============================================================================
00790 #ifdef HAVE_IFPACK_EPETRAEXT
00791 int Ifpack_Chebyshev::
00792 CG(const int MaximumIterations, 
00793    double& lambda_min, double& lambda_max)
00794 {
00795   IFPACK_CHK_ERR(-1);// NTS: This always seems to yield errors in AztecOO, ergo,
00796                      // I turned it off.
00797 
00798   if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
00799   
00800 #ifdef HAVE_IFPACK_AZTECOO
00801   Epetra_Vector x(Operator_->OperatorDomainMap());
00802   Epetra_Vector y(Operator_->OperatorRangeMap());
00803   x.Random();
00804   y.PutScalar(0.0);
00805   Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
00806   
00807   AztecOO solver(LP);
00808   solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
00809   solver.SetAztecOption(AZ_output, AZ_none);
00810 
00811   solver.SetPrecOperator(&*InvBlockDiagonal_);
00812   solver.Iterate(MaximumIterations, 1e-10);
00813 
00814   const double* status = solver.GetAztecStatus();
00815 
00816   lambda_min = status[AZ_lambda_min];
00817   lambda_max = status[AZ_lambda_max];
00818   
00819   return(0);
00820 #else
00821   cout << "You need to configure IFPACK with support for AztecOO" << endl;
00822   cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
00823   cout << "in your configure script." << endl;
00824   IFPACK_CHK_ERR(-1);
00825 #endif
00826 }
00827 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:34 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3