Ifpack_PointRelaxation.cpp

00001 #include "Ifpack_ConfigDefs.h"
00002 #ifdef HAVE_IFPACK_TEUCHOS
00003 #include <iomanip>
00004 #include "Epetra_Operator.h"
00005 #include "Epetra_RowMatrix.h"
00006 #include "Epetra_CrsMatrix.h"
00007 #include "Epetra_Comm.h"
00008 #include "Epetra_Map.h"
00009 #include "Epetra_MultiVector.h"
00010 #include "Epetra_Vector.h"
00011 #include "Epetra_Time.h"
00012 #include "Epetra_Import.h"
00013 #include "Ifpack_Preconditioner.h"
00014 #include "Ifpack_PointRelaxation.h"
00015 #include "Ifpack_Utils.h"
00016 #include "Ifpack_Condest.h"
00017 
00018 static const int IFPACK_JACOBI = 0;
00019 static const int IFPACK_GS = 1;
00020 static const int IFPACK_SGS = 2;
00021 
00022 //==============================================================================
00023 Ifpack_PointRelaxation::
00024 Ifpack_PointRelaxation(const Epetra_RowMatrix* Matrix) :
00025   IsInitialized_(false),
00026   IsComputed_(false),
00027   NumInitialize_(0),
00028   NumCompute_(0),
00029   NumApplyInverse_(0),
00030   InitializeTime_(0.0),
00031   ComputeTime_(0.0),
00032   ApplyInverseTime_(0.0),
00033   ComputeFlops_(0.0),
00034   ApplyInverseFlops_(0.0),
00035   NumSweeps_(1),
00036   DampingFactor_(1.0),
00037   UseTranspose_(false),
00038   Condest_(-1.0),
00039   ComputeCondest_(false),
00040   PrecType_(IFPACK_JACOBI),
00041   MinDiagonalValue_(0.0),
00042   NumMyRows_(0),
00043   NumMyNonzeros_(0),
00044   NumGlobalRows_(0),
00045   NumGlobalNonzeros_(0),
00046   Matrix_(Matrix),
00047   Importer_(0),
00048   Diagonal_(0),
00049   Time_(0),
00050   IsParallel_(false),
00051   ZeroStartingSolution_(true)
00052 {
00053 }
00054 
00055 //==============================================================================
00056 Ifpack_PointRelaxation::~Ifpack_PointRelaxation()
00057 {
00058   if (Diagonal_)
00059     delete Diagonal_;
00060   if (Time_)
00061     delete Time_;
00062   if (Importer_)
00063     delete Importer_;
00064 }
00065 
00066 //==============================================================================
00067 int Ifpack_PointRelaxation::SetParameters(Teuchos::ParameterList& List)
00068 {
00069 
00070   string PT;
00071   if (PrecType_ == IFPACK_JACOBI)
00072     PT = "Jacobi";
00073   else if (PrecType_ == IFPACK_GS)
00074     PT = "Gauss-Seidel";
00075   else if (PrecType_ == IFPACK_SGS)
00076     PT = "symmetric Gauss-Seidel";
00077 
00078   PT = List.get("relaxation: type", PT);
00079 
00080   if (PT == "Jacobi")
00081     PrecType_ = IFPACK_JACOBI;
00082   else if (PT == "Gauss-Seidel")
00083     PrecType_ = IFPACK_GS;
00084   else if (PT == "symmetric Gauss-Seidel")
00085     PrecType_ = IFPACK_SGS;
00086   else {
00087     IFPACK_CHK_ERR(-2);
00088   }
00089   
00090   NumSweeps_            = List.get("relaxation: sweeps",NumSweeps_);
00091   DampingFactor_        = List.get("relaxation: damping factor", 
00092                                    DampingFactor_);
00093   MinDiagonalValue_     = List.get("relaxation: min diagonal value", 
00094                                    MinDiagonalValue_);
00095   ZeroStartingSolution_ = List.get("relaxation: zero starting solution", 
00096                                    ZeroStartingSolution_);
00097 
00098   SetLabel();
00099 
00100   return(0);
00101 }
00102 
00103 //==============================================================================
00104 const Epetra_Comm& Ifpack_PointRelaxation::Comm() const
00105 {
00106   return(Matrix_->Comm());
00107 }
00108 
00109 //==============================================================================
00110 const Epetra_Map& Ifpack_PointRelaxation::OperatorDomainMap() const
00111 {
00112   return(Matrix_->OperatorDomainMap());
00113 }
00114 
00115 //==============================================================================
00116 const Epetra_Map& Ifpack_PointRelaxation::OperatorRangeMap() const
00117 {
00118   return(Matrix_->OperatorRangeMap());
00119 }
00120 
00121 //==============================================================================
00122 int Ifpack_PointRelaxation::
00123 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00124 {
00125 
00126   if (IsComputed() == false)
00127     IFPACK_CHK_ERR(-3);
00128 
00129   if (X.NumVectors() != Y.NumVectors())
00130     IFPACK_CHK_ERR(-2);
00131 
00132   IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
00133   return(0);
00134 }
00135 
00136 //==============================================================================
00137 int Ifpack_PointRelaxation::Initialize()
00138 {
00139   IsInitialized_ = false;
00140 
00141   if (Matrix_ == 0)
00142     IFPACK_CHK_ERR(-2);
00143 
00144   if (Time_ == 0)
00145     Time_ = new Epetra_Time(Comm());
00146 
00147   if (Matrix().NumGlobalRows() != Matrix().NumGlobalCols())
00148     IFPACK_CHK_ERR(-2); // only square matrices
00149 
00150   NumMyRows_ = Matrix_->NumMyRows();
00151   NumMyNonzeros_ = Matrix_->NumMyNonzeros();
00152   NumGlobalRows_ = Matrix_->NumGlobalRows();
00153   NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros();
00154 
00155   if (Comm().NumProc() != 1)
00156     IsParallel_ = true;
00157   else
00158     IsParallel_ = false;
00159 
00160   ++NumInitialize_;
00161   InitializeTime_ += Time_->ElapsedTime();
00162   IsInitialized_ = true;
00163   return(0);
00164 }
00165 
00166 //==============================================================================
00167 int Ifpack_PointRelaxation::Compute()
00168 {
00169   if (!IsInitialized())
00170     IFPACK_CHK_ERR(Initialize());
00171 
00172   Time_->ResetStartTime();
00173 
00174   // reset values
00175   IsComputed_ = false;
00176   Condest_ = -1.0;
00177 
00178   if (NumSweeps_ <= 0)
00179     IFPACK_CHK_ERR(-2); // at least one application
00180   
00181   if (Diagonal_)
00182     delete Diagonal_;
00183   Diagonal_ = new Epetra_Vector(Matrix().RowMatrixRowMap());
00184 
00185   if (Diagonal_ == 0)
00186     IFPACK_CHK_ERR(-5);
00187 
00188   IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*Diagonal_));
00189 
00190   // check diagonal elements, store the inverses, and verify that
00191   // no zeros are around. If an element is zero, then by default
00192   // its inverse is zero as well (that is, the row is ignored).
00193   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00194     double& diag = (*Diagonal_)[i];
00195     if (IFPACK_ABS(diag) < MinDiagonalValue_)
00196       diag = MinDiagonalValue_;
00197     if (diag != 0.0)
00198       diag = 1.0 / diag;
00199   }
00200   ComputeFlops_ += NumMyRows_;
00201 
00202 #if 0
00203   // some methods require the inverse of the diagonal, compute it
00204   // now and store it in `Diagonal_'
00205   if ((PrecType_ == IFPACK_JACOBI) || (PrecType_ == IFPACK_GS)) {
00206     Diagonal_->Reciprocal(*Diagonal_);
00207     // update flops
00208     ComputeFlops_ += NumMyRows_;
00209   }
00210 #endif
00211 
00212   // We need to import data from external processors. Here I create an
00213   // Epetra_Import object because I cannot assume that Matrix_ has one.
00214   // This is a bit of waste of resources (but the code is more robust).
00215   // Note that I am doing some strange stuff to set the components of Y
00216   // from Y2 (to save some time).
00217   //
00218   if (IsParallel_ && ((PrecType_ == IFPACK_GS) || (PrecType_ == IFPACK_SGS))) {
00219     Importer_ = new Epetra_Import(Matrix().RowMatrixColMap(),
00220                                   Matrix().RowMatrixRowMap());
00221     if (Importer_ == 0) IFPACK_CHK_ERR(-5);
00222   }
00223 
00224   ++NumCompute_;
00225   ComputeTime_ += Time_->ElapsedTime();
00226   IsComputed_ = true;
00227 
00228   return(0);
00229 }
00230 
00231 //==============================================================================
00232 ostream& Ifpack_PointRelaxation::Print(ostream & os) const
00233 {
00234 
00235   double MyMinVal, MyMaxVal;
00236   double MinVal, MaxVal;
00237 
00238   if (IsComputed_) {
00239     Diagonal_->MinValue(&MyMinVal);
00240     Diagonal_->MaxValue(&MyMaxVal);
00241     Comm().MinAll(&MyMinVal,&MinVal,1);
00242     Comm().MinAll(&MyMaxVal,&MaxVal,1);
00243   }
00244 
00245   if (!Comm().MyPID()) {
00246     os << endl;
00247     os << "================================================================================" << endl;
00248     os << "Ifpack_PointRelaxation" << endl;
00249     os << "Sweeps         = " << NumSweeps_ << endl;
00250     os << "damping factor = " << DampingFactor_ << endl;
00251     if (PrecType_ == IFPACK_JACOBI)
00252       os << "Type           = Jacobi" << endl;
00253     else if (PrecType_ == IFPACK_GS)
00254       os << "Type           = Gauss-Seidel" << endl;
00255     else if (PrecType_ == IFPACK_SGS)
00256       os << "Type           = symmetric Gauss-Seidel" << endl;
00257     if (ZeroStartingSolution_) 
00258       os << "Using zero starting solution" << endl;
00259     else
00260       os << "Using input starting solution" << endl;
00261     os << "Condition number estimate = " << Condest() << endl;
00262     os << "Global number of rows            = " << Matrix_->NumGlobalRows() << endl;
00263     if (IsComputed_) {
00264       os << "Minimum value on stored diagonal = " << MinVal << endl;
00265       os << "Maximum value on stored diagonal = " << MaxVal << endl;
00266     }
00267     os << endl;
00268     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00269     os << "-----           -------   --------------       ------------     --------" << endl;
00270     os << "Initialize()    "   << std::setw(5) << NumInitialize_ 
00271        << "  " << std::setw(15) << InitializeTime_ 
00272        << "              0.0              0.0" << endl;
00273     os << "Compute()       "   << std::setw(5) << NumCompute_ 
00274        << "  " << std::setw(15) << ComputeTime_
00275        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_;
00276     if (ComputeTime_ != 0.0)
00277       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
00278     else
00279       os << "  " << std::setw(15) << 0.0 << endl;
00280     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse_ 
00281        << "  " << std::setw(15) << ApplyInverseTime_
00282        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
00283     if (ApplyInverseTime_ != 0.0)
00284       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
00285     else
00286       os << "  " << std::setw(15) << 0.0 << endl;
00287     os << "================================================================================" << endl;
00288     os << endl;
00289   }
00290 
00291   return(os);
00292 }
00293 
00294 //==============================================================================
00295 double Ifpack_PointRelaxation::
00296 Condest(const Ifpack_CondestType CT, 
00297         const int MaxIters, const double Tol,
00298     Epetra_RowMatrix* Matrix)
00299 {
00300   if (!IsComputed()) // cannot compute right now
00301     return(-1.0);
00302 
00303   // always computes it. Call Condest() with no parameters to get
00304   // the previous estimate.
00305   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);
00306 
00307   return(Condest_);
00308 }
00309 
00310 //==============================================================================
00311 void Ifpack_PointRelaxation::SetLabel()
00312 {
00313   string PT;
00314   if (PrecType_ == IFPACK_JACOBI)
00315     PT = "Jacobi";
00316   else if (PrecType_ == IFPACK_GS)
00317     PT = "GS";
00318   else if (PrecType_ == IFPACK_SGS)
00319     PT = "SGS";
00320 
00321   Label_ = "IFPACK (" + PT + ", sweeps=" + Ifpack_toString(NumSweeps_)
00322     + ", damping=" + Ifpack_toString(DampingFactor_) + ")";
00323 }
00324 
00325 //==============================================================================
00326 // Note that Ifpack_PointRelaxation and Jacobi is much faster than
00327 // Ifpack_AdditiveSchwarz<Ifpack_PointRelaxation> (because of the
00328 // way the matrix-vector produce is performed).
00329 //
00330 // Another ML-related observation is that the starting solution (in Y)
00331 // is NOT supposed to be zero. This may slow down the application of just
00332 // one sweep of Jacobi.
00333 //
00334 int Ifpack_PointRelaxation::
00335 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00336 {
00337   if (!IsComputed())
00338     IFPACK_CHK_ERR(-3);
00339 
00340   if (X.NumVectors() != Y.NumVectors())
00341     IFPACK_CHK_ERR(-2);
00342 
00343   Time_->ResetStartTime();
00344 
00345   // AztecOO gives X and Y pointing to the same memory location,
00346   // need to create an auxiliary vector, Xcopy
00347   const Epetra_MultiVector* Xcopy;
00348   if (X.Pointers()[0] == Y.Pointers()[0])
00349     Xcopy = new Epetra_MultiVector(X);
00350   else
00351     Xcopy = &X;
00352     
00353   if (ZeroStartingSolution_)
00354     Y.PutScalar(0.0);
00355 
00356   // Flops are updated in each of the following. 
00357   switch (PrecType_) {
00358   case IFPACK_JACOBI:
00359     IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
00360     break;
00361   case IFPACK_GS:
00362     IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
00363     break;
00364   case IFPACK_SGS:
00365     IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
00366     break;
00367   default:
00368     IFPACK_CHK_ERR(-1); // something wrong
00369   }
00370 
00371   if (Xcopy != &X)
00372     delete Xcopy;
00373 
00374   ++NumApplyInverse_;
00375   ApplyInverseTime_ += Time_->ElapsedTime();
00376   return(0);
00377 }
00378 
00379 //==============================================================================
00380 // This preconditioner can be much slower than AztecOO and ML versions
00381 // if the matrix-vector product requires all ExtractMyRowCopy() 
00382 // (as done through Ifpack_AdditiveSchwarz).
00383 int Ifpack_PointRelaxation::
00384 ApplyInverseJacobi(const Epetra_MultiVector& RHS, Epetra_MultiVector& LHS) const
00385 {
00386 
00387   int NumVectors = LHS.NumVectors();
00388   Epetra_MultiVector* A_times_LHS;
00389   A_times_LHS = new Epetra_MultiVector(LHS.Map(),NumVectors);
00390   if (A_times_LHS == 0) IFPACK_CHK_ERR(-5);
00391 
00392   for (int j = 0; j < NumSweeps_ ; j++) {
00393 
00394     IFPACK_CHK_ERR(Apply(LHS,*A_times_LHS));
00395     IFPACK_CHK_ERR(A_times_LHS->Update(1.0,RHS,-1.0));
00396     for (int v = 0 ; v < NumVectors ; ++v)
00397       IFPACK_CHK_ERR(LHS(v)->Multiply(DampingFactor_, *(*A_times_LHS)(v), 
00398                                      *Diagonal_, 1.0));
00399 
00400   }
00401   delete A_times_LHS;
00402 
00403   // Flops:
00404   // - matrix vector              (2 * NumGlobalNonzeros_)
00405   // - update                     (2 * NumGlobalRows_)
00406   // - Multiply:
00407   //   - DampingFactor            (NumGlobalRows_)
00408   //   - Diagonal                 (NumGlobalRows_)
00409   //   - A + B                    (NumGlobalRows_)
00410   //   - 1.0                      (NumGlobalRows_)
00411   ApplyInverseFlops_ += NumVectors * (6 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
00412 
00413   return(0);
00414 }
00415 
00416 //==============================================================================
00417 int Ifpack_PointRelaxation::
00418 ApplyInverseGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00419 {
00420   int NumVectors = X.NumVectors();
00421 
00422   int Length = Matrix().MaxNumEntries();
00423   vector<int> Indices(Length);
00424   vector<double> Values(Length);
00425 
00426   Epetra_MultiVector* Y2;
00427   if (IsParallel_)
00428     Y2 = new Epetra_MultiVector(Importer_->TargetMap(), NumVectors);
00429   else
00430     Y2 = &Y;
00431                         
00432   // extract views (for nicer and faster code)
00433   double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
00434   X.ExtractView(&x_ptr);
00435   Y.ExtractView(&y_ptr);
00436   Y2->ExtractView(&y2_ptr);
00437   Diagonal_->ExtractView(&d_ptr);
00438 
00439   for (int j = 0; j < NumSweeps_ ; j++) {
00440 
00441     // data exchange is here, once per sweep
00442     if (IsParallel_)
00443       IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00444 
00445     // FIXME: do I really need this code below?
00446     if (NumVectors == 1) {
00447 
00448       double* y0_ptr = y_ptr[0];
00449       double* y20_ptr = y2_ptr[0];
00450       double* x0_ptr = x_ptr[0];
00451 
00452       for (int i = 0 ; i < NumMyRows_ ; ++i) {
00453 
00454         int NumEntries;
00455         int col;
00456         IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
00457                                                  &Values[0], &Indices[0]));
00458 
00459         double dtemp = 0.0;
00460         for (int k = 0 ; k < NumEntries ; ++k) {
00461 
00462           col = Indices[k];
00463           dtemp += Values[k] * y20_ptr[col];
00464         }
00465 
00466         y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
00467       }
00468       // using Export() sounded quite expensive
00469       if (IsParallel_)
00470         for (int i = 0 ; i < NumMyRows_ ; ++i)
00471           y0_ptr[i] = y20_ptr[i];
00472 
00473     }
00474     else {
00475 
00476       for (int i = 0 ; i < NumMyRows_ ; ++i) {
00477 
00478         int NumEntries;
00479         int col;
00480         IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
00481                                                  &Values[0], &Indices[0]));
00482 
00483         for (int m = 0 ; m < NumVectors ; ++m) {
00484 
00485           double dtemp = 0.0;
00486           for (int k = 0 ; k < NumEntries ; ++k) {
00487 
00488             col = Indices[k];
00489             dtemp += Values[k] * y2_ptr[m][col];
00490           }
00491 
00492           y2_ptr[m][i] += DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
00493         }
00494         // using Export() sounded quite expensive
00495       }
00496 
00497       if (IsParallel_)
00498         for (int m = 0 ; m < NumVectors ; ++m) 
00499           for (int i = 0 ; i < NumMyRows_ ; ++i)
00500             y_ptr[m][i] = y2_ptr[m][i];
00501 
00502     }
00503   }
00504 
00505   if (IsParallel_)
00506     delete Y2;
00507 
00508   ApplyInverseFlops_ += NumVectors * (4 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
00509 
00510   return(0);
00511 }
00512 
00513 //==============================================================================
00514 int Ifpack_PointRelaxation::
00515 ApplyInverseSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00516 {
00517   const Epetra_CrsMatrix* CrsMatrix = dynamic_cast<const Epetra_CrsMatrix*>(Matrix_);
00518   // try to pick the best option; performances may be improved
00519   // if several sweeps are used.
00520   if (CrsMatrix != 0)
00521   {
00522     if (CrsMatrix->StorageOptimized())
00523       return(ApplyInverseSGS_FastCrsMatrix(CrsMatrix, X, Y));
00524     else
00525       return(ApplyInverseSGS_CrsMatrix(CrsMatrix, X, Y));
00526   }
00527   else
00528     return(ApplyInverseSGS_RowMatrix(X, Y));
00529 }
00530 
00531 //==============================================================================
00532 int Ifpack_PointRelaxation::
00533 ApplyInverseSGS_RowMatrix(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00534 {
00535   int NumVectors = X.NumVectors();
00536   int Length = Matrix().MaxNumEntries();
00537   vector<int> Indices(Length);
00538   vector<double> Values(Length);
00539 
00540   Epetra_MultiVector* Y2;
00541   if (IsParallel_) {
00542     Y2 = new Epetra_MultiVector(Importer_->TargetMap(), NumVectors);
00543   }
00544   else
00545     Y2 = &Y;
00546 
00547   double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
00548   X.ExtractView(&x_ptr);
00549   Y.ExtractView(&y_ptr);
00550   Y2->ExtractView(&y2_ptr);
00551   Diagonal_->ExtractView(&d_ptr);
00552   
00553   for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
00554     
00555     // only one data exchange per sweep
00556     if (IsParallel_)
00557       IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00558 
00559     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00560 
00561       int NumEntries;
00562       int col;
00563       double diag = d_ptr[i];
00564 
00565       IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
00566                                                &Values[0], &Indices[0]));
00567 
00568       for (int m = 0 ; m < NumVectors ; ++m) {
00569 
00570         double dtemp = 0.0;
00571 
00572         for (int k = 0 ; k < NumEntries ; ++k) {
00573 
00574           col = Indices[k];
00575           dtemp += Values[k] * y2_ptr[m][col];
00576         }
00577 
00578         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00579       }
00580     }
00581 
00582     for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {
00583 
00584       int NumEntries;
00585       int col;
00586       double diag = d_ptr[i];
00587 
00588       IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
00589                                                &Values[0], &Indices[0]));
00590 
00591       for (int m = 0 ; m < NumVectors ; ++m) {
00592 
00593         double dtemp = 0.0;
00594         for (int k = 0 ; k < NumEntries ; ++k) {
00595 
00596           col = Indices[k];
00597           dtemp += Values[k] * y2_ptr[m][col];
00598         }
00599 
00600         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00601 
00602       }
00603     }
00604 
00605     if (IsParallel_)
00606       for (int m = 0 ; m < NumVectors ; ++m) 
00607         for (int i = 0 ; i < NumMyRows_ ; ++i)
00608           y_ptr[m][i] = y2_ptr[m][i];
00609   }
00610 
00611   if (IsParallel_)
00612     delete Y2;
00613 
00614   ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
00615   return(0);
00616 }
00617 
00618 //==============================================================================
00619 int Ifpack_PointRelaxation::
00620 ApplyInverseSGS_CrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00621 {
00622   int NumVectors = X.NumVectors();
00623 
00624   int* Indices;
00625   double* Values;
00626 
00627   Epetra_MultiVector* Y2;
00628   if (IsParallel_) {
00629     Y2 = new Epetra_MultiVector(Importer_->TargetMap(), NumVectors);
00630   }
00631   else
00632     Y2 = &Y;
00633 
00634   double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
00635   X.ExtractView(&x_ptr);
00636   Y.ExtractView(&y_ptr);
00637   Y2->ExtractView(&y2_ptr);
00638   Diagonal_->ExtractView(&d_ptr);
00639   
00640   for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
00641     
00642     // only one data exchange per sweep
00643     if (IsParallel_)
00644       IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00645 
00646     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00647 
00648       int NumEntries;
00649       int col;
00650       double diag = d_ptr[i];
00651 
00652       IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
00653 
00654       for (int m = 0 ; m < NumVectors ; ++m) {
00655 
00656         double dtemp = 0.0;
00657 
00658         for (int k = 0; k < NumEntries; ++k) {
00659 
00660           col = Indices[k];
00661           dtemp += Values[k] * y2_ptr[m][col];
00662         }
00663 
00664         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00665       }
00666     }
00667 
00668     for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {
00669 
00670       int NumEntries;
00671       int col;
00672       double diag = d_ptr[i];
00673 
00674       IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
00675 
00676       for (int m = 0 ; m < NumVectors ; ++m) {
00677 
00678         double dtemp = 0.0;
00679         for (int k = 0; k < NumEntries; ++k) {
00680 
00681           col = Indices[k];
00682           dtemp += Values[k] * y2_ptr[m][col];
00683         }
00684 
00685         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00686 
00687       }
00688     }
00689 
00690     if (IsParallel_)
00691       for (int m = 0 ; m < NumVectors ; ++m) 
00692         for (int i = 0 ; i < NumMyRows_ ; ++i)
00693           y_ptr[m][i] = y2_ptr[m][i];
00694   }
00695 
00696   if (IsParallel_)
00697     delete Y2;
00698 
00699   ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
00700   return(0);
00701 }
00702 
00703 //==============================================================================
00704 // this requires Epetra_CrsMatrix + OptimizeStorage()
00705 int Ifpack_PointRelaxation::
00706 ApplyInverseSGS_FastCrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00707 {
00708   int* IndexOffset;
00709   int* Indices;
00710   double* Values;
00711   IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
00712 
00713   int NumVectors = X.NumVectors();
00714 
00715   Epetra_MultiVector* Y2;
00716   if (IsParallel_) {
00717     Y2 = new Epetra_MultiVector(Importer_->TargetMap(), NumVectors);
00718   }
00719   else
00720     Y2 = &Y;
00721 
00722   double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
00723   X.ExtractView(&x_ptr);
00724   Y.ExtractView(&y_ptr);
00725   Y2->ExtractView(&y2_ptr);
00726   Diagonal_->ExtractView(&d_ptr);
00727   
00728   for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
00729     
00730     // only one data exchange per sweep
00731     if (IsParallel_)
00732       IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00733 
00734     for (int i = 0 ; i < NumMyRows_ ; ++i) {
00735 
00736       int col;
00737       double diag = d_ptr[i];
00738 
00739       // no need to extract row i
00740       //IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
00741 
00742       for (int m = 0 ; m < NumVectors ; ++m) {
00743 
00744         double dtemp = 0.0;
00745 
00746         for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
00747 
00748           col = Indices[k];
00749           dtemp += Values[k] * y2_ptr[m][col];
00750         }
00751 
00752         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00753       }
00754     }
00755 
00756     for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {
00757 
00758       int col;
00759       double diag = d_ptr[i];
00760 
00761       // no need to extract row i
00762       //IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
00763 
00764       for (int m = 0 ; m < NumVectors ; ++m) {
00765 
00766         double dtemp = 0.0;
00767         for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
00768 
00769           col = Indices[k];
00770           dtemp += Values[k] * y2_ptr[m][col];
00771         }
00772 
00773         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00774 
00775       }
00776     }
00777 
00778     if (IsParallel_)
00779       for (int m = 0 ; m < NumVectors ; ++m) 
00780         for (int i = 0 ; i < NumMyRows_ ; ++i)
00781           y_ptr[m][i] = y2_ptr[m][i];
00782   }
00783 
00784   if (IsParallel_)
00785     delete Y2;
00786 
00787   ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
00788   return(0);
00789 }
00790 
00791 #endif

Generated on Thu Sep 18 12:37:08 2008 for IFPACK by doxygen 1.3.9.1