Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_IHSS.cpp
Go to the documentation of this file.
00001 #include "Ifpack_IHSS.h"
00002 #include "Ifpack.h"
00003 #include "Ifpack_Utils.h"
00004 #include "Teuchos_ParameterList.hpp"
00005 #include "Teuchos_RefCountPtr.hpp"
00006 
00007 
00008 
00009 using Teuchos::RefCountPtr;
00010 using Teuchos::rcp;
00011 
00012 
00013 #ifdef HAVE_IFPACK_EPETRAEXT
00014 #include "EpetraExt_MatrixMatrix.h"
00015 
00016 
00017 Ifpack_IHSS::Ifpack_IHSS(Epetra_RowMatrix* A):
00018   IsInitialized_(false),
00019   IsComputed_(false),
00020   Label_(),
00021   EigMaxIters_(10),
00022   EigRatio_(30.0),
00023   LambdaMax_(-1.0),
00024   Alpha_(-1.0),
00025   NumSweeps_(1),
00026 
00027   Time_(A->Comm())
00028 {
00029   Epetra_CrsMatrix *Acrs=dynamic_cast<Epetra_CrsMatrix*>(A);
00030   TEST_FOR_EXCEPT(!Acrs) 
00031   A_=rcp(Acrs,false);
00032 }
00033 
00034 void Ifpack_IHSS::Destroy(){
00035 }
00036 
00037 
00038 
00039 int Ifpack_IHSS::Initialize(){
00040   EigMaxIters_          = List_.get("ihss: eigenvalue max iterations",EigMaxIters_);
00041   EigRatio_             = List_.get("ihss: ratio eigenvalue", EigRatio_);
00042   NumSweeps_            = List_.get("ihss: sweeps",NumSweeps_);
00043 
00044   // Counters
00045   IsInitialized_=true;
00046   NumInitialize_++;
00047   return 0;
00048 }
00049 
00050 int Ifpack_IHSS::SetParameters(Teuchos::ParameterList& parameterlist){
00051   List_=parameterlist;
00052   return 0;
00053 }
00054 
00055 
00056 int Ifpack_IHSS::Compute(){
00057   if(!IsInitialized_) Initialize();
00058 
00059   int rv;
00060   Ifpack Factory;
00061   Epetra_CrsMatrix *Askew=0,*Aherm=0;
00062   Ifpack_Preconditioner *Pskew=0, *Pherm=0;
00063   Time_.ResetStartTime();
00064 
00065   // Create Aherm (w/o diagonal)
00066   rv=EpetraExt::MatrixMatrix::Add(*A_,false,.5,*A_,true,.5,Aherm);
00067   Aherm->FillComplete();
00068   if(rv) IFPACK_CHK_ERR(-1); 
00069 
00070   // Grab Aherm's diagonal
00071   Epetra_Vector avec(Aherm->RowMap());
00072   IFPACK_CHK_ERR(Aherm->ExtractDiagonalCopy(avec));
00073 
00074 
00075   // Compute alpha using the Bai, Golub & Ng 2003 formula, not the more multigrid-appropriate Hamilton, Benzi and Haber 2007.
00076   //  PowerMethod(Aherm, EigMaxIters_,LambdaMax_);
00077   //  Alpha_=LambdaMax_ / sqrt(EigRatio_);
00078 
00079   // Try something more Hamilton inspired, using the maximum diagonal value of Aherm.
00080   avec.MaxValue(&Alpha_);
00081 
00082   // Add alpha to the diagonal of Aherm
00083   for(int i=0;i<Aherm->NumMyRows();i++) avec[i]+=Alpha_;   
00084   IFPACK_CHK_ERR(Aherm->ReplaceDiagonalValues(avec));
00085   Aherm_=rcp(Aherm);
00086 
00087   // Compute Askew (and add diagonal)
00088   Askew=new Epetra_CrsMatrix(Copy,A_->RowMap(),0);
00089   rv=EpetraExt::MatrixMatrix::Add(*A_,false,.5,*A_,true,-.5,Askew);
00090   if(rv) IFPACK_CHK_ERR(-2); 
00091   for(int i=0;i<Askew->NumMyRows();i++) {
00092     int gid=Askew->GRID(i);
00093     Askew->InsertGlobalValues(gid,1,&Alpha_,&gid);
00094   }
00095   Askew->FillComplete();
00096   Askew_=rcp(Askew);
00097 
00098   // Compute preconditioner for Aherm
00099   Teuchos::ParameterList PLh=List_.sublist("ihss: hermetian list");
00100   string htype=List_.get("ihss: hermetian type","ILU");
00101   Pherm= Factory.Create(htype, Aherm);
00102   Pherm->SetParameters(PLh);
00103   IFPACK_CHK_ERR(Pherm->Compute());
00104   Pherm_=rcp(Pherm);
00105 
00106   // Compute preconditoner for Askew 
00107   Teuchos::ParameterList PLs=List_.sublist("ihss: skew hermetian list");
00108   string stype=List_.get("ihss: skew hermetian type","ILU");
00109   Pskew= Factory.Create(stype, Askew);
00110   Pskew->SetParameters(PLs);
00111   IFPACK_CHK_ERR(Pskew->Compute());
00112   Pskew_=rcp(Pskew);
00113   
00114   // Label
00115   sprintf(Label_, "IFPACK IHSS (H,S)=(%s/%s)",htype.c_str(),stype.c_str()); 
00116 
00117   // Counters
00118   IsComputed_=true;
00119   NumCompute_++;
00120   ComputeTime_ += Time_.ElapsedTime();
00121   return 0;
00122 }
00123 
00124 
00125 int Ifpack_IHSS::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{  
00126   if(!IsComputed_) return -1;
00127   Time_.ResetStartTime();
00128   bool initial_guess_is_zero=false;  
00129 
00130   // AztecOO gives X and Y pointing to the same memory location,
00131   // need to create an auxiliary vector, Xcopy
00132   Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
00133   Epetra_MultiVector Temp(X);
00134   if (X.Pointers()[0] == Y.Pointers()[0]){
00135     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00136     // Since the user didn't give us anything better, our initial guess is zero.
00137     Y.Scale(0.0);
00138     initial_guess_is_zero=true;
00139   }
00140   else
00141     Xcopy = Teuchos::rcp( &X, false );
00142 
00143   Epetra_MultiVector T1(Y),T2(*Xcopy);
00144 
00145   // Note: Since Aherm and Askew are actually (aI+H) and (aI+S) accordingly (to save memory), 
00146   // the application thereof needs to be a little different.
00147   // The published algorithm is:
00148   // temp = (aI+H)^{-1} [ (aI-S) y + x ]
00149   // y = (aI+S)^{-1} [ (aI-H) temp + x ]
00150   //
00151   // But we're doing:
00152   // temp = (aI+H)^{-1} [ 2a y - Shat y + x ]
00153   // y = (aI+S)^{-1} [ 2 a temp - Hhat temp + x ]
00154 
00155   for(int i=0;i<NumSweeps_;i++){
00156     // temp = (aI+H)^{-1} [ 2a y - Shat y + x ]
00157     if(!initial_guess_is_zero || i >0 ){
00158       Askew_->Apply(Y,T1);
00159       T2.Update(2*Alpha_,Y,-1,T1,1);
00160     }
00161     Pherm_->ApplyInverse(T2,Y);
00162     
00163     // y = (aI+S)^{-1} [ 2 a temp - Hhat temp + x ]
00164     Aherm_->Apply(Y,T1);
00165     T2.Scale(1.0,*Xcopy);
00166     T2.Update(2*Alpha_,Y,-1,T1,1.0);
00167     Pskew_->ApplyInverse(T2,Y);  
00168   }
00169 
00170   // Counter update
00171   NumApplyInverse_++;
00172   ApplyInverseTime_ += Time_.ElapsedTime();
00173   return 0;
00174 }
00175 
00176 
00177 ostream& Ifpack_IHSS::Print(ostream& os) const{
00178   os<<"Ifpack_IHSS"<<endl;
00179   os<<"-Hermetian preconditioner"<<endl;
00180   os<<"-Skew Hermetian preconditioner"<<endl;
00181   os<<endl;
00182   return os;
00183 }
00184 
00185 
00186 double Ifpack_IHSS::Condest(const Ifpack_CondestType CT, 
00187                              const int MaxIters,
00188                              const double Tol,
00189                              Epetra_RowMatrix* Matrix_in){
00190   return -1.0;
00191 }
00192 
00193 
00194 int Ifpack_IHSS::PowerMethod(Epetra_Operator * Op,const int MaximumIterations,  double& lambda_max)
00195 {
00196 
00197   // this is a simple power method
00198   lambda_max = 0.0;
00199   double RQ_top, RQ_bottom, norm;
00200   Epetra_Vector x(Op->OperatorDomainMap());
00201   Epetra_Vector y(Op->OperatorRangeMap());
00202   x.Random();
00203   x.Norm2(&norm);
00204   if (norm == 0.0) IFPACK_CHK_ERR(-1);
00205 
00206   x.Scale(1.0 / norm);
00207 
00208   for (int iter = 0; iter < MaximumIterations; ++iter)
00209   {
00210     Op->Apply(x, y);
00211     IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
00212     IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
00213     lambda_max = RQ_top / RQ_bottom;
00214     IFPACK_CHK_ERR(y.Norm2(&norm));
00215     if (norm == 0.0) IFPACK_CHK_ERR(-1);
00216     IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
00217   }
00218 
00219   return(0);
00220 }
00221 
00222 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines