Ifpack_Chebyshev.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #ifdef HAVE_IFPACK_TEUCHOS
00003 #include <iomanip>
00004 #include "Epetra_Operator.h"
00005 #include "Epetra_RowMatrix.h"
00006 #include "Epetra_Comm.h"
00007 #include "Epetra_Map.h"
00008 #include "Epetra_MultiVector.h"
00009 #include "Epetra_Vector.h"
00010 #include "Epetra_Time.h"
00011 #include "Ifpack_Chebyshev.h"
00012 #include "Ifpack_Utils.h"
00013 #include "Ifpack_Condest.h"
00014 #ifdef HAVE_IFPACK_AZTECOO
00015 #include "Ifpack_DiagPreconditioner.h"
00016 #include "AztecOO.h"
00017 #endif
00018 
00019 //==============================================================================
00020 // NOTE: any change to the default values should be committed to the other
00021 //       constructor as well.
00022 Ifpack_Chebyshev::
00023 Ifpack_Chebyshev(const Epetra_Operator* Operator) :
00024   IsInitialized_(false),
00025   IsComputed_(false),
00026   NumInitialize_(0),
00027   NumCompute_(0),
00028   NumApplyInverse_(0),
00029   InitializeTime_(0.0),
00030   ComputeTime_(0.0),
00031   ApplyInverseTime_(0.0),
00032   ComputeFlops_(0.0),
00033   ApplyInverseFlops_(0.0),
00034   PolyDegree_(1),
00035   UseTranspose_(false),
00036   Condest_(-1.0),
00037   ComputeCondest_(false),
00038   EigRatio_(30.0),
00039   Label_(),
00040   LambdaMin_(0.0),
00041   LambdaMax_(100.0),
00042   MinDiagonalValue_(0.0),
00043   NumMyRows_(0),
00044   NumMyNonzeros_(0),
00045   NumGlobalRows_(0),
00046   NumGlobalNonzeros_(0),
00047   Operator_(Operator),
00048   Matrix_(0),
00049   InvDiagonal_(0),
00050   IsRowMatrix_(false), 
00051   Time_(0),
00052   ZeroStartingSolution_(true)
00053 {
00054 }
00055 
00056 //==============================================================================
00057 // NOTE: This constructor has been introduced because SWIG does not appear
00058 //       to appreciate dynamic_cast. An instruction of type
00059 //       Matrix_ = dynamic_cast<const Epetra_RowMatrix*> in the
00060 //       other construction does not work in PyTrilinos -- of course
00061 //       it does in any C++ code (for an Epetra_Operator that is also
00062 //       an Epetra_RowMatrix).
00063 //
00064 // FIXME: move declarations into a separate method?
00065 Ifpack_Chebyshev::
00066 Ifpack_Chebyshev(const Epetra_RowMatrix* Operator) :
00067   IsInitialized_(false),
00068   IsComputed_(false),
00069   NumInitialize_(0),
00070   NumCompute_(0),
00071   NumApplyInverse_(0),
00072   InitializeTime_(0.0),
00073   ComputeTime_(0.0),
00074   ApplyInverseTime_(0.0),
00075   ComputeFlops_(0.0),
00076   ApplyInverseFlops_(0.0),
00077   PolyDegree_(1),
00078   UseTranspose_(false),
00079   Condest_(-1.0),
00080   ComputeCondest_(false),
00081   EigRatio_(30.0),
00082   Label_(),
00083   LambdaMin_(0.0),
00084   LambdaMax_(100.0),
00085   MinDiagonalValue_(0.0),
00086   NumMyRows_(0),
00087   NumMyNonzeros_(0),
00088   NumGlobalRows_(0),
00089   NumGlobalNonzeros_(0),
00090   Operator_(Operator),
00091   Matrix_(Operator),
00092   InvDiagonal_(0),
00093   IsRowMatrix_(true), 
00094   Time_(0),
00095   ZeroStartingSolution_(true)
00096 {
00097 }
00098 
00099 //==============================================================================
00100 Ifpack_Chebyshev::~Ifpack_Chebyshev()
00101 {
00102   if (InvDiagonal_)
00103     delete InvDiagonal_;
00104   if (Time_)
00105     delete Time_;
00106 }
00107 
00108 //==============================================================================
00109 int Ifpack_Chebyshev::SetParameters(Teuchos::ParameterList& List)
00110 {
00111 
00112   EigRatio_             = List.get("chebyshev: ratio eigenvalue", EigRatio_);
00113   LambdaMin_            = List.get("chebyshev: min eigenvalue", LambdaMin_);
00114   LambdaMax_            = List.get("chebyshev: max eigenvalue", LambdaMax_);
00115   PolyDegree_           = List.get("chebyshev: degree",PolyDegree_);
00116   MinDiagonalValue_     = List.get("chebyshev: min diagonal value", 
00117                                    MinDiagonalValue_);
00118   ZeroStartingSolution_ = List.get("chebyshev: zero starting solution", 
00119                                    ZeroStartingSolution_);
00120 
00121   Epetra_Vector* ID     = List.get("chebyshev: operator inv diagonal", 
00122                                    (Epetra_Vector*)0);
00123 
00124   if (ID != 0) 
00125   {
00126     if (InvDiagonal_) delete InvDiagonal_;
00127     InvDiagonal_ = new Epetra_Vector(*ID);
00128   }
00129 
00130   SetLabel();
00131 
00132   return(0);
00133 }
00134 
00135 //==============================================================================
00136 const Epetra_Comm& Ifpack_Chebyshev::Comm() const
00137 {
00138   return(Operator_->Comm());
00139 }
00140 
00141 //==============================================================================
00142 const Epetra_Map& Ifpack_Chebyshev::OperatorDomainMap() const
00143 {
00144   return(Operator_->OperatorDomainMap());
00145 }
00146 
00147 //==============================================================================
00148 const Epetra_Map& Ifpack_Chebyshev::OperatorRangeMap() const
00149 {
00150   return(Operator_->OperatorRangeMap());
00151 }
00152 
00153 //==============================================================================
00154 int Ifpack_Chebyshev::
00155 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00156 {
00157   if (IsComputed() == false)
00158     IFPACK_CHK_ERR(-3);
00159 
00160   if (X.NumVectors() != Y.NumVectors())
00161     IFPACK_CHK_ERR(-2);
00162 
00163   if (IsRowMatrix_)
00164   {
00165     IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
00166   }
00167   else
00168   {
00169     IFPACK_CHK_ERR(Operator_->Apply(X,Y));
00170   }
00171 
00172   return(0);
00173 }
00174 
00175 //==============================================================================
00176 int Ifpack_Chebyshev::Initialize()
00177 {
00178   IsInitialized_ = false;
00179 
00180   if (Operator_ == 0)
00181     IFPACK_CHK_ERR(-2);
00182 
00183   if (Time_ == 0)
00184     Time_ = new Epetra_Time(Comm());
00185 
00186   if (IsRowMatrix_)
00187   {
00188     if (Matrix().NumGlobalRows() != Matrix().NumGlobalCols())
00189       IFPACK_CHK_ERR(-2); // only square matrices
00190 
00191     NumMyRows_ = Matrix_->NumMyRows();
00192     NumMyNonzeros_ = Matrix_->NumMyNonzeros();
00193     NumGlobalRows_ = Matrix_->NumGlobalRows();
00194     NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros();
00195   }
00196   else
00197   {
00198     if (Operator_->OperatorDomainMap().NumGlobalElements() !=       
00199         Operator_->OperatorRangeMap().NumGlobalElements())
00200       IFPACK_CHK_ERR(-2); // only square operators
00201   }
00202 
00203   ++NumInitialize_;
00204   InitializeTime_ += Time_->ElapsedTime();
00205   IsInitialized_ = true;
00206   return(0);
00207 }
00208 
00209 //==============================================================================
00210 int Ifpack_Chebyshev::Compute()
00211 {
00212   if (!IsInitialized())
00213     IFPACK_CHK_ERR(Initialize());
00214 
00215   Time_->ResetStartTime();
00216 
00217   // reset values
00218   IsComputed_ = false;
00219   Condest_ = -1.0;
00220 
00221   if (PolyDegree_ <= 0)
00222     IFPACK_CHK_ERR(-2); // at least one application
00223   
00224   if (IsRowMatrix_ && InvDiagonal_ == 0)
00225   {
00226     InvDiagonal_ = new Epetra_Vector(Matrix().Map());
00227 
00228     if (InvDiagonal_ == 0)
00229       IFPACK_CHK_ERR(-5);
00230 
00231     IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*InvDiagonal_));
00232 
00233     // Inverse diagonal elements
00234     // Replace zeros with 1.0
00235     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00236       double diag = (*InvDiagonal_)[i];
00237       if (IFPACK_ABS(diag) < MinDiagonalValue_)
00238         (*InvDiagonal_)[i] = MinDiagonalValue_;
00239       else
00240         (*InvDiagonal_)[i] = 1.0 / diag;
00241     }
00242   }
00243   // otherwise the inverse of the diagonal has been given by the user
00244 
00245   ComputeFlops_ += NumMyRows_;
00246 
00247   ++NumCompute_;
00248   ComputeTime_ += Time_->ElapsedTime();
00249   IsComputed_ = true;
00250 
00251   return(0);
00252 }
00253 
00254 //==============================================================================
00255 ostream& Ifpack_Chebyshev::Print(ostream & os) const
00256 {
00257 
00258   double MyMinVal, MyMaxVal;
00259   double MinVal, MaxVal;
00260 
00261   if (IsComputed_) {
00262     InvDiagonal_->MinValue(&MyMinVal);
00263     InvDiagonal_->MaxValue(&MyMaxVal);
00264     Comm().MinAll(&MyMinVal,&MinVal,1);
00265     Comm().MinAll(&MyMaxVal,&MaxVal,1);
00266   }
00267 
00268   if (!Comm().MyPID()) {
00269     os << endl;
00270     os << "================================================================================" << endl;
00271     os << "Ifpack_Chebyshev" << endl;
00272     os << "Degree of polynomial      = " << PolyDegree_ << endl;
00273     os << "Condition number estimate = " << Condest() << endl;
00274     os << "Global number of rows     = " << Operator_->OperatorRangeMap().NumGlobalElements() << endl;
00275     if (IsComputed_) {
00276       os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
00277       os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
00278     }
00279     if (ZeroStartingSolution_) 
00280       os << "Using zero starting solution" << endl;
00281     else
00282       os << "Using input starting solution" << endl;
00283     os << endl;
00284     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00285     os << "-----           -------   --------------       ------------     --------" << endl;
00286     os << "Initialize()    "   << std::setw(5) << NumInitialize_ 
00287        << "  " << std::setw(15) << InitializeTime_ 
00288        << "              0.0              0.0" << endl;
00289     os << "Compute()       "   << std::setw(5) << NumCompute_ 
00290        << "  " << std::setw(15) << ComputeTime_
00291        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_;
00292     if (ComputeTime_ != 0.0)
00293       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
00294     else
00295       os << "  " << std::setw(15) << 0.0 << endl;
00296     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse_ 
00297        << "  " << std::setw(15) << ApplyInverseTime_
00298        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
00299     if (ApplyInverseTime_ != 0.0)
00300       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
00301     else
00302       os << "  " << std::setw(15) << 0.0 << endl;
00303     os << "================================================================================" << endl;
00304     os << endl;
00305   }
00306 
00307   return(os);
00308 }
00309 
00310 //==============================================================================
00311 double Ifpack_Chebyshev::
00312 Condest(const Ifpack_CondestType CT, 
00313         const int MaxIters, const double Tol,
00314     Epetra_RowMatrix* Matrix)
00315 {
00316   if (!IsComputed()) // cannot compute right now
00317     return(-1.0);
00318 
00319   // always computes it. Call Condest() with no parameters to get
00320   // the previous estimate.
00321   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);
00322 
00323   return(Condest_);
00324 }
00325 
00326 //==============================================================================
00327 void Ifpack_Chebyshev::SetLabel()
00328 {
00329   Label_ = "IFPACK (Chebyshev polynomial), degree=" + Ifpack_toString(PolyDegree_);
00330 }
00331 
00332 //==============================================================================
00333 int Ifpack_Chebyshev::
00334 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00335 {
00336   if (!IsComputed())
00337     IFPACK_CHK_ERR(-3);
00338 
00339   if (PolyDegree_ == 0)
00340     return 0;
00341 
00342   int nVec = X.NumVectors();
00343   int len = X.MyLength();
00344   if (nVec != Y.NumVectors())
00345     IFPACK_CHK_ERR(-2);
00346 
00347   Time_->ResetStartTime();
00348 
00349   // AztecOO gives X and Y pointing to the same memory location,
00350   // need to create an auxiliary vector, Xcopy
00351   const Epetra_MultiVector* Xcopy;
00352   if (X.Pointers()[0] == Y.Pointers()[0])
00353     Xcopy = new Epetra_MultiVector(X);
00354   else
00355     Xcopy = &X;
00356 
00357   double **xPtr = 0, **yPtr = 0;
00358   Xcopy->ExtractView(&xPtr);
00359   Y.ExtractView(&yPtr);
00360 
00361   //--- Do a quick solve when the matrix is identity
00362   double *invDiag = InvDiagonal_->Values();
00363   if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
00364     if (nVec == 1) {
00365       double *yPointer = yPtr[0], *xPointer = xPtr[0];
00366       for (int i = 0; i < len; ++i)
00367         yPointer[i] = xPointer[i]*invDiag[i];
00368     }
00369     else {
00370       int i, k;
00371       for (i = 0; i < len; ++i) {
00372         double coeff = invDiag[i];
00373         for (k = 0; k < nVec; ++k)
00374           yPtr[k][i] = xPtr[k][i] * coeff;
00375       }
00376     } // if (nVec == 1)
00377     return 0;
00378   } // if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_))
00379 
00380   //--- Initialize coefficients
00381   // Note that delta stores the inverse of ML_Cheby::delta
00382   double alpha = LambdaMax_ / EigRatio_;
00383   double beta = 1.1 * LambdaMax_;
00384   double delta = 2.0 / (beta - alpha);
00385   double theta = 0.5 * (beta + alpha);
00386   double s1 = theta * delta;
00387 
00388   //--- Define vectors
00389   // In ML_Cheby, V corresponds to pAux and W to dk
00390   Epetra_MultiVector V(X);
00391   Epetra_MultiVector W(X);
00392 
00393   double *vPointer = V.Values(), *wPointer = W.Values();
00394 
00395   double oneOverTheta = 1.0/theta;
00396   int i, j, k;
00397 
00398   // Do the smoothing when block scaling is turned OFF
00399   // --- Treat the initial guess
00400   if (ZeroStartingSolution_ == false) {
00401     Operator_->Apply(Y, V);
00402     // Compute W = invDiag * ( X - V )/ Theta
00403     if (nVec == 1) {
00404       double *xPointer = xPtr[0];
00405       for (i = 0; i < len; ++i)
00406         wPointer[i] = invDiag[i] * (xPointer[i] - vPointer[i]) * oneOverTheta;
00407     }
00408     else {
00409       for (i = 0; i < len; ++i) {
00410         double coeff = invDiag[i]*oneOverTheta;
00411         double *wi = wPointer + i, *vi = vPointer + i;
00412         for (k = 0; k < nVec; ++k) {
00413           *wi = (xPtr[k][i] - (*vi)) * coeff;
00414           wi = wi + len; vi = vi + len;
00415         }
00416       }
00417     } // if (nVec == 1)
00418     // Update the vector Y
00419     Y.Update(1.0, W, 1.0);
00420   }
00421   else {
00422     // Compute W = invDiag * X / Theta
00423     if (nVec == 1) {
00424       double *xPointer = xPtr[0];
00425       for (i = 0; i < len; ++i)
00426         wPointer[i] = invDiag[i] * xPointer[i] * oneOverTheta;
00427       memcpy(yPtr[0], wPointer, len*sizeof(double));
00428     }
00429     else {
00430       for (i = 0; i < len; ++i) {
00431         double coeff = invDiag[i]*oneOverTheta;
00432         double *wi = wPointer + i;
00433         for (k = 0; k < nVec; ++k) {
00434           *wi = xPtr[k][i] * coeff;
00435           wi = wi + len;
00436         }
00437       }
00438       for (k = 0; k < nVec; ++k)
00439         memcpy(yPtr[k], wPointer + k*len, len*sizeof(double));
00440     } // if (nVec == 1)
00441   } // if (ZeroStartingSolution_ == false)
00442 
00443   //--- Apply the polynomial
00444   double rhok = 1.0/s1, rhokp1;
00445   double dtemp1, dtemp2;
00446   int degreeMinusOne = PolyDegree_ - 1;
00447   if (nVec == 1) {
00448     double *xPointer = xPtr[0];
00449     for (k = 0; k < degreeMinusOne; ++k) {
00450       Operator_->Apply(Y, V);
00451       rhokp1 = 1.0 / (2.0*s1 - rhok);
00452       dtemp1 = rhokp1 * rhok;
00453       dtemp2 = 2.0 * rhokp1 * delta;
00454       rhok = rhokp1;
00455       // Compute W = dtemp1 * W
00456       W.Scale(dtemp1);
00457       // Compute W = W + dtemp2 * invDiag * ( X - V )
00458       for (i = 0; i < len; ++i)
00459         wPointer[i] += dtemp2* invDiag[i] * (xPointer[i] - vPointer[i]);
00460       // Update the vector Y
00461       Y.Update(1.0, W, 1.0);
00462     } // for (k = 0; k < degreeMinusOne; ++k)
00463   }
00464   else {
00465     for (k = 0; k < degreeMinusOne; ++k) {
00466       Operator_->Apply(Y, V);
00467       rhokp1 = 1.0 / (2.0*s1 - rhok);
00468       dtemp1 = rhokp1 * rhok;
00469       dtemp2 = 2.0 * rhokp1 * delta;
00470       rhok = rhokp1;
00471       // Compute W = dtemp1 * W
00472       W.Scale(dtemp1);
00473       // Compute W = W + dtemp2 * invDiag * ( X - V )
00474       for (i = 0; i < len; ++i) {
00475         double coeff = invDiag[i]*dtemp2;
00476         double *wi = wPointer + i, *vi = vPointer + i;
00477         for (j = 0; j < nVec; ++j) {
00478           *wi += (xPtr[j][i] - (*vi)) * coeff;
00479           wi = wi + len; vi = vi + len;
00480         }
00481       }
00482       // Update the vector Y
00483       Y.Update(1.0, W, 1.0);
00484     } // for (k = 0; k < degreeMinusOne; ++k)
00485   } // if (nVec == 1)
00486 
00487   // Flops are updated in each of the following. 
00488 
00489   if (Xcopy != &X)
00490     delete Xcopy;
00491 
00492   ++NumApplyInverse_;
00493   ApplyInverseTime_ += Time_->ElapsedTime();
00494   return(0);
00495 }
00496 
00497 //==============================================================================
00498 int Ifpack_Chebyshev::
00499 PowerMethod(const Epetra_Operator& Operator, 
00500             const Epetra_Vector& InvPointDiagonal, 
00501             const int MaximumIterations, 
00502             double& lambda_max)
00503 {
00504   // this is a simple power method
00505   lambda_max = 0.0;
00506   double RQ_top, RQ_bottom, norm;
00507   Epetra_Vector x(Operator.OperatorDomainMap());
00508   Epetra_Vector y(Operator.OperatorRangeMap());
00509   x.Random();
00510   x.Norm2(&norm);
00511   if (norm == 0.0) IFPACK_CHK_ERR(-1);
00512 
00513   x.Scale(1.0 / norm);
00514 
00515   for (int iter = 0; iter < MaximumIterations; ++iter)
00516   {
00517     Operator.Apply(x, y);
00518     IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
00519     IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
00520     IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
00521     lambda_max = RQ_top / RQ_bottom;
00522     IFPACK_CHK_ERR(y.Norm2(&norm));
00523     if (norm == 0.0) IFPACK_CHK_ERR(-1);
00524     IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
00525   }
00526 
00527   return(0);
00528 }
00529 
00530 //==============================================================================
00531 int Ifpack_Chebyshev::
00532 CG(const Epetra_Operator& Operator, 
00533    const Epetra_Vector& InvPointDiagonal, 
00534    const int MaximumIterations, 
00535    double& lambda_min, double& lambda_max)
00536 {
00537 #ifdef HAVE_IFPACK_AZTECOO
00538   Epetra_Vector x(Operator.OperatorDomainMap());
00539   Epetra_Vector y(Operator.OperatorRangeMap());
00540   x.Random();
00541   y.PutScalar(0.0);
00542 
00543   Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
00544   AztecOO solver(LP);
00545   solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
00546   solver.SetAztecOption(AZ_output, AZ_none);
00547 
00548   Ifpack_DiagPreconditioner diag(Operator.OperatorDomainMap(),
00549                                  Operator.OperatorRangeMap(),
00550                                  InvPointDiagonal);
00551   solver.SetPrecOperator(&diag);
00552   solver.Iterate(MaximumIterations, 1e-10);
00553 
00554   const double* status = solver.GetAztecStatus();
00555 
00556   lambda_min = status[AZ_lambda_min];
00557   lambda_max = status[AZ_lambda_max];
00558 
00559   return(0);
00560 #else
00561   cout << "You need to configure IFPACK with support for AztecOO" << endl;
00562   cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
00563   cout << "in your configure script." << endl;
00564   IFPACK_CHK_ERR(-1);
00565 #endif
00566 }
00567 
00568 #endif

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1