Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_SORa.cpp
Go to the documentation of this file.
00001 #include "Ifpack_SORa.h"
00002 #include "Ifpack.h"
00003 #include "Ifpack_Utils.h"
00004 #include "Teuchos_ParameterList.hpp"
00005 #include "Teuchos_RefCountPtr.hpp"
00006 #include "Epetra_Import.h"
00007 #include "Epetra_Export.h"
00008 #include "Epetra_IntVector.h"
00009 
00010 using Teuchos::RefCountPtr;
00011 using Teuchos::rcp;
00012 
00013 
00014 
00015 #ifdef HAVE_IFPACK_EPETRAEXT
00016 #include "EpetraExt_MatrixMatrix.h"
00017 #include "EpetraExt_RowMatrixOut.h"
00018 #include "EpetraExt_MultiVectorOut.h"
00019 #include "EpetraExt_Transpose_RowMatrix.h"
00020 
00021 
00022 #define ABS(x) ((x)>=0?(x):-(x))
00023 #define MIN(x,y) ((x)<(y)?(x):(y))
00024 #define MAX(x,y) ((x)>(y)?(x):(y))
00025 
00026 // Useful functions horked from ML
00027 int* FindLocalDiricheltRowsFromOnesAndZeros(const Epetra_CrsMatrix & Matrix, int &numBCRows);
00028 void Apply_BCsToMatrixRowsAndColumns(const int *dirichletRows, int numBCRows,const Epetra_IntVector &dirichletColumns,const Epetra_CrsMatrix & Matrix);
00029 Epetra_IntVector * FindLocalDirichletColumnsFromRows(const int *dirichletRows, int numBCRows,const Epetra_CrsMatrix & Matrix);
00030 
00031 Ifpack_SORa::Ifpack_SORa(Epetra_RowMatrix* A):
00032   IsInitialized_(false),
00033   IsComputed_(false),
00034   Label_(),
00035   Alpha_(1.5),
00036   Gamma_(1.0),
00037   NumSweeps_(1),
00038   IsParallel_(false),
00039   HaveOAZBoundaries_(false),
00040   UseInterprocDamping_(false),
00041   UseGlobalDamping_(false),
00042   LambdaMax_(-1.0),
00043   Time_(A->Comm())
00044 {
00045   Epetra_CrsMatrix *Acrs=dynamic_cast<Epetra_CrsMatrix*>(A);
00046   if(Acrs) Acrs_=rcp(Acrs,false);
00047   A_=rcp(A,false);
00048 }
00049 
00050 void Ifpack_SORa::Destroy(){
00051 }
00052 
00053 
00054 
00055 int Ifpack_SORa::Initialize(){
00056   Alpha_            = List_.get("sora: alpha", Alpha_);
00057   Gamma_            = List_.get("sora: gamma",Gamma_);
00058   NumSweeps_        = List_.get("sora: sweeps",NumSweeps_);
00059   HaveOAZBoundaries_= List_.get("sora: oaz boundaries", HaveOAZBoundaries_);
00060   UseInterprocDamping_ = List_.get("sora: use interproc damping",UseInterprocDamping_);
00061   UseGlobalDamping_ = List_.get("sora: use global damping",UseGlobalDamping_);
00062 
00063   if (A_->Comm().NumProc() != 1) IsParallel_ = true;
00064   else {
00065     IsParallel_ = false;    
00066     // Don't use interproc damping, for obvious reasons
00067     UseInterprocDamping_=false;
00068   }
00069 
00070   // Counters
00071   IsInitialized_=true;
00072   NumInitialize_++;
00073   return 0;
00074 }
00075 
00076 int Ifpack_SORa::SetParameters(Teuchos::ParameterList& parameterlist){
00077   List_=parameterlist;
00078   return 0;
00079 }
00080 
00081 
00082 int Ifpack_SORa::Compute(){
00083   if(!IsInitialized_) Initialize();
00084   Epetra_Map *RowMap=const_cast<Epetra_Map*>(&A_->RowMatrixRowMap());
00085   Epetra_Vector Adiag(*RowMap);
00086   Epetra_CrsMatrix *Askew2=0, *Aherm2=0,*W=0;
00087   int *rowptr_s,*colind_s,*rowptr_h,*colind_h;
00088   double *vals_s,*vals_h;
00089   bool RowMatrixMode=(Acrs_==Teuchos::null);
00090 
00091   // Label
00092   sprintf(Label_, "IFPACK SORa (alpha=%5.2f gamma=%5.2f)",Alpha_,Gamma_); 
00093 
00094   if(RowMatrixMode){
00095     if(!A_->Comm().MyPID()) cout<<"SORa: RowMatrix mode enabled"<<endl;
00096     // RowMatrix mode, build Acrs_ the hard way.
00097     Epetra_RowMatrix *Arow=&*A_;
00098     Epetra_Map *ColMap=const_cast<Epetra_Map*>(&A_->RowMatrixColMap());
00099 
00100     int Nmax=A_->MaxNumEntries();
00101     int length;
00102     vector<int> indices(Nmax);
00103     vector<double> values(Nmax); 
00104     Epetra_CrsMatrix *Acrs=new Epetra_CrsMatrix(Copy,*RowMap,Nmax);
00105 
00106     for(int i=0;i<Arow->NumMyRows();i++){
00107       Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices[0]);
00108       for(int j=0;j<length;j++)
00109   indices[j]=ColMap->GID(indices[j]);
00110       Acrs->InsertGlobalValues(RowMap->GID(i),length,&values[0],&indices[0]);
00111     }
00112     Acrs->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
00113     Acrs_=rcp(Acrs,true);
00114   }
00115 
00116   // Create Askew2
00117   // Note: Here I let EpetraExt do the thinking for me.  Since it gets all the maps correct for the E + F^T stencil.
00118   // There are probably more efficient ways to do this but this method has the bonus of allowing code reuse.
00119   IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,-1,Askew2));
00120   Askew2->FillComplete();
00121   
00122   // Create Aherm2
00123   IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,1,Aherm2));
00124   Aherm2->FillComplete();
00125 
00126   int nnz2=Askew2->NumMyNonzeros();
00127   int N=Askew2->NumMyRows();
00128 
00129 
00130   // Grab pointers
00131   IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
00132   IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));
00133 
00134   // Sanity checking: Make sure the sparsity patterns are the same
00135 #define SANITY_CHECK
00136 #ifdef SANITY_CHECK
00137   for(int i=0;i<N;i++)
00138     if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
00139   for(int i=0;i<nnz2;i++)
00140     if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
00141 #endif
00142 
00143   // Dirichlet Detection & Nuking of Aherm2 and Askew2
00144   // Note: Relies on Aherm2/Askew2 having identical sparsity patterns (see sanity check above)
00145   if(HaveOAZBoundaries_){ 
00146     int numBCRows;
00147     int* dirRows=FindLocalDiricheltRowsFromOnesAndZeros(*Acrs_,numBCRows);
00148     Epetra_IntVector* dirCols=FindLocalDirichletColumnsFromRows(dirRows,numBCRows,*Aherm2);
00149     Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Aherm2);
00150     Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Askew2);
00151     delete [] dirRows;
00152     delete dirCols;
00153   }
00154 
00155   // Grab diagonal of A
00156   A_->ExtractDiagonalCopy(Adiag);
00157 
00158   // Allocate the diagonal for W
00159   Epetra_Vector *Wdiag = new Epetra_Vector(*RowMap);
00160 
00161   // Build the W matrix (lower triangle only)
00162   // Note: Relies on EpetraExt giving me identical sparsity patterns for both Askew2 and Aherm2 (see sanity check above)
00163   int maxentries=Askew2->MaxNumEntries();
00164   int* gids=new int [maxentries];
00165   double* newvals=new double[maxentries];
00166   W=new Epetra_CrsMatrix(Copy,*RowMap,0);
00167   for(int i=0;i<N;i++){
00168     // Build the - (1+alpha)/2 E - (1-alpha)/2 F part of the W matrix
00169     int rowgid=Acrs_->GRID(i);
00170     double c_data=0.0;
00171     double ipdamp=0.0;
00172     int idx=0;
00173 
00174     for(int j=rowptr_s[i];j<rowptr_s[i+1];j++){      
00175       int colgid=Askew2->GCID(colind_s[j]);
00176       c_data+=fabs(vals_s[j]);
00177       if(rowgid>colgid){
00178   // Rely on the fact that off-diagonal entries are always numbered last, dropping the entry entirely.
00179   if(colind_s[j] < N) {       
00180     gids[idx]=colgid;
00181     newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
00182     idx++;
00183   }
00184   else{
00185     ipdamp+=fabs(vals_h[j]/2 + Alpha_ * vals_s[j]/2);
00186   }
00187       }   
00188     }
00189     if(idx>0)
00190       IFPACK_CHK_ERR(W->InsertGlobalValues(rowgid,idx,newvals,gids));
00191 
00192     // Do the diagonal
00193     double w_val= c_data*Alpha_*Gamma_/4 + Adiag[Acrs_->LRID(rowgid)];
00194     if(UseInterprocDamping_) w_val+=ipdamp;
00195 
00196     W->InsertGlobalValues(rowgid,1,&w_val,&rowgid);
00197     IFPACK_CHK_ERR(Wdiag->ReplaceGlobalValues(1,&w_val,&rowgid));
00198   }
00199   W->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());      
00200   W_=rcp(W);
00201   Wdiag_=rcp(Wdiag);
00202 
00203   // Mark as computed
00204   IsComputed_=true;
00205 
00206   // Global damping, if wanted
00207   if(UseGlobalDamping_) {
00208     PowerMethod(10,LambdaMax_);
00209     if(!A_->Comm().MyPID()) printf("SORa: Global damping parameter = %6.4e (lmax=%6.4e)\n",GetOmega(),LambdaMax_);
00210   }
00211 
00212   // Cleanup
00213   delete [] gids;
00214   delete [] newvals;
00215   delete Aherm2;
00216   delete Askew2;
00217   if(RowMatrixMode) {
00218     Acrs_=Teuchos::null;
00219     transposer2_=Teuchos::null;
00220   }
00221 
00222   // Counters
00223   NumCompute_++;
00224   ComputeTime_ += Time_.ElapsedTime();
00225   return 0;
00226 }
00227 
00228 
00229 
00230 int Ifpack_SORa::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{  
00231   if(!IsComputed_) return -1;
00232   Time_.ResetStartTime();
00233   bool initial_guess_is_zero=false;  
00234   int NumMyRows=W_->NumMyRows();
00235   int NumVectors = X.NumVectors();
00236   Epetra_MultiVector Temp(A_->RowMatrixRowMap(),NumVectors);
00237 
00238   double omega=GetOmega();
00239 
00240   // need to create an auxiliary vector, Xcopy
00241   Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
00242   if (X.Pointers()[0] == Y.Pointers()[0]){
00243     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00244     // Since the user didn't give us anything better, our initial guess is zero.
00245     Y.Scale(0.0);
00246     initial_guess_is_zero=true;
00247   }
00248   else
00249     Xcopy = Teuchos::rcp( &X, false ); 
00250 
00251   Teuchos::RefCountPtr< Epetra_MultiVector > T2;
00252   // Note: Assuming that the matrix has an importer.  Ifpack_PointRelaxation doesn't do this, but given that 
00253   // I have a CrsMatrix, I'm probably OK.
00254   // Note: This is the lazy man's version sacrificing a few extra flops for avoiding if statements to determine 
00255   // if things are on or off processor.
00256   // Note: T2 must be zero'd out
00257   if (IsParallel_ && W_->Importer())  T2 = Teuchos::rcp( new Epetra_MultiVector(W_->Importer()->TargetMap(),NumVectors,true));
00258   else T2 = Teuchos::rcp( new Epetra_MultiVector(A_->RowMatrixRowMap(),NumVectors,true));
00259 
00260   // Pointer grabs
00261   int* rowptr,*colind;
00262   double *values;
00263   double **t_ptr,** y_ptr, ** t2_ptr, **x_ptr,*d_ptr;
00264   T2->ExtractView(&t2_ptr);
00265   Y.ExtractView(&y_ptr);
00266   Temp.ExtractView(&t_ptr);
00267   Xcopy->ExtractView(&x_ptr);
00268   Wdiag_->ExtractView(&d_ptr);
00269   IFPACK_CHK_ERR(W_->ExtractCrsDataPointers(rowptr,colind,values));
00270 
00271 
00272   for(int i=0; i<NumSweeps_; i++){
00273     // Calculate b-Ax 
00274     if(!initial_guess_is_zero  || i > 0) {      
00275       A_->Apply(Y,Temp);
00276       Temp.Update(1.0,*Xcopy,-1.0);
00277     }
00278     else 
00279       Temp.Update(1.0,*Xcopy,0.0);
00280 
00281     // Note: The off-processor entries of T2 never get touched (they're always zero) and the other entries are updated 
00282     // in this sweep before they are used, so we don't need to reset T2 to zero here.
00283 
00284     // Do backsolve & update
00285     // x = x  + W^{-1} (b - A x)
00286     for(int j=0; j<NumMyRows; j++){
00287       double diag=d_ptr[j];
00288       for (int m=0 ; m<NumVectors; m++) {
00289   double dtmp=0.0;
00290   // Note: Since the diagonal is in the matrix, we need to zero that entry of T2 here to make sure it doesn't contribute.
00291   t2_ptr[m][j]=0.0;
00292   for(int k=rowptr[j];k<rowptr[j+1];k++){
00293     dtmp+= values[k]*t2_ptr[m][colind[k]];
00294   }
00295   // Yes, we need to update both of these.
00296   t2_ptr[m][j] = (t_ptr[m][j]- dtmp)/diag;     
00297   y_ptr[m][j] += omega*t2_ptr[m][j];
00298       }
00299     }
00300   }
00301 
00302   // Counter update
00303   NumApplyInverse_++;
00304   ApplyInverseTime_ += Time_.ElapsedTime();
00305   return 0;
00306 }
00307 
00308 
00309 ostream& Ifpack_SORa::Print(ostream& os) const{
00310   os<<"Ifpack_SORa"<<endl;
00311   os<<" alpha = "<<Alpha_<<endl;
00312   os<<" gamma = "<<Gamma_<<endl;
00313   os<<endl;
00314   return os;
00315 }
00316 
00317 
00318 double Ifpack_SORa::Condest(const Ifpack_CondestType CT, 
00319                              const int MaxIters,
00320                              const double Tol,
00321                              Epetra_RowMatrix* Matrix_in){
00322   return -1.0;
00323 }
00324 
00325 
00326 
00327 
00328 
00329 // ============================================================================
00330 inline int* FindLocalDiricheltRowsFromOnesAndZeros(const Epetra_CrsMatrix & Matrix, int &numBCRows){
00331   int *dirichletRows = new int[Matrix.NumMyRows()];
00332   numBCRows = 0;
00333   for (int i=0; i<Matrix.NumMyRows(); i++) {
00334     int numEntries, *cols;
00335     double *vals;
00336     int ierr = Matrix.ExtractMyRowView(i,numEntries,vals,cols);
00337     if (ierr == 0) {
00338       int nz=0;
00339       for (int j=0; j<numEntries; j++) if (vals[j]!=0) nz++;      
00340       if (nz==1) dirichletRows[numBCRows++] = i;           
00341       // EXPERIMENTAL: Treat Special Inflow Boundaries as Dirichlet Boundaries
00342       if(nz==2) dirichletRows[numBCRows++] = i;        
00343     }/*end if*/
00344   }/*end for*/
00345   return dirichletRows;
00346 }/*end FindLocalDiricheltRowsFromOnesAndZeros*/
00347 
00348 
00349 // ====================================================================== 
00351 inline Epetra_IntVector * FindLocalDirichletColumnsFromRows(const int *dirichletRows, int numBCRows,const Epetra_CrsMatrix & Matrix){
00352   // Zero the columns corresponding to the Dirichlet rows.  Completely ignore the matrix maps.
00353   
00354   // Build rows
00355   Epetra_IntVector ZeroRows(Matrix.RowMap());
00356   Epetra_IntVector *ZeroCols=new Epetra_IntVector(Matrix.ColMap());
00357   ZeroRows.PutValue(0);  
00358   ZeroCols->PutValue(0);  
00359 
00360   // Easy Case: We're all on one processor
00361   if(Matrix.RowMap().SameAs(Matrix.ColMap())){
00362     for (int i=0; i < numBCRows; i++)
00363       (*ZeroCols)[dirichletRows[i]]=1;
00364     return ZeroCols;
00365   }
00366 
00367   // Flag the rows which are zero locally
00368   for (int i=0; i < numBCRows; i++)
00369     ZeroRows[dirichletRows[i]]=1;
00370 
00371   // Boundary exchange to move the data  
00372   if(Matrix.RowMap().SameAs(Matrix.DomainMap())){
00373     // Use A's Importer if we have one
00374     ZeroCols->Import(ZeroRows,*Matrix.Importer(),Insert);     
00375   }
00376   else{
00377     // Use custom importer if we don't
00378     Epetra_Import Importer(Matrix.ColMap(),Matrix.RowMap());
00379     ZeroCols->Import(ZeroRows,Importer,Insert);      
00380   }
00381   return ZeroCols;
00382 }
00383 
00384 
00385 // ====================================================================== 
00386 inline void Apply_BCsToMatrixRowsAndColumns(const int *dirichletRows, int numBCRows,const Epetra_IntVector &dirichletColumns,const Epetra_CrsMatrix & Matrix){
00387   /* This function zeros out rows & columns of Matrix.
00388      Comments: The graph of Matrix is unchanged.
00389   */
00390   // Nuke the rows
00391   for(int i=0;i<numBCRows;i++){
00392     int numEntries, *cols;
00393     double *vals;
00394     Matrix.ExtractMyRowView(dirichletRows[i],numEntries,vals,cols);
00395     for (int j=0; j<numEntries; j++) vals[j]=0.0; 
00396   }/*end for*/
00397 
00398   // Nuke the columns
00399   for (int i=0; i < Matrix.NumMyRows(); i++) {
00400     int numEntries;
00401     double *vals;
00402     int *cols;
00403     Matrix.ExtractMyRowView(i,numEntries,vals,cols);
00404     for (int j=0; j < numEntries; j++) {
00405       if (dirichletColumns[ cols[j] ] > 0)  vals[j] = 0.0;
00406     }/*end for*/
00407   }/*end for*/
00408 }/* end Apply_BCsToMatrixColumns */
00409 
00410 
00411 
00412 
00413 
00414 int Ifpack_SORa::
00415 PowerMethod(const int MaximumIterations,  double& lambda_max)
00416 {
00417   // this is a simple power method
00418   lambda_max = 0.0;
00419   double RQ_top, RQ_bottom, norm;
00420   Epetra_Vector x(A_->OperatorDomainMap());
00421   Epetra_Vector y(A_->OperatorRangeMap());
00422   Epetra_Vector z(A_->OperatorRangeMap());
00423   x.Random();
00424   x.Norm2(&norm);
00425   if (norm == 0.0) IFPACK_CHK_ERR(-1);
00426 
00427   x.Scale(1.0 / norm);
00428 
00429   // Only do 1 sweep per PM step, turn off global damping
00430   int NumSweepsBackup=NumSweeps_;
00431   bool UseGlobalDampingBackup=UseGlobalDamping_;
00432   NumSweeps_=1;UseGlobalDamping_=false;
00433 
00434   for (int iter = 0; iter < MaximumIterations; ++iter)
00435   {
00436     y.PutScalar(0.0);
00437     A_->Apply(x, z);
00438     ApplyInverse(z,y);
00439     y.Update(1.0,x,-1.0);
00440 
00441     // Compute the Rayleigh Quotient
00442     y.Dot(x, &RQ_top);
00443     x.Dot(x, &RQ_bottom);
00444     lambda_max = RQ_top / RQ_bottom;   
00445     y.Norm2(&norm);
00446     if (norm == 0.0) IFPACK_CHK_ERR(-1);
00447     x.Update(1.0 / norm, y, 0.0);
00448 
00449   }
00450 
00451   // Enable if we want to prevent over-relaxation
00452   //  lambda_max=MAX(lambda_max,1.0);
00453 
00454   // Reset parameters
00455   NumSweeps_=NumSweepsBackup;
00456   UseGlobalDamping_=UseGlobalDampingBackup;
00457 
00458   return(0);
00459 }
00460 
00461 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines