Ifpack_BlockRelaxation.h

00001 #ifndef IFPACK_BLOCKPRECONDITIONER_H
00002 #define IFPACK_BLOCKPRECONDITIONER_H
00003 
00004 #include "Ifpack_ConfigDefs.h"
00005 #ifdef HAVE_IFPACK_TEUCHOS
00006 #include "Ifpack_Preconditioner.h" 
00007 #include "Ifpack_Partitioner.h"
00008 #include "Ifpack_LinearPartitioner.h"
00009 #include "Ifpack_GreedyPartitioner.h"
00010 #include "Ifpack_METISPartitioner.h"
00011 #include "Ifpack_EquationPartitioner.h"
00012 #include "Ifpack_UserPartitioner.h"
00013 #include "Ifpack_Graph_Epetra_RowMatrix.h"
00014 #include "Ifpack_DenseContainer.h" 
00015 #include "Ifpack_Utils.h" 
00016 #include "Teuchos_ParameterList.hpp"
00017 #include "Epetra_RowMatrix.h"
00018 #include "Epetra_MultiVector.h"
00019 #include "Epetra_Vector.h"
00020 #include "Epetra_Time.h"
00021 #include "Epetra_Import.h"
00022 
00023 static const int IFPACK_JACOBI = 0;
00024 static const int IFPACK_GS = 1;
00025 static const int IFPACK_SGS = 2;
00026 
00027 
00029 
00092 template<typename T>
00093 class Ifpack_BlockRelaxation : public Ifpack_Preconditioner {
00094 
00095 public:
00096 
00098 
00099 
00104   Ifpack_BlockRelaxation(const Epetra_RowMatrix* Matrix);
00105 
00106   virtual ~Ifpack_BlockRelaxation();
00107 
00109 
00111 
00113 
00121   virtual int Apply(const Epetra_MultiVector& X, 
00122             Epetra_MultiVector& Y) const;
00123 
00125 
00134   virtual int ApplyInverse(const Epetra_MultiVector& X, 
00135                Epetra_MultiVector& Y) const;
00136 
00138   virtual double NormInf() const
00139   {
00140     return(-1.0);
00141   }
00143 
00145 
00146   virtual int SetUseTranspose(bool UseTranspose)
00147   {
00148     if (UseTranspose)
00149       IFPACK_CHK_ERR(-98); // FIXME: can I work with the transpose?
00150     return(0);
00151   }
00152 
00153   virtual const char* Label() const;
00154  
00156   virtual bool UseTranspose() const
00157   {
00158     return(false);
00159   }
00160 
00162   virtual bool HasNormInf() const
00163   {
00164     return(false);
00165   }
00166 
00168   virtual const Epetra_Comm & Comm() const;
00169 
00171   virtual const Epetra_Map & OperatorDomainMap() const;
00172 
00174   virtual const Epetra_Map & OperatorRangeMap() const;
00176 
00178   int NumLocalBlocks() const 
00179   {
00180     return(NumLocalBlocks_);
00181   }
00182 
00184   virtual bool IsInitialized() const
00185   {
00186     return(IsInitialized_);
00187   }
00188 
00190   virtual bool IsComputed() const
00191   {
00192     return(IsComputed_);
00193   }
00194 
00196   virtual int SetParameters(Teuchos::ParameterList& List);
00197 
00199   virtual int Initialize();
00200 
00202   virtual int Compute();
00203 
00204   virtual const Epetra_RowMatrix& Matrix() const
00205   {
00206     return(*Matrix_);
00207   }
00208 
00209   virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
00210                          const int MaxIters = 1550,
00211                          const double Tol = 1e-9,
00212              Epetra_RowMatrix* Matrix = 0)
00213   {
00214     return(-1.0);
00215   }
00216 
00217   virtual double Condest() const
00218   {
00219     return(-1.0);
00220   }
00221 
00222   std::ostream& Print(std::ostream& os) const;
00223 
00225   virtual int NumInitialize() const
00226   {
00227     return(NumInitialize_);
00228   }
00229 
00231   virtual int NumCompute() const
00232   {
00233     return(NumCompute_);
00234   }
00235 
00237   virtual int NumApplyInverse() const
00238   {
00239     return(NumApplyInverse_);
00240   }
00241 
00243   virtual double InitializeTime() const
00244   {
00245     return(InitializeTime_);
00246   }
00247 
00249   virtual double ComputeTime() const
00250   {
00251     return(ComputeTime_);
00252   }
00253 
00255   virtual double ApplyInverseTime() const
00256   {
00257     return(ApplyInverseTime_);
00258   }
00259 
00261   virtual double InitializeFlops() const
00262   {
00263     if (Containers_.size() == 0)
00264       return(0.0);
00265 
00266     // the total number of flops is computed each time InitializeFlops() is
00267     // called. This is becase I also have to add the contribution from each
00268     // container.
00269     double total = InitializeFlops_;
00270     for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
00271       total += Containers_[i]->InitializeFlops();
00272     return(total);
00273   }
00274 
00275   virtual double ComputeFlops() const
00276   {
00277     if (Containers_.size() == 0)
00278       return(0.0);
00279     
00280     double total = ComputeFlops_;
00281     for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
00282       total += Containers_[i]->ComputeFlops();
00283     return(total);
00284   }
00285 
00286   virtual double ApplyInverseFlops() const
00287   {
00288     if (Containers_.size() == 0)
00289       return(0.0);
00290 
00291     double total = ApplyInverseFlops_;
00292     for (unsigned int i = 0 ; i < Containers_.size() ; ++i) {
00293       total += Containers_[i]->ApplyInverseFlops();
00294     }
00295     return(total);
00296   }
00297 
00298 private:
00299 
00301   Ifpack_BlockRelaxation(const Ifpack_BlockRelaxation& rhs);
00302 
00304   Ifpack_BlockRelaxation & operator=(const Ifpack_BlockRelaxation& rhs)
00305   {
00306     return(*this);
00307   }
00308 
00309   virtual int ApplyInverseJacobi(const Epetra_MultiVector& X, 
00310                                  Epetra_MultiVector& Y) const;
00311 
00312   virtual int DoJacobi(const Epetra_MultiVector& X, 
00313                                   Epetra_MultiVector& Y) const;
00314 
00315   virtual int ApplyInverseGS(const Epetra_MultiVector& X, 
00316                              Epetra_MultiVector& Y) const;
00317 
00318   virtual int DoGaussSeidel(Epetra_MultiVector& X, 
00319                             Epetra_MultiVector& Y) const;
00320 
00321   virtual int ApplyInverseSGS(const Epetra_MultiVector& X, 
00322                               Epetra_MultiVector& Y) const;
00323 
00324   virtual int DoSGS(const Epetra_MultiVector& X,
00325                     Epetra_MultiVector& Xtmp,
00326                     Epetra_MultiVector& Y) const;
00327 
00328   int ExtractSubmatrices();
00329 
00330   // @{ Initializations, timing and flops
00331 
00333   bool IsInitialized_;
00335   bool IsComputed_;
00337   int NumInitialize_;
00339   int NumCompute_;
00341   mutable int NumApplyInverse_;
00343   double InitializeTime_;
00345   double ComputeTime_;
00347   mutable double ApplyInverseTime_;
00349   double InitializeFlops_;
00351   double ComputeFlops_;
00353   mutable double ApplyInverseFlops_;
00354   // @}
00355 
00356   // @{ Settings
00358   int NumSweeps_;
00360   double DampingFactor_;
00362   int NumLocalBlocks_;
00364   Teuchos::ParameterList List_;
00365   // @}
00366 
00367   // @{ Other data
00370   const Epetra_RowMatrix* Matrix_;
00371   mutable vector<T*> Containers_;
00373   Ifpack_Partitioner* Partitioner_;
00374   string PartitionerType_;
00375   int PrecType_;
00377   string Label_;
00379   bool ZeroStartingSolution_;
00380   Ifpack_Graph* Graph_;
00382   Epetra_Vector* W_;
00383   // Level of overlap among blocks (for Jacobi only).
00384   int OverlapLevel_;
00385   mutable Epetra_Time Time_;
00386   bool IsParallel_;
00387   Epetra_Import* Importer_;
00388   // @}
00389 
00390 }; // class Ifpack_BlockRelaxation
00391 
00392 //==============================================================================
00393 template<typename T>
00394 Ifpack_BlockRelaxation<T>::
00395 Ifpack_BlockRelaxation(const Epetra_RowMatrix* Matrix) :
00396   IsInitialized_(false),
00397   IsComputed_(false),
00398   NumInitialize_(0),
00399   NumCompute_(0),
00400   NumApplyInverse_(0),
00401   InitializeTime_(0.0),
00402   ComputeTime_(0.0),
00403   ApplyInverseTime_(0.0),
00404   InitializeFlops_(0.0),
00405   ComputeFlops_(0.0),
00406   ApplyInverseFlops_(0.0),
00407   NumSweeps_(1),
00408   DampingFactor_(1.0),
00409   NumLocalBlocks_(1),
00410   Matrix_(Matrix),
00411   Containers_(0),
00412   Partitioner_(0),
00413   PartitionerType_("greedy"),
00414   PrecType_(IFPACK_JACOBI),
00415   ZeroStartingSolution_(true),
00416   Graph_(0),
00417   W_(0),
00418   OverlapLevel_(0),
00419   Time_(Comm()),
00420   IsParallel_(false),
00421   Importer_(0)
00422 {
00423   if (Matrix->Comm().NumProc() != 1)
00424     IsParallel_ = true;
00425 }
00426 
00427 //==============================================================================
00428 template<typename T>
00429 Ifpack_BlockRelaxation<T>::~Ifpack_BlockRelaxation()
00430 {
00431   for (int i = 0 ; i < NumLocalBlocks() ; ++i)
00432     delete Containers_[i];
00433   if (Partitioner_)
00434     delete Partitioner_;
00435   if (Graph_)
00436     delete Graph_;
00437   if (W_)
00438     delete W_;
00439   if (Importer_)
00440     delete Importer_;
00441 }
00442 
00443 //==============================================================================
00444 template<typename T>
00445 const char* Ifpack_BlockRelaxation<T>::Label() const
00446 {
00447   return(Label_.c_str());
00448 }
00449 
00450 //==============================================================================
00451 template<typename T>
00452 int Ifpack_BlockRelaxation<T>::
00453 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00454 {
00455   IFPACK_RETURN(Matrix().Apply(X,Y));
00456 }
00457 
00458 //==============================================================================
00459 template<typename T>
00460 const Epetra_Comm& Ifpack_BlockRelaxation<T>::
00461 Comm() const
00462 {
00463   return(Matrix().Comm());
00464 }
00465 
00466 //==============================================================================
00467 template<typename T>
00468 const Epetra_Map& Ifpack_BlockRelaxation<T>::
00469 OperatorDomainMap() const
00470 {
00471   return(Matrix().OperatorDomainMap());
00472 }
00473 
00474 //==============================================================================
00475 template<typename T>
00476 const Epetra_Map& Ifpack_BlockRelaxation<T>::
00477 OperatorRangeMap() const
00478 {
00479   return(Matrix().OperatorRangeMap());
00480 }
00481 
00482 //==============================================================================
00483 template<typename T>
00484 int Ifpack_BlockRelaxation<T>::ExtractSubmatrices()
00485 {
00486 
00487   if (Partitioner_ == 0)
00488     IFPACK_CHK_ERR(-3);
00489 
00490   NumLocalBlocks_ = Partitioner_->NumLocalParts();
00491 
00492   Containers_.resize(NumLocalBlocks());
00493 
00494   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00495 
00496     int rows = Partitioner_->NumRowsInPart(i);
00497     Containers_[i] = new T(rows);
00498     
00499     //Ifpack_DenseContainer* DC = 0;
00500     //DC = dynamic_cast<Ifpack_DenseContainer*>(Containers_[i]);
00501 
00502     if (Containers_[i] == 0)
00503       IFPACK_CHK_ERR(-5);
00504     
00505     IFPACK_CHK_ERR(Containers_[i]->SetParameters(List_));
00506     IFPACK_CHK_ERR(Containers_[i]->Initialize());
00507     // flops in Initialize() will be computed on-the-fly in method InitializeFlops().
00508 
00509     // set "global" ID of each partitioner row
00510     for (int j = 0 ; j < rows ; ++j) {
00511       int LRID = (*Partitioner_)(i,j);
00512       Containers_[i]->ID(j) = LRID;
00513     }
00514 
00515     IFPACK_CHK_ERR(Containers_[i]->Compute(*Matrix_));
00516     // flops in Compute() will be computed on-the-fly in method ComputeFlops().
00517 
00518   }
00519 
00520   return(0);
00521 }
00522 
00523 //==============================================================================
00524 template<typename T>
00525 int Ifpack_BlockRelaxation<T>::Compute()
00526 {
00527 
00528   if (!IsInitialized())
00529     IFPACK_CHK_ERR(Initialize());
00530 
00531   Time_.ResetStartTime();
00532 
00533   IsComputed_ = false;
00534 
00535   if (Matrix().NumGlobalRows() != Matrix().NumGlobalCols())
00536     IFPACK_CHK_ERR(-2); // only square matrices
00537 
00538   IFPACK_CHK_ERR(ExtractSubmatrices());
00539   
00540   if (IsParallel_ && PrecType_ != IFPACK_JACOBI) {
00541     // not needed by Jacobi (done by matvec)
00542     Importer_ = new Epetra_Import(Matrix().RowMatrixColMap(),
00543                                   Matrix().RowMatrixRowMap());
00544 
00545     if (Importer_ == 0) IFPACK_CHK_ERR(-5);
00546   }
00547   IsComputed_ = true;
00548   ComputeTime_ += Time_.ElapsedTime();
00549   ++NumCompute_;
00550 
00551   return(0);
00552 
00553 }
00554 
00555 //==============================================================================
00556 template<typename T>
00557 int Ifpack_BlockRelaxation<T>::
00558 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00559 {
00560   if (!IsComputed())
00561     IFPACK_CHK_ERR(-3);
00562 
00563   if (X.NumVectors() != Y.NumVectors())
00564     IFPACK_CHK_ERR(-2);
00565 
00566   Time_.ResetStartTime();
00567 
00568   // AztecOO gives X and Y pointing to the same memory location,
00569   // need to create an auxiliary vector, Xcopy
00570   const Epetra_MultiVector* Xcopy;
00571   if (X.Pointers()[0] == Y.Pointers()[0])
00572     Xcopy = new Epetra_MultiVector(X);
00573   else
00574     Xcopy = &X;
00575 
00576   switch (PrecType_) {
00577   case IFPACK_JACOBI:
00578     IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
00579     break;
00580   case IFPACK_GS:
00581     IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
00582     break;
00583   case IFPACK_SGS:
00584     IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
00585     break;
00586   }
00587 
00588   if (Xcopy != &X)
00589     delete Xcopy;
00590 
00591   ApplyInverseTime_ += Time_.ElapsedTime();
00592   ++NumApplyInverse_;
00593 
00594   return(0);
00595 }
00596 
00597 //==============================================================================
00598 // This method in general will not work with AztecOO if used
00599 // outside Ifpack_AdditiveSchwarz and OverlapLevel_ != 0
00600 //
00601 template<typename T>
00602 int Ifpack_BlockRelaxation<T>::
00603 ApplyInverseJacobi(const Epetra_MultiVector& X, 
00604                    Epetra_MultiVector& Y) const
00605 {
00606 
00607   if (ZeroStartingSolution_)
00608     Y.PutScalar(0.0);
00609 
00610   // do not compute the residual in this case
00611   if (NumSweeps_ == 1 && ZeroStartingSolution_) {
00612     IFPACK_RETURN(DoJacobi(X,Y));
00613   }
00614 
00615   Epetra_MultiVector AX(Y);
00616 
00617   for (int j = 0; j < NumSweeps_ ; j++) {
00618     IFPACK_CHK_ERR(Apply(Y,AX));
00619     ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalNonzeros();
00620     IFPACK_CHK_ERR(AX.Update(1.0,X,-1.0));
00621     ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalRows();
00622     IFPACK_CHK_ERR(DoJacobi(AX,Y));
00623     // flops counted in DoJacobi()
00624   }
00625 
00626 
00627   return(0);
00628 }
00629 
00630 //==============================================================================
00631 template<typename T>
00632 int Ifpack_BlockRelaxation<T>::
00633 DoJacobi(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00634 {
00635   int NumVectors = X.NumVectors();
00636 
00637   if (OverlapLevel_ == 0) {
00638 
00639     for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00640      
00641       // may happen that a partition is empty
00642       if (Containers_[i]->NumRows() == 0) 
00643         continue;
00644 
00645       int LID;
00646 
00647       // extract RHS from X
00648       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00649         LID = Containers_[i]->ID(j);
00650         for (int k = 0 ; k < NumVectors ; ++k) {
00651           Containers_[i]->RHS(j,k) = X[k][LID];
00652         }
00653       }
00654 
00655       // apply the inverse of each block. NOTE: flops occurred
00656       // in ApplyInverse() of each block are summed up in method
00657       // ApplyInverseFlops().
00658       IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00659 
00660       // copy back into solution vector Y
00661       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00662         LID = Containers_[i]->ID(j);
00663         for (int k = 0 ; k < NumVectors ; ++k) {
00664           Y[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00665         }
00666       }
00667 
00668     }
00669     // NOTE: flops for ApplyInverse() of each block are summed up
00670     // in method ApplyInverseFlops()
00671     ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00672 
00673   }
00674   else {
00675 
00676     for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00677 
00678       // may happen that a partition is empty
00679       if (Containers_[i]->NumRows() == 0) 
00680         continue;
00681 
00682       int LID;
00683 
00684       // extract RHS from X
00685       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00686         LID = Containers_[i]->ID(j);
00687         for (int k = 0 ; k < NumVectors ; ++k) {
00688           Containers_[i]->RHS(j,k) = (*W_)[LID] * X[k][LID];
00689         }
00690       }
00691 
00692       // apply the inverse of each block
00693       IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00694 
00695       // copy back into solution vector Y
00696       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00697         LID = Containers_[i]->ID(j);
00698         for (int k = 0 ; k < NumVectors ; ++k) {
00699           Y[k][LID] += DampingFactor_ * (*W_)[LID] * Containers_[i]->LHS(j,k);
00700         }
00701       }
00702 
00703     }
00704     // NOTE: flops for ApplyInverse() of each block are summed up
00705     // in method ApplyInverseFlops()
00706     // NOTE: do not count for simplicity the flops due to overlapping rows
00707     ApplyInverseFlops_ += NumVectors * 4 * Matrix_->NumGlobalRows();
00708   }
00709 
00710   return(0);
00711 }
00712 
00713 //==============================================================================
00714 template<typename T>
00715 int Ifpack_BlockRelaxation<T>::
00716 ApplyInverseGS(const Epetra_MultiVector& X, 
00717                Epetra_MultiVector& Y) const
00718 {
00719 
00720   if (ZeroStartingSolution_)
00721     Y.PutScalar(0.0);
00722 
00723   Epetra_MultiVector Xcopy(X);
00724   for (int j = 0; j < NumSweeps_ ; j++) {
00725     IFPACK_CHK_ERR(DoGaussSeidel(Xcopy,Y));
00726     if (j != NumSweeps_ - 1)
00727       Xcopy = X;
00728   }
00729 
00730   return(0);
00731 
00732 }
00733 
00734 //==============================================================================
00735 template<typename T>
00736 int Ifpack_BlockRelaxation<T>::
00737 DoGaussSeidel(Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00738 {
00739 
00740   // cycle over all local subdomains
00741 
00742   int Length = Matrix().MaxNumEntries();
00743   vector<int> Indices(Length);
00744   vector<double> Values(Length);
00745 
00746   int NumMyRows = Matrix().NumMyRows();
00747   int NumVectors = X.NumVectors();
00748 
00749   // an additonal vector is needed by parallel computations
00750   // (note that applications through Ifpack_AdditiveSchwarz
00751   // are always seen are serial)
00752   Epetra_MultiVector* Y2;
00753   if (IsParallel_)
00754     Y2 = new Epetra_MultiVector(Importer_->TargetMap(), NumVectors);
00755   else
00756     Y2 = &Y;
00757 
00758   double** y_ptr;
00759   double** y2_ptr;
00760   Y.ExtractView(&y_ptr);
00761   Y2->ExtractView(&y2_ptr);
00762 
00763   // data exchange is here, once per sweep
00764   if (IsParallel_)
00765     IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00766 
00767   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00768 
00769     // may happen that a partition is empty
00770     if (Containers_[i]->NumRows() == 0) 
00771       continue;
00772 
00773     int LID;
00774 
00775     // update from previous block
00776 
00777     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00778       LID = Containers_[i]->ID(j);
00779 
00780       int NumEntries;
00781       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00782                                                &Values[0], &Indices[0]));
00783 
00784       for (int k = 0 ; k < NumEntries ; ++k) {
00785         int col = Indices[k];
00786 
00787           for (int kk = 0 ; kk < NumVectors ; ++kk) {
00788             X[kk][LID] -= Values[k] * y2_ptr[kk][col];
00789           }
00790       }
00791     }
00792 
00793     // solve with this block
00794 
00795     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00796       LID = Containers_[i]->ID(j);
00797       for (int k = 0 ; k < NumVectors ; ++k) {
00798         Containers_[i]->RHS(j,k) = X[k][LID];
00799       }
00800     }
00801 
00802     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00803     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00804 
00805     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00806       LID = Containers_[i]->ID(j);
00807       for (int k = 0 ; k < NumVectors ; ++k) {
00808         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00809       }
00810     }
00811 
00812   }
00813 
00814   // operations for all getrow()'s
00815   // NOTE: flops for ApplyInverse() of each block are summed up
00816   // in method ApplyInverseFlops()
00817   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00818   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00819 
00820   // Attention: this is delicate... Not all combinations
00821   // of Y2 and Y will always work (tough for ML it should be ok)
00822   if (IsParallel_)
00823     for (int m = 0 ; m < NumVectors ; ++m) 
00824       for (int i = 0 ; i < NumMyRows ; ++i)
00825         y_ptr[m][i] = y2_ptr[m][i];
00826 
00827   if (IsParallel_)
00828     delete Y2;
00829 
00830   return(0);
00831 }
00832 
00833 //==============================================================================
00834 template<typename T>
00835 int Ifpack_BlockRelaxation<T>::
00836 ApplyInverseSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00837 {
00838 
00839   if (ZeroStartingSolution_)
00840     Y.PutScalar(0.0);
00841 
00842   Epetra_MultiVector Xcopy(X);
00843   for (int j = 0; j < NumSweeps_ ; j++) {
00844     IFPACK_CHK_ERR(DoSGS(X,Xcopy,Y));
00845     if (j != NumSweeps_ - 1)
00846       Xcopy = X;
00847   }
00848   return(0);
00849 }
00850 
00851 //==============================================================================
00852 template<typename T>
00853 int Ifpack_BlockRelaxation<T>::
00854 DoSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Xcopy, 
00855       Epetra_MultiVector& Y) const
00856 {
00857 
00858   int NumMyRows = Matrix().NumMyRows();
00859   int NumVectors = X.NumVectors();
00860 
00861   int Length = Matrix().MaxNumEntries();
00862   vector<int> Indices;
00863   vector<double> Values;
00864   Indices.resize(Length);
00865   Values.resize(Length);
00866 
00867   // an additonal vector is needed by parallel computations
00868   // (note that applications through Ifpack_AdditiveSchwarz
00869   // are always seen are serial)
00870   Epetra_MultiVector* Y2;
00871   if (IsParallel_)
00872     Y2 = new Epetra_MultiVector(Importer_->TargetMap(), NumVectors);
00873   else
00874     Y2 = &Y;
00875 
00876   double** y_ptr;
00877   double** y2_ptr;
00878   Y.ExtractView(&y_ptr);
00879   Y2->ExtractView(&y2_ptr);
00880 
00881   // data exchange is here, once per sweep
00882   if (IsParallel_)
00883     IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00884 
00885   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00886 
00887     // may happen that a partition is empty
00888     if (Containers_[i]->NumRows() == 0) 
00889       continue;
00890 
00891     int LID;
00892 
00893     // update from previous block
00894 
00895     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00896       LID = Containers_[i]->ID(j);
00897 
00898       int NumEntries;
00899       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00900                                                &Values[0], &Indices[0]));
00901 
00902       for (int k = 0 ; k < NumEntries ; ++k) {
00903         int col = Indices[k];
00904 
00905         for (int kk = 0 ; kk < NumVectors ; ++kk) {
00906           Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
00907         }
00908       }
00909     }
00910 
00911     // solve with this block
00912 
00913     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00914       LID = Containers_[i]->ID(j);
00915       for (int k = 0 ; k < NumVectors ; ++k) {
00916         Containers_[i]->RHS(j,k) = Xcopy[k][LID];
00917       }
00918     }
00919 
00920     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00921     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00922 
00923     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00924       LID = Containers_[i]->ID(j);
00925       for (int k = 0 ; k < NumVectors ; ++k) {
00926         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00927       }
00928     }
00929   }
00930 
00931   // operations for all getrow()'s
00932   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00933   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00934 
00935   Xcopy = X;
00936 
00937   for (int i = NumLocalBlocks() - 1; i >=0 ; --i) {
00938 
00939     if (Containers_[i]->NumRows() == 0) 
00940       continue;
00941 
00942     int LID;
00943 
00944     // update from previous block
00945 
00946     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00947       LID = Containers_[i]->ID(j);
00948 
00949       int NumEntries;
00950       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00951                                                &Values[0], &Indices[0]));
00952 
00953       for (int k = 0 ; k < NumEntries ; ++k) {
00954         int col = Indices[k];
00955 
00956           for (int kk = 0 ; kk < NumVectors ; ++kk) {
00957             Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
00958           }
00959       }
00960     }
00961 
00962     // solve with this block
00963 
00964     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00965       LID = Containers_[i]->ID(j);
00966       for (int k = 0 ; k < NumVectors ; ++k) {
00967         Containers_[i]->RHS(j,k) = Xcopy[k][LID];
00968       }
00969     }
00970 
00971     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00972     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00973 
00974     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00975       LID = Containers_[i]->ID(j);
00976       for (int k = 0 ; k < NumVectors ; ++k) {
00977         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00978       }
00979     }
00980   }
00981 
00982   // operations for all getrow()'s
00983   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00984   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00985 
00986   // Attention: this is delicate... Not all combinations
00987   // of Y2 and Y will always work (tough for ML it should be ok)
00988   if (IsParallel_)
00989     for (int m = 0 ; m < NumVectors ; ++m) 
00990       for (int i = 0 ; i < NumMyRows ; ++i)
00991         y_ptr[m][i] = y2_ptr[m][i];
00992 
00993   if (IsParallel_)
00994     delete Y2;
00995 
00996   return(0);
00997 }
00998 
00999 //==============================================================================
01000 template<typename T>
01001 ostream& Ifpack_BlockRelaxation<T>::Print(ostream & os) const
01002 {
01003 
01004   string PT;
01005   if (PrecType_ == IFPACK_JACOBI)
01006     PT = "Jacobi";
01007   else if (PrecType_ == IFPACK_GS)
01008     PT = "Gauss-Seidel";
01009   else if (PrecType_ == IFPACK_SGS)
01010     PT = "symmetric Gauss-Seidel";
01011 
01012   if (!Comm().MyPID()) {
01013     os << endl;
01014     os << "================================================================================" << endl;
01015     os << "Ifpack_BlockRelaxation, " << PT << endl;
01016     os << "Sweeps = " << NumSweeps_ << endl;
01017     os << "Damping factor = " << DampingFactor_;
01018     if (ZeroStartingSolution_) 
01019       os << ", using zero starting solution" << endl;
01020     else
01021       os << ", using input starting solution" << endl;
01022     os << "Number of local blocks = " << Partitioner_->NumLocalParts() << endl;
01023     //os << "Condition number estimate = " << Condest_ << endl;
01024     os << "Global number of rows            = " << Matrix_->NumGlobalRows() << endl;
01025     os << endl;
01026     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
01027     os << "-----           -------   --------------       ------------     --------" << endl;
01028     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
01029        << "  " << std::setw(15) << InitializeTime() 
01030        << "  " << std::setw(15) << 1.0e-6 * InitializeFlops();
01031     if (InitializeTime() != 0.0)
01032       os << "  " << std::setw(15) << 1.0e-6 * InitializeFlops() / InitializeTime() << endl;
01033     else
01034       os << "  " << std::setw(15) << 0.0 << endl;
01035     os << "Compute()       "   << std::setw(5) << NumCompute() 
01036        << "  " << std::setw(15) << ComputeTime()
01037        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
01038     if (ComputeTime() != 0.0) 
01039       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
01040     else
01041       os << "  " << std::setw(15) << 0.0 << endl;
01042     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
01043        << "  " << std::setw(15) << ApplyInverseTime()
01044        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
01045     if (ApplyInverseTime() != 0.0) 
01046       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
01047     else
01048       os << "  " << std::setw(15) << 0.0 << endl;
01049     os << "================================================================================" << endl;
01050     os << endl;
01051   }
01052 
01053   return(os);
01054 }
01055 
01056 //==============================================================================
01057 template<typename T>
01058 int Ifpack_BlockRelaxation<T>::SetParameters(Teuchos::ParameterList& List)
01059 {
01060 
01061   string PT;
01062   if (PrecType_ == IFPACK_JACOBI)
01063     PT = "Jacobi";
01064   else if (PrecType_ == IFPACK_GS)
01065     PT = "Gauss-Seidel";
01066   else if (PrecType_ == IFPACK_SGS)
01067     PT = "symmetric Gauss-Seidel";
01068 
01069   PT = List.get("relaxation: type", PT);
01070 
01071   if (PT == "Jacobi") {
01072     PrecType_ = IFPACK_JACOBI;
01073   }
01074   else if (PT == "Gauss-Seidel") {
01075     PrecType_ = IFPACK_GS;
01076   }
01077   else if (PT == "symmetric Gauss-Seidel") {
01078     PrecType_ = IFPACK_SGS;
01079   } else {
01080     cerr << "Option `relaxation: type' has an incorrect value ("
01081       << PT << ")'" << endl;
01082     cerr << "(file " << __FILE__ << ", line " << __LINE__ << ")" << endl;
01083     exit(EXIT_FAILURE);
01084   }
01085 
01086   NumSweeps_            = List.get("relaxation: sweeps", NumSweeps_);
01087   DampingFactor_        = List.get("relaxation: damping factor", 
01088                                    DampingFactor_);
01089   ZeroStartingSolution_ = List.get("relaxation: zero starting solution", 
01090                                    ZeroStartingSolution_);
01091   PartitionerType_      = List.get("partitioner: type", 
01092                                    PartitionerType_);
01093   NumLocalBlocks_       = List.get("partitioner: local parts", 
01094                                    NumLocalBlocks_);
01095   // only Jacobi can work with overlap among local domains,
01096   OverlapLevel_         = List.get("partitioner: overlap", 
01097                                    OverlapLevel_);
01098 
01099   // check parameters
01100   if (PrecType_ != IFPACK_JACOBI)
01101     OverlapLevel_ = 0;
01102   if (NumLocalBlocks_ < 0)
01103     NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
01104   // other checks are performed in Partitioner_
01105   
01106   // copy the list as each subblock's constructor will
01107   // require it later
01108   List_ = List;
01109 
01110   // set the label
01111   string PT2;
01112   if (PrecType_ == IFPACK_JACOBI)
01113     PT2 = "BJ";
01114   else if (PrecType_ == IFPACK_GS)
01115     PT2 = "BGS";
01116   else if (PrecType_ == IFPACK_SGS)
01117     PT2 = "BSGS";
01118   Label_ = "IFPACK (" + PT2 + ", sweeps=" 
01119     + Ifpack_toString(NumSweeps_) + ", damping="
01120     + Ifpack_toString(DampingFactor_) + ", blocks="
01121     + Ifpack_toString(NumLocalBlocks()) + ")";
01122 
01123   return(0);
01124 }
01125 
01126 //==============================================================================
01127 template<typename T>
01128 int Ifpack_BlockRelaxation<T>::Initialize()
01129 {
01130   IsInitialized_ = false;
01131   Time_.ResetStartTime();
01132 
01133   if (Partitioner_)
01134     delete Partitioner_;
01135   if (Graph_)
01136     delete Graph_;
01137 
01138   Graph_ = new Ifpack_Graph_Epetra_RowMatrix(&Matrix());
01139   if (Graph_ == 0) IFPACK_CHK_ERR(-5);
01140 
01141   if (PartitionerType_ == "linear")
01142     Partitioner_ = new Ifpack_LinearPartitioner(Graph_);
01143   else if (PartitionerType_ == "greedy")
01144     Partitioner_ = new Ifpack_GreedyPartitioner(Graph_);
01145   else if (PartitionerType_ == "metis")
01146     Partitioner_ = new Ifpack_METISPartitioner(Graph_);
01147   else if (PartitionerType_ == "equation")
01148     Partitioner_ = new Ifpack_EquationPartitioner(Graph_);
01149   else if (PartitionerType_ == "user")
01150     Partitioner_ = new Ifpack_UserPartitioner(Graph_);
01151   else
01152     IFPACK_CHK_ERR(-2);
01153 
01154   if (Partitioner_ == 0) IFPACK_CHK_ERR(-5);
01155 
01156   // need to partition the graph of A
01157   IFPACK_CHK_ERR(Partitioner_->SetParameters(List_));
01158   IFPACK_CHK_ERR(Partitioner_->Compute());
01159 
01160   // get actual number of partitions
01161   NumLocalBlocks_ = Partitioner_->NumLocalParts();
01162 
01163   // weight of each vertex
01164   if (W_)
01165     delete W_;
01166   W_ = new Epetra_Vector(Matrix().RowMatrixRowMap());
01167   W_->PutScalar(0.0);
01168 
01169   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
01170 
01171     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
01172       int LID = (*Partitioner_)(i,j);
01173       (*W_)[LID]++;
01174     }
01175   }
01176   W_->Reciprocal(*W_);
01177 
01178   InitializeTime_ += Time_.ElapsedTime();
01179   IsInitialized_ = true;
01180   ++NumInitialize_;
01181 
01182   return(0);
01183 }
01184 
01185 //==============================================================================
01186 #endif // HAVE_IFPACK_TEUCHOS
01187 #endif // IFPACK_BLOCKPRECONDITIONER_H

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