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 
00008 
00009 using Teuchos::RefCountPtr;
00010 using Teuchos::rcp;
00011 
00012 
00013 #ifdef HAVE_IFPACK_EPETRAEXT
00014 #include "EpetraExt_MatrixMatrix.h"
00015 #include "EpetraExt_RowMatrixOut.h"
00016 #include "EpetraExt_MultiVectorOut.h"
00017 
00018 
00019 #define ABS(x) ((x)>=0?(x):-(x))
00020 
00021 Ifpack_SORa::Ifpack_SORa(Epetra_RowMatrix* A):
00022   IsInitialized_(false),
00023   IsComputed_(false),
00024   Label_(),
00025   Alpha_(1.5),
00026   Gamma_(1.0),
00027   NumSweeps_(1),
00028   IsParallel_(false),
00029   Time_(A->Comm())
00030 {
00031   Epetra_CrsMatrix *Acrs=dynamic_cast<Epetra_CrsMatrix*>(A);
00032   TEST_FOR_EXCEPT(!Acrs) 
00033   A_=rcp(Acrs,false);
00034 }
00035 
00036 void Ifpack_SORa::Destroy(){
00037 }
00038 
00039 
00040 
00041 int Ifpack_SORa::Initialize(){
00042   Alpha_            = List_.get("sora: alpha", Alpha_);
00043   Gamma_            = List_.get("sora: gamma",Gamma_);
00044   NumSweeps_        = List_.get("sora: sweeps",NumSweeps_);
00045 
00046   if (A_->Comm().NumProc() != 1) IsParallel_ = true;
00047   else IsParallel_ = false;
00048 
00049   // Counters
00050   IsInitialized_=true;
00051   NumInitialize_++;
00052   return 0;
00053 }
00054 
00055 int Ifpack_SORa::SetParameters(Teuchos::ParameterList& parameterlist){
00056   List_=parameterlist;
00057   return 0;
00058 }
00059 
00060 
00061 int Ifpack_SORa::Compute(){
00062   if(!IsInitialized_) Initialize();
00063 
00064   
00065   Epetra_CrsMatrix *Askew2=0, *Aherm2=0,*W=0;
00066   int *rowptr_s,*colind_s,*rowptr_h,*colind_h;
00067   double *vals_s,*vals_h;
00068   Epetra_Vector Adiag(A_->RowMap());
00069 
00070   // Label
00071   sprintf(Label_, "IFPACK SORa (alpha=%5.2f gamma=%5.2f)",Alpha_,Gamma_); 
00072 
00073   // Create Askew2
00074   // Note: Here I let EpetraExt do the thinking for me.  Since it gets all the maps correct for the E + F^T stencil.
00075   // There are probably more efficient ways to do this but this method has the bonus of allowing code reuse.
00076   IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*A_,false,1,*A_,true,-1,Askew2));
00077   Askew2->FillComplete();
00078   int nnz2=Askew2->NumMyNonzeros();
00079   int N=Askew2->NumMyRows();
00080 
00081   // Create Aherm2
00082   IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*A_,false,1,*A_,true,1,Aherm2));
00083   Aherm2->FillComplete();
00084 
00085   // Grab pointers
00086   IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
00087   IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));
00088 
00089   // Sanity checking: Make sure the sparsity patterns are the same
00090 #define SANITY_CHECK
00091 #ifdef SANITY_CHECK
00092   for(int i=0;i<N;i++)
00093     if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
00094   for(int i=0;i<nnz2;i++)
00095     if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
00096 #endif
00097 
00098   // Grab diagonal of A
00099   A_->ExtractDiagonalCopy(Adiag);
00100 
00101   // Allocate the diagonal for W
00102   Epetra_Vector *Wdiag = new Epetra_Vector(A_->RowMap());
00103 
00104   // Build the W matrix (strict lower triangle only)
00105   // Note: Relies on EpetraExt giving me identical sparsity patterns for both Askew2 and Aherm2 (see sanity check above)
00106   int maxentries=Askew2->MaxNumEntries();
00107   int* gids=new int [maxentries];
00108   double* newvals=new double[maxentries];
00109   W=new Epetra_CrsMatrix(Copy,A_->RowMap(),0);
00110   for(int i=0;i<N;i++){
00111     // Build the - (1+alpha)/2 E - (1-alpha)/2 F part of the W matrix
00112     int rowgid=A_->GRID(i);
00113     double c_data=0.0;
00114     int idx=0;
00115     for(int j=rowptr_s[i];j<rowptr_s[i+1];j++){      
00116       int colgid=Askew2->GCID(colind_s[j]);
00117       c_data+=fabs(vals_s[j]);
00118       if(rowgid>colgid){
00119   gids[idx]=colgid;
00120   newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
00121   idx++;
00122       }
00123     }
00124     IFPACK_CHK_ERR(W->InsertGlobalValues(rowgid,idx,newvals,gids));
00125     // Do the diagonal
00126     double w_val= c_data*Alpha_*Gamma_/4 + Adiag[A_->LRID(rowgid)];
00127 
00128     // Note: I drop a zero on the diagonal just to make sure that my rowmap is a subset of my column map.
00129     double zero=0.0;
00130     W->InsertGlobalValues(rowgid,1,&zero,&rowgid);//HAQ
00131     IFPACK_CHK_ERR(Wdiag->ReplaceGlobalValues(1,&w_val,&rowgid));
00132   }
00133   W->FillComplete(A_->DomainMap(),A_->RangeMap());      
00134   W_=rcp(W);
00135   Wdiag_=rcp(Wdiag);
00136 
00137   // Cleanup
00138   delete [] gids;
00139   delete [] newvals;
00140   delete Aherm2;
00141   delete Askew2;
00142 
00143   // Counters
00144   IsComputed_=true;
00145   NumCompute_++;
00146   ComputeTime_ += Time_.ElapsedTime();
00147   return 0;
00148 }
00149 
00150 
00151 int Ifpack_SORa::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{  
00152   if(!IsComputed_) return -1;
00153   Time_.ResetStartTime();
00154   bool initial_guess_is_zero=false;  
00155   int NumMyRows=W_->NumMyRows();
00156   int NumVectors = X.NumVectors();
00157 
00158   // need to create an auxiliary vector, Xcopy
00159   Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
00160   if (X.Pointers()[0] == Y.Pointers()[0]){
00161     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00162     // Since the user didn't give us anything better, our initial guess is zero.
00163     Y.Scale(0.0);
00164     initial_guess_is_zero=true;
00165   }
00166   else
00167     Xcopy = Teuchos::rcp( &X, false );
00168 
00169   Epetra_MultiVector Temp(Y);
00170   
00171   Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
00172   // Note: Assuming that the matrix has an importer.  Ifpack_PointRelaxation doesn't do this, but given that 
00173   // I have a CrsMatrix, I'm probably OK.
00174   if (IsParallel_)  Y2 = Teuchos::rcp( new Epetra_MultiVector(W_->Importer()->TargetMap(),NumVectors));
00175   else Y2 = Teuchos::rcp( &Y, false );
00176 
00177   // Pointer grabs
00178   int* rowptr,*colind;
00179   double *values;
00180   double **t_ptr,** y_ptr, ** y2_ptr, **x_ptr,*d_ptr;
00181   Y2->ExtractView(&y2_ptr);
00182   Y.ExtractView(&y_ptr);
00183   Temp.ExtractView(&t_ptr);
00184   Xcopy->ExtractView(&x_ptr);
00185   Wdiag_->ExtractView(&d_ptr);
00186   IFPACK_CHK_ERR(W_->ExtractCrsDataPointers(rowptr,colind,values));
00187 
00188   for(int i=0; i<NumSweeps_; i++){
00189     // Import Y2 if parallel
00190     if(IsParallel_)
00191       IFPACK_CHK_ERR(Y2->Import(Y,*W_->Importer(),Insert));
00192     
00193     // Calculate Ax 
00194     if(!initial_guess_is_zero  || i > 0) A_->Apply(Y,Temp);
00195     else Temp.Scale(0.0);
00196 
00197     // Do backsolve & update
00198     // x = x  + W^{-1} (b - A x)
00199     // NTS: This works since W_ has only the strict lower triangle and Wdiag_ has the diagonal
00200     for(int j=0; j<NumMyRows; j++){
00201       double diag=d_ptr[j];
00202       for (int m=0 ; m<NumVectors; m++) {
00203   double dtmp=0.0;
00204   for(int k=rowptr[j];k<rowptr[j+1];k++){
00205     dtmp+= values[k]*y2_ptr[m][colind[k]];
00206   }
00207   // Yes, we need to update both of these.
00208   y2_ptr[m][j] += (x_ptr[m][j] - t_ptr[m][j] - dtmp)/diag;     
00209   y_ptr[m][j] = y2_ptr[m][j];
00210       }
00211     }
00212   }
00213   
00214   // Counter update
00215   NumApplyInverse_++;
00216   ApplyInverseTime_ += Time_.ElapsedTime();
00217   return 0;
00218 }
00219 
00220 
00221 ostream& Ifpack_SORa::Print(ostream& os) const{
00222   os<<"Ifpack_SORa"<<endl;
00223   os<<" alpha = "<<Alpha_<<endl;
00224   os<<" gamma = "<<Gamma_<<endl;
00225   os<<endl;
00226   return os;
00227 }
00228 
00229 
00230 double Ifpack_SORa::Condest(const Ifpack_CondestType CT, 
00231                              const int MaxIters,
00232                              const double Tol,
00233                              Epetra_RowMatrix* Matrix_in){
00234   return -1.0;
00235 }
00236 
00237 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:35 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3