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

Generated on Tue Oct 20 12:48:53 2009 for IFPACK by doxygen 1.4.7