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