Ifpack_BlockRelaxation.h

00001 #ifndef IFPACK_BLOCKPRECONDITIONER_H
00002 #define IFPACK_BLOCKPRECONDITIONER_H
00003 
00004 #include "Ifpack_ConfigDefs.h"
00005 #include "Ifpack_Preconditioner.h" 
00006 #include "Ifpack_Partitioner.h"
00007 #include "Ifpack_LinearPartitioner.h"
00008 #include "Ifpack_GreedyPartitioner.h"
00009 #include "Ifpack_METISPartitioner.h"
00010 #include "Ifpack_EquationPartitioner.h"
00011 #include "Ifpack_UserPartitioner.h"
00012 #include "Ifpack_Graph_Epetra_RowMatrix.h"
00013 #include "Ifpack_DenseContainer.h" 
00014 #include "Ifpack_Utils.h" 
00015 #include "Teuchos_ParameterList.hpp"
00016 #include "Teuchos_RefCountPtr.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_in)
00147   {
00148     if (UseTranspose_in)
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_in = 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   Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_;
00371   mutable std::vector<Teuchos::RefCountPtr<T> > Containers_;
00373   Teuchos::RefCountPtr<Ifpack_Partitioner> Partitioner_;
00374   string PartitionerType_;
00375   int PrecType_;
00377   string Label_;
00379   bool ZeroStartingSolution_;
00380   Teuchos::RefCountPtr<Ifpack_Graph> Graph_;
00382   Teuchos::RefCountPtr<Epetra_Vector> W_;
00383   // Level of overlap among blocks (for Jacobi only).
00384   int OverlapLevel_;
00385   mutable Epetra_Time Time_;
00386   bool IsParallel_;
00387   Teuchos::RefCountPtr<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_in) :
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_(Teuchos::rcp(Matrix_in,false)),
00411   PartitionerType_("greedy"),
00412   PrecType_(IFPACK_JACOBI),
00413   ZeroStartingSolution_(true),
00414   OverlapLevel_(0),
00415   Time_(Comm()),
00416   IsParallel_(false)
00417 {
00418   if (Matrix_in->Comm().NumProc() != 1)
00419     IsParallel_ = true;
00420 }
00421 
00422 //==============================================================================
00423 template<typename T>
00424 Ifpack_BlockRelaxation<T>::~Ifpack_BlockRelaxation()
00425 {
00426 }
00427 
00428 //==============================================================================
00429 template<typename T>
00430 const char* Ifpack_BlockRelaxation<T>::Label() const
00431 {
00432   return(Label_.c_str());
00433 }
00434 
00435 //==============================================================================
00436 template<typename T>
00437 int Ifpack_BlockRelaxation<T>::
00438 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00439 {
00440   IFPACK_RETURN(Matrix().Apply(X,Y));
00441 }
00442 
00443 //==============================================================================
00444 template<typename T>
00445 const Epetra_Comm& Ifpack_BlockRelaxation<T>::
00446 Comm() const
00447 {
00448   return(Matrix().Comm());
00449 }
00450 
00451 //==============================================================================
00452 template<typename T>
00453 const Epetra_Map& Ifpack_BlockRelaxation<T>::
00454 OperatorDomainMap() const
00455 {
00456   return(Matrix().OperatorDomainMap());
00457 }
00458 
00459 //==============================================================================
00460 template<typename T>
00461 const Epetra_Map& Ifpack_BlockRelaxation<T>::
00462 OperatorRangeMap() const
00463 {
00464   return(Matrix().OperatorRangeMap());
00465 }
00466 
00467 //==============================================================================
00468 template<typename T>
00469 int Ifpack_BlockRelaxation<T>::ExtractSubmatrices()
00470 {
00471 
00472   if (Partitioner_ == Teuchos::null)
00473     IFPACK_CHK_ERR(-3);
00474 
00475   NumLocalBlocks_ = Partitioner_->NumLocalParts();
00476 
00477   Containers_.resize(NumLocalBlocks());
00478 
00479   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00480 
00481     int rows = Partitioner_->NumRowsInPart(i);
00482     Containers_[i] = Teuchos::rcp( new T(rows) );
00483     
00484     //Ifpack_DenseContainer* DC = 0;
00485     //DC = dynamic_cast<Ifpack_DenseContainer*>(Containers_[i]);
00486 
00487     if (Containers_[i] == Teuchos::null)
00488       IFPACK_CHK_ERR(-5);
00489     
00490     IFPACK_CHK_ERR(Containers_[i]->SetParameters(List_));
00491     IFPACK_CHK_ERR(Containers_[i]->Initialize());
00492     // flops in Initialize() will be computed on-the-fly in method InitializeFlops().
00493 
00494     // set "global" ID of each partitioner row
00495     for (int j = 0 ; j < rows ; ++j) {
00496       int LRID = (*Partitioner_)(i,j);
00497       Containers_[i]->ID(j) = LRID;
00498     }
00499 
00500     IFPACK_CHK_ERR(Containers_[i]->Compute(*Matrix_));
00501     // flops in Compute() will be computed on-the-fly in method ComputeFlops().
00502 
00503   }
00504 
00505   return(0);
00506 }
00507 
00508 //==============================================================================
00509 template<typename T>
00510 int Ifpack_BlockRelaxation<T>::Compute()
00511 {
00512 
00513   if (!IsInitialized())
00514     IFPACK_CHK_ERR(Initialize());
00515 
00516   Time_.ResetStartTime();
00517 
00518   IsComputed_ = false;
00519 
00520   if (Matrix().NumGlobalRows() != Matrix().NumGlobalCols())
00521     IFPACK_CHK_ERR(-2); // only square matrices
00522 
00523   IFPACK_CHK_ERR(ExtractSubmatrices());
00524   
00525   if (IsParallel_ && PrecType_ != IFPACK_JACOBI) {
00526     // not needed by Jacobi (done by matvec)
00527     Importer_ = Teuchos::rcp( new Epetra_Import(Matrix().RowMatrixColMap(),
00528                                                 Matrix().RowMatrixRowMap()) );
00529 
00530     if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
00531   }
00532   IsComputed_ = true;
00533   ComputeTime_ += Time_.ElapsedTime();
00534   ++NumCompute_;
00535 
00536   return(0);
00537 
00538 }
00539 
00540 //==============================================================================
00541 template<typename T>
00542 int Ifpack_BlockRelaxation<T>::
00543 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00544 {
00545   if (!IsComputed())
00546     IFPACK_CHK_ERR(-3);
00547 
00548   if (X.NumVectors() != Y.NumVectors())
00549     IFPACK_CHK_ERR(-2);
00550 
00551   Time_.ResetStartTime();
00552 
00553   // AztecOO gives X and Y pointing to the same memory location,
00554   // need to create an auxiliary vector, Xcopy
00555   Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
00556   if (X.Pointers()[0] == Y.Pointers()[0])
00557     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00558   else
00559     Xcopy = Teuchos::rcp( &X, false );
00560 
00561   switch (PrecType_) {
00562   case IFPACK_JACOBI:
00563     IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
00564     break;
00565   case IFPACK_GS:
00566     IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
00567     break;
00568   case IFPACK_SGS:
00569     IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
00570     break;
00571   }
00572 
00573   ApplyInverseTime_ += Time_.ElapsedTime();
00574   ++NumApplyInverse_;
00575 
00576   return(0);
00577 }
00578 
00579 //==============================================================================
00580 // This method in general will not work with AztecOO if used
00581 // outside Ifpack_AdditiveSchwarz and OverlapLevel_ != 0
00582 //
00583 template<typename T>
00584 int Ifpack_BlockRelaxation<T>::
00585 ApplyInverseJacobi(const Epetra_MultiVector& X, 
00586                    Epetra_MultiVector& Y) const
00587 {
00588 
00589   if (ZeroStartingSolution_)
00590     Y.PutScalar(0.0);
00591 
00592   // do not compute the residual in this case
00593   if (NumSweeps_ == 1 && ZeroStartingSolution_) {
00594     IFPACK_RETURN(DoJacobi(X,Y));
00595   }
00596 
00597   Epetra_MultiVector AX(Y);
00598 
00599   for (int j = 0; j < NumSweeps_ ; j++) {
00600     IFPACK_CHK_ERR(Apply(Y,AX));
00601     ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalNonzeros();
00602     IFPACK_CHK_ERR(AX.Update(1.0,X,-1.0));
00603     ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalRows();
00604     IFPACK_CHK_ERR(DoJacobi(AX,Y));
00605     // flops counted in DoJacobi()
00606   }
00607 
00608 
00609   return(0);
00610 }
00611 
00612 //==============================================================================
00613 template<typename T>
00614 int Ifpack_BlockRelaxation<T>::
00615 DoJacobi(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00616 {
00617   int NumVectors = X.NumVectors();
00618 
00619   if (OverlapLevel_ == 0) {
00620 
00621     for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00622      
00623       // may happen that a partition is empty
00624       if (Containers_[i]->NumRows() == 0) 
00625         continue;
00626 
00627       int LID;
00628 
00629       // extract RHS from X
00630       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00631         LID = Containers_[i]->ID(j);
00632         for (int k = 0 ; k < NumVectors ; ++k) {
00633           Containers_[i]->RHS(j,k) = X[k][LID];
00634         }
00635       }
00636 
00637       // apply the inverse of each block. NOTE: flops occurred
00638       // in ApplyInverse() of each block are summed up in method
00639       // ApplyInverseFlops().
00640       IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00641 
00642       // copy back into solution vector Y
00643       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00644         LID = Containers_[i]->ID(j);
00645         for (int k = 0 ; k < NumVectors ; ++k) {
00646           Y[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00647         }
00648       }
00649 
00650     }
00651     // NOTE: flops for ApplyInverse() of each block are summed up
00652     // in method ApplyInverseFlops()
00653     ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00654 
00655   }
00656   else {
00657 
00658     for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00659 
00660       // may happen that a partition is empty
00661       if (Containers_[i]->NumRows() == 0) 
00662         continue;
00663 
00664       int LID;
00665 
00666       // extract RHS from X
00667       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00668         LID = Containers_[i]->ID(j);
00669         for (int k = 0 ; k < NumVectors ; ++k) {
00670           Containers_[i]->RHS(j,k) = (*W_)[LID] * X[k][LID];
00671         }
00672       }
00673 
00674       // apply the inverse of each block
00675       IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00676 
00677       // copy back into solution vector Y
00678       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00679         LID = Containers_[i]->ID(j);
00680         for (int k = 0 ; k < NumVectors ; ++k) {
00681           Y[k][LID] += DampingFactor_ * (*W_)[LID] * Containers_[i]->LHS(j,k);
00682         }
00683       }
00684 
00685     }
00686     // NOTE: flops for ApplyInverse() of each block are summed up
00687     // in method ApplyInverseFlops()
00688     // NOTE: do not count for simplicity the flops due to overlapping rows
00689     ApplyInverseFlops_ += NumVectors * 4 * Matrix_->NumGlobalRows();
00690   }
00691 
00692   return(0);
00693 }
00694 
00695 //==============================================================================
00696 template<typename T>
00697 int Ifpack_BlockRelaxation<T>::
00698 ApplyInverseGS(const Epetra_MultiVector& X, 
00699                Epetra_MultiVector& Y) const
00700 {
00701 
00702   if (ZeroStartingSolution_)
00703     Y.PutScalar(0.0);
00704 
00705   Epetra_MultiVector Xcopy(X);
00706   for (int j = 0; j < NumSweeps_ ; j++) {
00707     IFPACK_CHK_ERR(DoGaussSeidel(Xcopy,Y));
00708     if (j != NumSweeps_ - 1)
00709       Xcopy = X;
00710   }
00711 
00712   return(0);
00713 
00714 }
00715 
00716 //==============================================================================
00717 template<typename T>
00718 int Ifpack_BlockRelaxation<T>::
00719 DoGaussSeidel(Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00720 {
00721 
00722   // cycle over all local subdomains
00723 
00724   int Length = Matrix().MaxNumEntries();
00725   std::vector<int> Indices(Length);
00726   std::vector<double> Values(Length);
00727 
00728   int NumMyRows = Matrix().NumMyRows();
00729   int NumVectors = X.NumVectors();
00730 
00731   // an additonal vector is needed by parallel computations
00732   // (note that applications through Ifpack_AdditiveSchwarz
00733   // are always seen are serial)
00734   Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
00735   if (IsParallel_)
00736     Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
00737   else
00738     Y2 = Teuchos::rcp( &Y, false );
00739 
00740   double** y_ptr;
00741   double** y2_ptr;
00742   Y.ExtractView(&y_ptr);
00743   Y2->ExtractView(&y2_ptr);
00744 
00745   // data exchange is here, once per sweep
00746   if (IsParallel_)
00747     IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00748 
00749   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00750 
00751     // may happen that a partition is empty
00752     if (Containers_[i]->NumRows() == 0) 
00753       continue;
00754 
00755     int LID;
00756 
00757     // update from previous block
00758 
00759     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00760       LID = Containers_[i]->ID(j);
00761 
00762       int NumEntries;
00763       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00764                                                &Values[0], &Indices[0]));
00765 
00766       for (int k = 0 ; k < NumEntries ; ++k) {
00767         int col = Indices[k];
00768 
00769           for (int kk = 0 ; kk < NumVectors ; ++kk) {
00770             X[kk][LID] -= Values[k] * y2_ptr[kk][col];
00771           }
00772       }
00773     }
00774 
00775     // solve with this block
00776 
00777     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00778       LID = Containers_[i]->ID(j);
00779       for (int k = 0 ; k < NumVectors ; ++k) {
00780         Containers_[i]->RHS(j,k) = X[k][LID];
00781       }
00782     }
00783 
00784     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00785     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00786 
00787     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00788       LID = Containers_[i]->ID(j);
00789       for (int k = 0 ; k < NumVectors ; ++k) {
00790         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00791       }
00792     }
00793 
00794   }
00795 
00796   // operations for all getrow()'s
00797   // NOTE: flops for ApplyInverse() of each block are summed up
00798   // in method ApplyInverseFlops()
00799   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00800   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00801 
00802   // Attention: this is delicate... Not all combinations
00803   // of Y2 and Y will always work (tough for ML it should be ok)
00804   if (IsParallel_)
00805     for (int m = 0 ; m < NumVectors ; ++m) 
00806       for (int i = 0 ; i < NumMyRows ; ++i)
00807         y_ptr[m][i] = y2_ptr[m][i];
00808 
00809   return(0);
00810 }
00811 
00812 //==============================================================================
00813 template<typename T>
00814 int Ifpack_BlockRelaxation<T>::
00815 ApplyInverseSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00816 {
00817 
00818   if (ZeroStartingSolution_)
00819     Y.PutScalar(0.0);
00820 
00821   Epetra_MultiVector Xcopy(X);
00822   for (int j = 0; j < NumSweeps_ ; j++) {
00823     IFPACK_CHK_ERR(DoSGS(X,Xcopy,Y));
00824     if (j != NumSweeps_ - 1)
00825       Xcopy = X;
00826   }
00827   return(0);
00828 }
00829 
00830 //==============================================================================
00831 template<typename T>
00832 int Ifpack_BlockRelaxation<T>::
00833 DoSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Xcopy, 
00834       Epetra_MultiVector& Y) const
00835 {
00836 
00837   int NumMyRows = Matrix().NumMyRows();
00838   int NumVectors = X.NumVectors();
00839 
00840   int Length = Matrix().MaxNumEntries();
00841   std::vector<int> Indices;
00842   std::vector<double> Values;
00843   Indices.resize(Length);
00844   Values.resize(Length);
00845 
00846   // an additonal vector is needed by parallel computations
00847   // (note that applications through Ifpack_AdditiveSchwarz
00848   // are always seen are serial)
00849   Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
00850   if (IsParallel_)
00851     Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
00852   else
00853     Y2 = Teuchos::rcp( &Y, false );
00854 
00855   double** y_ptr;
00856   double** y2_ptr;
00857   Y.ExtractView(&y_ptr);
00858   Y2->ExtractView(&y2_ptr);
00859 
00860   // data exchange is here, once per sweep
00861   if (IsParallel_)
00862     IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00863 
00864   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00865 
00866     // may happen that a partition is empty
00867     if (Containers_[i]->NumRows() == 0) 
00868       continue;
00869 
00870     int LID;
00871 
00872     // update from previous block
00873 
00874     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00875       LID = Containers_[i]->ID(j);
00876 
00877       int NumEntries;
00878       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00879                                                &Values[0], &Indices[0]));
00880 
00881       for (int k = 0 ; k < NumEntries ; ++k) {
00882         int col = Indices[k];
00883 
00884         for (int kk = 0 ; kk < NumVectors ; ++kk) {
00885           Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
00886         }
00887       }
00888     }
00889 
00890     // solve with this block
00891 
00892     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00893       LID = Containers_[i]->ID(j);
00894       for (int k = 0 ; k < NumVectors ; ++k) {
00895         Containers_[i]->RHS(j,k) = Xcopy[k][LID];
00896       }
00897     }
00898 
00899     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00900     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00901 
00902     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00903       LID = Containers_[i]->ID(j);
00904       for (int k = 0 ; k < NumVectors ; ++k) {
00905         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00906       }
00907     }
00908   }
00909 
00910   // operations for all getrow()'s
00911   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00912   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00913 
00914   Xcopy = X;
00915 
00916   for (int i = NumLocalBlocks() - 1; i >=0 ; --i) {
00917 
00918     if (Containers_[i]->NumRows() == 0) 
00919       continue;
00920 
00921     int LID;
00922 
00923     // update from previous block
00924 
00925     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00926       LID = Containers_[i]->ID(j);
00927 
00928       int NumEntries;
00929       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00930                                                &Values[0], &Indices[0]));
00931 
00932       for (int k = 0 ; k < NumEntries ; ++k) {
00933         int col = Indices[k];
00934 
00935           for (int kk = 0 ; kk < NumVectors ; ++kk) {
00936             Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
00937           }
00938       }
00939     }
00940 
00941     // solve with this block
00942 
00943     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00944       LID = Containers_[i]->ID(j);
00945       for (int k = 0 ; k < NumVectors ; ++k) {
00946         Containers_[i]->RHS(j,k) = Xcopy[k][LID];
00947       }
00948     }
00949 
00950     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00951     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00952 
00953     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00954       LID = Containers_[i]->ID(j);
00955       for (int k = 0 ; k < NumVectors ; ++k) {
00956         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00957       }
00958     }
00959   }
00960 
00961   // operations for all getrow()'s
00962   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00963   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00964 
00965   // Attention: this is delicate... Not all combinations
00966   // of Y2 and Y will always work (tough for ML it should be ok)
00967   if (IsParallel_)
00968     for (int m = 0 ; m < NumVectors ; ++m) 
00969       for (int i = 0 ; i < NumMyRows ; ++i)
00970         y_ptr[m][i] = y2_ptr[m][i];
00971 
00972   return(0);
00973 }
00974 
00975 //==============================================================================
00976 template<typename T>
00977 ostream& Ifpack_BlockRelaxation<T>::Print(ostream & os) const
00978 {
00979 
00980   string PT;
00981   if (PrecType_ == IFPACK_JACOBI)
00982     PT = "Jacobi";
00983   else if (PrecType_ == IFPACK_GS)
00984     PT = "Gauss-Seidel";
00985   else if (PrecType_ == IFPACK_SGS)
00986     PT = "symmetric Gauss-Seidel";
00987 
00988   if (!Comm().MyPID()) {
00989     os << endl;
00990     os << "================================================================================" << endl;
00991     os << "Ifpack_BlockRelaxation, " << PT << endl;
00992     os << "Sweeps = " << NumSweeps_ << endl;
00993     os << "Damping factor = " << DampingFactor_;
00994     if (ZeroStartingSolution_) 
00995       os << ", using zero starting solution" << endl;
00996     else
00997       os << ", using input starting solution" << endl;
00998     os << "Number of local blocks = " << Partitioner_->NumLocalParts() << endl;
00999     //os << "Condition number estimate = " << Condest_ << endl;
01000     os << "Global number of rows            = " << Matrix_->NumGlobalRows() << endl;
01001     os << endl;
01002     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
01003     os << "-----           -------   --------------       ------------     --------" << endl;
01004     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
01005        << "  " << std::setw(15) << InitializeTime() 
01006        << "  " << std::setw(15) << 1.0e-6 * InitializeFlops();
01007     if (InitializeTime() != 0.0)
01008       os << "  " << std::setw(15) << 1.0e-6 * InitializeFlops() / InitializeTime() << endl;
01009     else
01010       os << "  " << std::setw(15) << 0.0 << endl;
01011     os << "Compute()       "   << std::setw(5) << NumCompute() 
01012        << "  " << std::setw(15) << ComputeTime()
01013        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
01014     if (ComputeTime() != 0.0) 
01015       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
01016     else
01017       os << "  " << std::setw(15) << 0.0 << endl;
01018     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
01019        << "  " << std::setw(15) << ApplyInverseTime()
01020        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
01021     if (ApplyInverseTime() != 0.0) 
01022       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
01023     else
01024       os << "  " << std::setw(15) << 0.0 << endl;
01025     os << "================================================================================" << endl;
01026     os << endl;
01027   }
01028 
01029   return(os);
01030 }
01031 
01032 //==============================================================================
01033 template<typename T>
01034 int Ifpack_BlockRelaxation<T>::SetParameters(Teuchos::ParameterList& List)
01035 {
01036 
01037   string PT;
01038   if (PrecType_ == IFPACK_JACOBI)
01039     PT = "Jacobi";
01040   else if (PrecType_ == IFPACK_GS)
01041     PT = "Gauss-Seidel";
01042   else if (PrecType_ == IFPACK_SGS)
01043     PT = "symmetric Gauss-Seidel";
01044 
01045   PT = List.get("relaxation: type", PT);
01046 
01047   if (PT == "Jacobi") {
01048     PrecType_ = IFPACK_JACOBI;
01049   }
01050   else if (PT == "Gauss-Seidel") {
01051     PrecType_ = IFPACK_GS;
01052   }
01053   else if (PT == "symmetric Gauss-Seidel") {
01054     PrecType_ = IFPACK_SGS;
01055   } else {
01056     cerr << "Option `relaxation: type' has an incorrect value ("
01057       << PT << ")'" << endl;
01058     cerr << "(file " << __FILE__ << ", line " << __LINE__ << ")" << endl;
01059     exit(EXIT_FAILURE);
01060   }
01061 
01062   NumSweeps_            = List.get("relaxation: sweeps", NumSweeps_);
01063   DampingFactor_        = List.get("relaxation: damping factor", 
01064                                    DampingFactor_);
01065   ZeroStartingSolution_ = List.get("relaxation: zero starting solution", 
01066                                    ZeroStartingSolution_);
01067   PartitionerType_      = List.get("partitioner: type", 
01068                                    PartitionerType_);
01069   NumLocalBlocks_       = List.get("partitioner: local parts", 
01070                                    NumLocalBlocks_);
01071   // only Jacobi can work with overlap among local domains,
01072   OverlapLevel_         = List.get("partitioner: overlap", 
01073                                    OverlapLevel_);
01074 
01075   // check parameters
01076   if (PrecType_ != IFPACK_JACOBI)
01077     OverlapLevel_ = 0;
01078   if (NumLocalBlocks_ < 0)
01079     NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
01080   // other checks are performed in Partitioner_
01081   
01082   // copy the list as each subblock's constructor will
01083   // require it later
01084   List_ = List;
01085 
01086   // set the label
01087   string PT2;
01088   if (PrecType_ == IFPACK_JACOBI)
01089     PT2 = "BJ";
01090   else if (PrecType_ == IFPACK_GS)
01091     PT2 = "BGS";
01092   else if (PrecType_ == IFPACK_SGS)
01093     PT2 = "BSGS";
01094   Label_ = "IFPACK (" + PT2 + ", sweeps=" 
01095     + Ifpack_toString(NumSweeps_) + ", damping="
01096     + Ifpack_toString(DampingFactor_) + ", blocks="
01097     + Ifpack_toString(NumLocalBlocks()) + ")";
01098 
01099   return(0);
01100 }
01101 
01102 //==============================================================================
01103 template<typename T>
01104 int Ifpack_BlockRelaxation<T>::Initialize()
01105 {
01106   IsInitialized_ = false;
01107   Time_.ResetStartTime();
01108 
01109   Graph_ = Teuchos::rcp( new Ifpack_Graph_Epetra_RowMatrix(Teuchos::rcp(&Matrix(),false)) );
01110   if (Graph_ == Teuchos::null) IFPACK_CHK_ERR(-5);
01111 
01112   if (PartitionerType_ == "linear")
01113     Partitioner_ = Teuchos::rcp( new Ifpack_LinearPartitioner(&*Graph_) );
01114   else if (PartitionerType_ == "greedy")
01115     Partitioner_ = Teuchos::rcp( new Ifpack_GreedyPartitioner(&*Graph_) );
01116   else if (PartitionerType_ == "metis")
01117     Partitioner_ = Teuchos::rcp( new Ifpack_METISPartitioner(&*Graph_) );
01118   else if (PartitionerType_ == "equation")
01119     Partitioner_ = Teuchos::rcp( new Ifpack_EquationPartitioner(&*Graph_) );
01120   else if (PartitionerType_ == "user")
01121     Partitioner_ = Teuchos::rcp( new Ifpack_UserPartitioner(&*Graph_) );
01122   else
01123     IFPACK_CHK_ERR(-2);
01124 
01125   if (Partitioner_ == Teuchos::null) IFPACK_CHK_ERR(-5);
01126 
01127   // need to partition the graph of A
01128   IFPACK_CHK_ERR(Partitioner_->SetParameters(List_));
01129   IFPACK_CHK_ERR(Partitioner_->Compute());
01130 
01131   // get actual number of partitions
01132   NumLocalBlocks_ = Partitioner_->NumLocalParts();
01133 
01134   // weight of each vertex
01135   W_ = Teuchos::rcp( new Epetra_Vector(Matrix().RowMatrixRowMap()) );
01136   W_->PutScalar(0.0);
01137 
01138   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
01139 
01140     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
01141       int LID = (*Partitioner_)(i,j);
01142       (*W_)[LID]++;
01143     }
01144   }
01145   W_->Reciprocal(*W_);
01146 
01147   InitializeTime_ += Time_.ElapsedTime();
01148   IsInitialized_ = true;
01149   ++NumInitialize_;
01150 
01151   return(0);
01152 }
01153 
01154 //==============================================================================
01155 #endif // IFPACK_BLOCKPRECONDITIONER_H
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:22 2011 for IFPACK by  doxygen 1.6.3