IFPACK Development
Ifpack_BlockRelaxation.h
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00006 //                 Copyright (2002) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 #ifndef IFPACK_BLOCKPRECONDITIONER_H
00045 #define IFPACK_BLOCKPRECONDITIONER_H
00046 
00047 #include "Ifpack_ConfigDefs.h"
00048 #include "Ifpack_Preconditioner.h" 
00049 #include "Ifpack_Partitioner.h"
00050 #include "Ifpack_LinePartitioner.h"
00051 #include "Ifpack_LinearPartitioner.h"
00052 #include "Ifpack_GreedyPartitioner.h"
00053 #include "Ifpack_METISPartitioner.h"
00054 #include "Ifpack_EquationPartitioner.h"
00055 #include "Ifpack_UserPartitioner.h"
00056 #include "Ifpack_Graph_Epetra_RowMatrix.h"
00057 #include "Ifpack_DenseContainer.h" 
00058 #include "Ifpack_Utils.h" 
00059 #include "Teuchos_ParameterList.hpp"
00060 #include "Teuchos_RefCountPtr.hpp"
00061 #include "Epetra_RowMatrix.h"
00062 #include "Epetra_MultiVector.h"
00063 #include "Epetra_Vector.h"
00064 #include "Epetra_Time.h"
00065 #include "Epetra_Import.h"
00066 
00067 static const int IFPACK_JACOBI = 0;
00068 static const int IFPACK_GS = 1;
00069 static const int IFPACK_SGS = 2;
00070 
00071 
00073 
00136 template<typename T>
00137 class Ifpack_BlockRelaxation : public Ifpack_Preconditioner {
00138 
00139 public:
00140 
00142 
00143 
00148   Ifpack_BlockRelaxation(const Epetra_RowMatrix* Matrix);
00149 
00150   virtual ~Ifpack_BlockRelaxation();
00151 
00153 
00155 
00157 
00165   virtual int Apply(const Epetra_MultiVector& X, 
00166             Epetra_MultiVector& Y) const;
00167 
00169 
00178   virtual int ApplyInverse(const Epetra_MultiVector& X, 
00179                Epetra_MultiVector& Y) const;
00180 
00182   virtual double NormInf() const
00183   {
00184     return(-1.0);
00185   }
00187 
00189 
00190   virtual int SetUseTranspose(bool UseTranspose_in)
00191   {
00192     if (UseTranspose_in)
00193       IFPACK_CHK_ERR(-98); // FIXME: can I work with the transpose?
00194     return(0);
00195   }
00196 
00197   virtual const char* Label() const;
00198  
00200   virtual bool UseTranspose() const
00201   {
00202     return(false);
00203   }
00204 
00206   virtual bool HasNormInf() const
00207   {
00208     return(false);
00209   }
00210 
00212   virtual const Epetra_Comm & Comm() const;
00213 
00215   virtual const Epetra_Map & OperatorDomainMap() const;
00216 
00218   virtual const Epetra_Map & OperatorRangeMap() const;
00220 
00222   int NumLocalBlocks() const 
00223   {
00224     return(NumLocalBlocks_);
00225   }
00226 
00228   virtual bool IsInitialized() const
00229   {
00230     return(IsInitialized_);
00231   }
00232 
00234   virtual bool IsComputed() const
00235   {
00236     return(IsComputed_);
00237   }
00238 
00240   virtual int SetParameters(Teuchos::ParameterList& List);
00241 
00243   virtual int Initialize();
00244 
00246   virtual int Compute();
00247 
00248   virtual const Epetra_RowMatrix& Matrix() const
00249   {
00250     return(*Matrix_);
00251   }
00252 
00253   virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
00254                          const int MaxIters = 1550,
00255                          const double Tol = 1e-9,
00256              Epetra_RowMatrix* Matrix_in = 0)
00257   {
00258     return(-1.0);
00259   }
00260 
00261   virtual double Condest() const
00262   {
00263     return(-1.0);
00264   }
00265 
00266   std::ostream& Print(std::ostream& os) const;
00267 
00269   virtual int NumInitialize() const
00270   {
00271     return(NumInitialize_);
00272   }
00273 
00275   virtual int NumCompute() const
00276   {
00277     return(NumCompute_);
00278   }
00279 
00281   virtual int NumApplyInverse() const
00282   {
00283     return(NumApplyInverse_);
00284   }
00285 
00287   virtual double InitializeTime() const
00288   {
00289     return(InitializeTime_);
00290   }
00291 
00293   virtual double ComputeTime() const
00294   {
00295     return(ComputeTime_);
00296   }
00297 
00299   virtual double ApplyInverseTime() const
00300   {
00301     return(ApplyInverseTime_);
00302   }
00303 
00305   virtual double InitializeFlops() const
00306   {
00307 #ifdef IFPACK_FLOPCOUNTERS
00308     if (Containers_.size() == 0)
00309       return(0.0);
00310 
00311     // the total number of flops is computed each time InitializeFlops() is
00312     // called. This is becase I also have to add the contribution from each
00313     // container.
00314     double total = InitializeFlops_;
00315     for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
00316       total += Containers_[i]->InitializeFlops();
00317     return(total);
00318 #else
00319     return(0.0);
00320 #endif
00321   }
00322 
00323   virtual double ComputeFlops() const
00324   {
00325 #ifdef IFPACK_FLOPCOUNTERS
00326     if (Containers_.size() == 0)
00327       return(0.0);
00328     
00329     double total = ComputeFlops_;
00330     for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
00331       total += Containers_[i]->ComputeFlops();
00332     return(total);
00333 #else
00334     return(0.0);
00335 #endif
00336   }
00337 
00338   virtual double ApplyInverseFlops() const
00339   {
00340 #ifdef IFPACK_FLOPCOUNTERS
00341     if (Containers_.size() == 0)
00342       return(0.0);
00343 
00344     double total = ApplyInverseFlops_;
00345     for (unsigned int i = 0 ; i < Containers_.size() ; ++i) {
00346       total += Containers_[i]->ApplyInverseFlops();
00347     }
00348     return(total);
00349 #else
00350     return(0.0);
00351 #endif
00352   }
00353 
00354 private:
00355 
00357   Ifpack_BlockRelaxation(const Ifpack_BlockRelaxation& rhs);
00358 
00360   Ifpack_BlockRelaxation & operator=(const Ifpack_BlockRelaxation& rhs)
00361   {
00362     return(*this);
00363   }
00364 
00365   virtual int ApplyInverseJacobi(const Epetra_MultiVector& X, 
00366                                  Epetra_MultiVector& Y) const;
00367 
00368   virtual int DoJacobi(const Epetra_MultiVector& X, 
00369                                   Epetra_MultiVector& Y) const;
00370 
00371   virtual int ApplyInverseGS(const Epetra_MultiVector& X, 
00372                              Epetra_MultiVector& Y) const;
00373 
00374   virtual int DoGaussSeidel(Epetra_MultiVector& X, 
00375                             Epetra_MultiVector& Y) const;
00376 
00377   virtual int ApplyInverseSGS(const Epetra_MultiVector& X, 
00378                               Epetra_MultiVector& Y) const;
00379 
00380   virtual int DoSGS(const Epetra_MultiVector& X,
00381                     Epetra_MultiVector& Xtmp,
00382                     Epetra_MultiVector& Y) const;
00383 
00384   int ExtractSubmatrices();
00385 
00386   // @{ Initializations, timing and flops
00387 
00389   bool IsInitialized_;
00391   bool IsComputed_;
00393   int NumInitialize_;
00395   int NumCompute_;
00397   mutable int NumApplyInverse_;
00399   double InitializeTime_;
00401   double ComputeTime_;
00403   mutable double ApplyInverseTime_;
00405   double InitializeFlops_;
00407   double ComputeFlops_;
00409   mutable double ApplyInverseFlops_;
00410   // @}
00411 
00412   // @{ Settings
00414   int NumSweeps_;
00416   double DampingFactor_;
00418   int NumLocalBlocks_;
00420   Teuchos::ParameterList List_;
00421   // @}
00422 
00423   // @{ Other data
00426   Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_;
00427   mutable std::vector<Teuchos::RefCountPtr<T> > Containers_;
00429   Teuchos::RefCountPtr<Ifpack_Partitioner> Partitioner_;
00430   string PartitionerType_;
00431   int PrecType_;
00433   string Label_;
00435   bool ZeroStartingSolution_;
00436   Teuchos::RefCountPtr<Ifpack_Graph> Graph_;
00438   Teuchos::RefCountPtr<Epetra_Vector> W_;
00439   // Level of overlap among blocks (for Jacobi only).
00440   int OverlapLevel_;
00441   mutable Epetra_Time Time_;
00442   bool IsParallel_;
00443   Teuchos::RefCountPtr<Epetra_Import> Importer_;
00444   // @}
00445   
00446 }; // class Ifpack_BlockRelaxation
00447 
00448 //==============================================================================
00449 template<typename T>
00450 Ifpack_BlockRelaxation<T>::
00451 Ifpack_BlockRelaxation(const Epetra_RowMatrix* Matrix_in) :
00452   IsInitialized_(false),
00453   IsComputed_(false),
00454   NumInitialize_(0),
00455   NumCompute_(0),
00456   NumApplyInverse_(0),
00457   InitializeTime_(0.0),
00458   ComputeTime_(0.0),
00459   ApplyInverseTime_(0.0),
00460   InitializeFlops_(0.0),
00461   ComputeFlops_(0.0),
00462   ApplyInverseFlops_(0.0),
00463   NumSweeps_(1),
00464   DampingFactor_(1.0),
00465   NumLocalBlocks_(1),
00466   Matrix_(Teuchos::rcp(Matrix_in,false)),
00467   PartitionerType_("greedy"),
00468   PrecType_(IFPACK_JACOBI),
00469   ZeroStartingSolution_(true),
00470   OverlapLevel_(0),
00471   Time_(Comm()),
00472   IsParallel_(false)
00473 {
00474   if (Matrix_in->Comm().NumProc() != 1)
00475     IsParallel_ = true;
00476 }
00477 
00478 //==============================================================================
00479 template<typename T>
00480 Ifpack_BlockRelaxation<T>::~Ifpack_BlockRelaxation()
00481 {
00482 }
00483 
00484 //==============================================================================
00485 template<typename T>
00486 const char* Ifpack_BlockRelaxation<T>::Label() const
00487 {
00488   return(Label_.c_str());
00489 }
00490 
00491 //==============================================================================
00492 template<typename T>
00493 int Ifpack_BlockRelaxation<T>::
00494 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00495 {
00496   int ierr = Matrix().Apply(X,Y);
00497   IFPACK_RETURN(ierr);
00498 }
00499 
00500 //==============================================================================
00501 template<typename T>
00502 const Epetra_Comm& Ifpack_BlockRelaxation<T>::
00503 Comm() const
00504 {
00505   return(Matrix().Comm());
00506 }
00507 
00508 //==============================================================================
00509 template<typename T>
00510 const Epetra_Map& Ifpack_BlockRelaxation<T>::
00511 OperatorDomainMap() const
00512 {
00513   return(Matrix().OperatorDomainMap());
00514 }
00515 
00516 //==============================================================================
00517 template<typename T>
00518 const Epetra_Map& Ifpack_BlockRelaxation<T>::
00519 OperatorRangeMap() const
00520 {
00521   return(Matrix().OperatorRangeMap());
00522 }
00523 
00524 //==============================================================================
00525 template<typename T>
00526 int Ifpack_BlockRelaxation<T>::ExtractSubmatrices()
00527 {
00528 
00529   if (Partitioner_ == Teuchos::null)
00530     IFPACK_CHK_ERR(-3);
00531 
00532   NumLocalBlocks_ = Partitioner_->NumLocalParts();
00533 
00534   Containers_.resize(NumLocalBlocks());
00535 
00536   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00537 
00538     int rows = Partitioner_->NumRowsInPart(i);
00539     Containers_[i] = Teuchos::rcp( new T(rows) );
00540     
00541     //Ifpack_DenseContainer* DC = 0;
00542     //DC = dynamic_cast<Ifpack_DenseContainer*>(Containers_[i]);
00543 
00544     if (Containers_[i] == Teuchos::null)
00545       IFPACK_CHK_ERR(-5);
00546     
00547     IFPACK_CHK_ERR(Containers_[i]->SetParameters(List_));
00548     IFPACK_CHK_ERR(Containers_[i]->Initialize());
00549     // flops in Initialize() will be computed on-the-fly in method InitializeFlops().
00550 
00551     // set "global" ID of each partitioner row
00552     for (int j = 0 ; j < rows ; ++j) {
00553       int LRID = (*Partitioner_)(i,j);
00554       Containers_[i]->ID(j) = LRID;
00555     }
00556 
00557     IFPACK_CHK_ERR(Containers_[i]->Compute(*Matrix_));
00558     // flops in Compute() will be computed on-the-fly in method ComputeFlops().
00559 
00560   }
00561 
00562   return(0);
00563 }
00564 
00565 //==============================================================================
00566 template<typename T>
00567 int Ifpack_BlockRelaxation<T>::Compute()
00568 {
00569 
00570   if (!IsInitialized())
00571     IFPACK_CHK_ERR(Initialize());
00572 
00573   Time_.ResetStartTime();
00574 
00575   IsComputed_ = false;
00576 
00577   if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
00578     IFPACK_CHK_ERR(-2); // only square matrices
00579 
00580   IFPACK_CHK_ERR(ExtractSubmatrices());
00581   
00582   if (IsParallel_ && PrecType_ != IFPACK_JACOBI) {
00583     // not needed by Jacobi (done by matvec)
00584     Importer_ = Teuchos::rcp( new Epetra_Import(Matrix().RowMatrixColMap(),
00585                                                 Matrix().RowMatrixRowMap()) );
00586 
00587     if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
00588   }
00589   IsComputed_ = true;
00590   ComputeTime_ += Time_.ElapsedTime();
00591   ++NumCompute_;
00592 
00593   return(0);
00594 
00595 }
00596 
00597 //==============================================================================
00598 template<typename T>
00599 int Ifpack_BlockRelaxation<T>::
00600 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00601 {
00602   if (!IsComputed())
00603     IFPACK_CHK_ERR(-3);
00604 
00605   if (X.NumVectors() != Y.NumVectors())
00606     IFPACK_CHK_ERR(-2);
00607 
00608   Time_.ResetStartTime();
00609 
00610   // AztecOO gives X and Y pointing to the same memory location,
00611   // need to create an auxiliary vector, Xcopy
00612   Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
00613   if (X.Pointers()[0] == Y.Pointers()[0])
00614     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00615   else
00616     Xcopy = Teuchos::rcp( &X, false );
00617 
00618   switch (PrecType_) {
00619   case IFPACK_JACOBI:
00620     IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
00621     break;
00622   case IFPACK_GS:
00623     IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
00624     break;
00625   case IFPACK_SGS:
00626     IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
00627     break;
00628   }
00629 
00630   ApplyInverseTime_ += Time_.ElapsedTime();
00631   ++NumApplyInverse_;
00632 
00633   return(0);
00634 }
00635 
00636 //==============================================================================
00637 // This method in general will not work with AztecOO if used
00638 // outside Ifpack_AdditiveSchwarz and OverlapLevel_ != 0
00639 //
00640 template<typename T>
00641 int Ifpack_BlockRelaxation<T>::
00642 ApplyInverseJacobi(const Epetra_MultiVector& X, 
00643                    Epetra_MultiVector& Y) const
00644 {
00645 
00646   if (ZeroStartingSolution_)
00647     Y.PutScalar(0.0);
00648 
00649   // do not compute the residual in this case
00650   if (NumSweeps_ == 1 && ZeroStartingSolution_) {
00651     int ierr = DoJacobi(X,Y);
00652     IFPACK_RETURN(ierr);
00653   }
00654 
00655   Epetra_MultiVector AX(Y);
00656 
00657   for (int j = 0; j < NumSweeps_ ; j++) {
00658     IFPACK_CHK_ERR(Apply(Y,AX));
00659     ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalNonzeros64();
00660     IFPACK_CHK_ERR(AX.Update(1.0,X,-1.0));
00661     ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalRows64();
00662     IFPACK_CHK_ERR(DoJacobi(AX,Y));
00663     // flops counted in DoJacobi()
00664   }
00665 
00666 
00667   return(0);
00668 }
00669 
00670 //==============================================================================
00671 template<typename T>
00672 int Ifpack_BlockRelaxation<T>::
00673 DoJacobi(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00674 {
00675   int NumVectors = X.NumVectors();
00676 
00677   if (OverlapLevel_ == 0) {
00678 
00679     for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00680      
00681       // may happen that a partition is empty
00682       if (Containers_[i]->NumRows() == 0) 
00683         continue;
00684 
00685       int LID;
00686 
00687       // extract RHS from X
00688       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00689         LID = Containers_[i]->ID(j);
00690         for (int k = 0 ; k < NumVectors ; ++k) {
00691           Containers_[i]->RHS(j,k) = X[k][LID];
00692         }
00693       }
00694 
00695       // apply the inverse of each block. NOTE: flops occurred
00696       // in ApplyInverse() of each block are summed up in method
00697       // ApplyInverseFlops().
00698       IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00699 
00700       // copy back into solution vector Y
00701       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00702         LID = Containers_[i]->ID(j);
00703         for (int k = 0 ; k < NumVectors ; ++k) {
00704           Y[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00705         }
00706       }
00707 
00708     }
00709     // NOTE: flops for ApplyInverse() of each block are summed up
00710     // in method ApplyInverseFlops()
00711 #ifdef IFPACK_FLOPCOUNTERS
00712     ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00713 #endif
00714 
00715   }
00716   else {
00717 
00718     for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00719 
00720       // may happen that a partition is empty
00721       if (Containers_[i]->NumRows() == 0) 
00722         continue;
00723 
00724       int LID;
00725 
00726       // extract RHS from X
00727       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00728         LID = Containers_[i]->ID(j);
00729         for (int k = 0 ; k < NumVectors ; ++k) {
00730           Containers_[i]->RHS(j,k) = (*W_)[LID] * X[k][LID];
00731         }
00732       }
00733 
00734       // apply the inverse of each block
00735       IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00736 
00737       // copy back into solution vector Y
00738       for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00739         LID = Containers_[i]->ID(j);
00740         for (int k = 0 ; k < NumVectors ; ++k) {
00741           Y[k][LID] += DampingFactor_ * (*W_)[LID] * Containers_[i]->LHS(j,k);
00742         }
00743       }
00744 
00745     }
00746     // NOTE: flops for ApplyInverse() of each block are summed up
00747     // in method ApplyInverseFlops()
00748     // NOTE: do not count for simplicity the flops due to overlapping rows
00749 #ifdef IFPACK_FLOPCOUNTERS
00750     ApplyInverseFlops_ += NumVectors * 4 * Matrix_->NumGlobalRows();
00751 #endif
00752   }
00753 
00754   return(0);
00755 }
00756 
00757 //==============================================================================
00758 template<typename T>
00759 int Ifpack_BlockRelaxation<T>::
00760 ApplyInverseGS(const Epetra_MultiVector& X, 
00761                Epetra_MultiVector& Y) const
00762 {
00763 
00764   if (ZeroStartingSolution_)
00765     Y.PutScalar(0.0);
00766 
00767   Epetra_MultiVector Xcopy(X);
00768   for (int j = 0; j < NumSweeps_ ; j++) {
00769     IFPACK_CHK_ERR(DoGaussSeidel(Xcopy,Y));
00770     if (j != NumSweeps_ - 1)
00771       Xcopy = X;
00772   }
00773 
00774   return(0);
00775 
00776 }
00777 
00778 //==============================================================================
00779 template<typename T>
00780 int Ifpack_BlockRelaxation<T>::
00781 DoGaussSeidel(Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00782 {
00783 
00784   // cycle over all local subdomains
00785 
00786   int Length = Matrix().MaxNumEntries();
00787   std::vector<int> Indices(Length);
00788   std::vector<double> Values(Length);
00789 
00790   int NumMyRows = Matrix().NumMyRows();
00791   int NumVectors = X.NumVectors();
00792 
00793   // an additonal vector is needed by parallel computations
00794   // (note that applications through Ifpack_AdditiveSchwarz
00795   // are always seen are serial)
00796   Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
00797   if (IsParallel_)
00798     Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
00799   else
00800     Y2 = Teuchos::rcp( &Y, false );
00801 
00802   double** y_ptr;
00803   double** y2_ptr;
00804   Y.ExtractView(&y_ptr);
00805   Y2->ExtractView(&y2_ptr);
00806 
00807   // data exchange is here, once per sweep
00808   if (IsParallel_)
00809     IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00810 
00811   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00812 
00813     // may happen that a partition is empty
00814     if (Containers_[i]->NumRows() == 0) 
00815       continue;
00816 
00817     int LID;
00818 
00819     // update from previous block
00820 
00821     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00822       LID = Containers_[i]->ID(j);
00823 
00824       int NumEntries;
00825       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00826                                                &Values[0], &Indices[0]));
00827 
00828       for (int k = 0 ; k < NumEntries ; ++k) {
00829         int col = Indices[k];
00830 
00831           for (int kk = 0 ; kk < NumVectors ; ++kk) {
00832             X[kk][LID] -= Values[k] * y2_ptr[kk][col];
00833           }
00834       }
00835     }
00836 
00837     // solve with this block
00838 
00839     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00840       LID = Containers_[i]->ID(j);
00841       for (int k = 0 ; k < NumVectors ; ++k) {
00842         Containers_[i]->RHS(j,k) = X[k][LID];
00843       }
00844     }
00845 
00846     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00847 #ifdef IFPACK_FLOPCOUNTERS
00848     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00849 #endif
00850 
00851     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00852       LID = Containers_[i]->ID(j);
00853       for (int k = 0 ; k < NumVectors ; ++k) {
00854         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00855       }
00856     }
00857 
00858   }
00859 
00860   // operations for all getrow()'s
00861   // NOTE: flops for ApplyInverse() of each block are summed up
00862   // in method ApplyInverseFlops()
00863 #ifdef IFPACK_FLOPCOUNTERS
00864   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00865   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00866 #endif
00867 
00868   // Attention: this is delicate... Not all combinations
00869   // of Y2 and Y will always work (tough for ML it should be ok)
00870   if (IsParallel_)
00871     for (int m = 0 ; m < NumVectors ; ++m) 
00872       for (int i = 0 ; i < NumMyRows ; ++i)
00873         y_ptr[m][i] = y2_ptr[m][i];
00874 
00875   return(0);
00876 }
00877 
00878 //==============================================================================
00879 template<typename T>
00880 int Ifpack_BlockRelaxation<T>::
00881 ApplyInverseSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00882 {
00883 
00884   if (ZeroStartingSolution_)
00885     Y.PutScalar(0.0);
00886 
00887   Epetra_MultiVector Xcopy(X);
00888   for (int j = 0; j < NumSweeps_ ; j++) {
00889     IFPACK_CHK_ERR(DoSGS(X,Xcopy,Y));
00890     if (j != NumSweeps_ - 1)
00891       Xcopy = X;
00892   }
00893   return(0);
00894 }
00895 
00896 //==============================================================================
00897 template<typename T>
00898 int Ifpack_BlockRelaxation<T>::
00899 DoSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Xcopy, 
00900       Epetra_MultiVector& Y) const
00901 {
00902 
00903   int NumMyRows = Matrix().NumMyRows();
00904   int NumVectors = X.NumVectors();
00905 
00906   int Length = Matrix().MaxNumEntries();
00907   std::vector<int> Indices;
00908   std::vector<double> Values;
00909   Indices.resize(Length);
00910   Values.resize(Length);
00911 
00912   // an additonal vector is needed by parallel computations
00913   // (note that applications through Ifpack_AdditiveSchwarz
00914   // are always seen are serial)
00915   Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
00916   if (IsParallel_)
00917     Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
00918   else
00919     Y2 = Teuchos::rcp( &Y, false );
00920 
00921   double** y_ptr;
00922   double** y2_ptr;
00923   Y.ExtractView(&y_ptr);
00924   Y2->ExtractView(&y2_ptr);
00925 
00926   // data exchange is here, once per sweep
00927   if (IsParallel_)
00928     IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00929 
00930   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00931 
00932     // may happen that a partition is empty
00933     if (Containers_[i]->NumRows() == 0) 
00934       continue;
00935 
00936     int LID;
00937 
00938     // update from previous block
00939 
00940     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00941       LID = Containers_[i]->ID(j);
00942 
00943       int NumEntries;
00944       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
00945                                                &Values[0], &Indices[0]));
00946 
00947       for (int k = 0 ; k < NumEntries ; ++k) {
00948         int col = Indices[k];
00949 
00950         for (int kk = 0 ; kk < NumVectors ; ++kk) {
00951           Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
00952         }
00953       }
00954     }
00955 
00956     // solve with this block
00957 
00958     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00959       LID = Containers_[i]->ID(j);
00960       for (int k = 0 ; k < NumVectors ; ++k) {
00961         Containers_[i]->RHS(j,k) = Xcopy[k][LID];
00962       }
00963     }
00964 
00965     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00966 #ifdef IFPACK_FLOPCOUNTERS
00967     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
00968 #endif
00969 
00970     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00971       LID = Containers_[i]->ID(j);
00972       for (int k = 0 ; k < NumVectors ; ++k) {
00973         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
00974       }
00975     }
00976   }
00977 
00978 #ifdef IFPACK_FLOPCOUNTERS
00979   // operations for all getrow()'s
00980   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00981   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00982 #endif
00983 
00984   Xcopy = X;
00985 
00986   for (int i = NumLocalBlocks() - 1; i >=0 ; --i) {
00987 
00988     if (Containers_[i]->NumRows() == 0) 
00989       continue;
00990 
00991     int LID;
00992 
00993     // update from previous block
00994 
00995     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
00996       LID = Containers_[i]->ID(j);
00997 
00998       int NumEntries;
00999       IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
01000                                                &Values[0], &Indices[0]));
01001 
01002       for (int k = 0 ; k < NumEntries ; ++k) {
01003         int col = Indices[k];
01004 
01005           for (int kk = 0 ; kk < NumVectors ; ++kk) {
01006             Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
01007           }
01008       }
01009     }
01010 
01011     // solve with this block
01012 
01013     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
01014       LID = Containers_[i]->ID(j);
01015       for (int k = 0 ; k < NumVectors ; ++k) {
01016         Containers_[i]->RHS(j,k) = Xcopy[k][LID];
01017       }
01018     }
01019 
01020     IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
01021 #ifdef IFPACK_FLOPCOUNTERS
01022     ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
01023 #endif
01024 
01025     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
01026       LID = Containers_[i]->ID(j);
01027       for (int k = 0 ; k < NumVectors ; ++k) {
01028         y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
01029       }
01030     }
01031   }
01032 
01033 #ifdef IFPACK_FLOPCOUNTERS
01034   // operations for all getrow()'s
01035   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
01036   ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
01037 #endif
01038 
01039   // Attention: this is delicate... Not all combinations
01040   // of Y2 and Y will always work (tough for ML it should be ok)
01041   if (IsParallel_)
01042     for (int m = 0 ; m < NumVectors ; ++m) 
01043       for (int i = 0 ; i < NumMyRows ; ++i)
01044         y_ptr[m][i] = y2_ptr[m][i];
01045 
01046   return(0);
01047 }
01048 
01049 //==============================================================================
01050 template<typename T>
01051 ostream& Ifpack_BlockRelaxation<T>::Print(ostream & os) const
01052 {
01053 
01054   string PT;
01055   if (PrecType_ == IFPACK_JACOBI)
01056     PT = "Jacobi";
01057   else if (PrecType_ == IFPACK_GS)
01058     PT = "Gauss-Seidel";
01059   else if (PrecType_ == IFPACK_SGS)
01060     PT = "symmetric Gauss-Seidel";
01061 
01062   if (!Comm().MyPID()) {
01063     os << endl;
01064     os << "================================================================================" << endl;
01065     os << "Ifpack_BlockRelaxation, " << PT << endl;
01066     os << "Sweeps = " << NumSweeps_ << endl;
01067     os << "Damping factor = " << DampingFactor_;
01068     if (ZeroStartingSolution_) 
01069       os << ", using zero starting solution" << endl;
01070     else
01071       os << ", using input starting solution" << endl;
01072     os << "Number of local blocks = " << Partitioner_->NumLocalParts() << endl;
01073     //os << "Condition number estimate = " << Condest_ << endl;
01074     os << "Global number of rows            = " << Matrix_->NumGlobalRows64() << endl;
01075     os << endl;
01076     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
01077     os << "-----           -------   --------------       ------------     --------" << endl;
01078     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
01079        << "  " << std::setw(15) << InitializeTime() 
01080        << "  " << std::setw(15) << 1.0e-6 * InitializeFlops();
01081     if (InitializeTime() != 0.0)
01082       os << "  " << std::setw(15) << 1.0e-6 * InitializeFlops() / InitializeTime() << endl;
01083     else
01084       os << "  " << std::setw(15) << 0.0 << endl;
01085     os << "Compute()       "   << std::setw(5) << NumCompute() 
01086        << "  " << std::setw(15) << ComputeTime()
01087        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
01088     if (ComputeTime() != 0.0) 
01089       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
01090     else
01091       os << "  " << std::setw(15) << 0.0 << endl;
01092     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
01093        << "  " << std::setw(15) << ApplyInverseTime()
01094        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
01095     if (ApplyInverseTime() != 0.0) 
01096       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
01097     else
01098       os << "  " << std::setw(15) << 0.0 << endl;
01099     os << "================================================================================" << endl;
01100     os << endl;
01101   }
01102 
01103   return(os);
01104 }
01105 
01106 //==============================================================================
01107 template<typename T>
01108 int Ifpack_BlockRelaxation<T>::SetParameters(Teuchos::ParameterList& List)
01109 {
01110 
01111   string PT;
01112   if (PrecType_ == IFPACK_JACOBI)
01113     PT = "Jacobi";
01114   else if (PrecType_ == IFPACK_GS)
01115     PT = "Gauss-Seidel";
01116   else if (PrecType_ == IFPACK_SGS)
01117     PT = "symmetric Gauss-Seidel";
01118 
01119   PT = List.get("relaxation: type", PT);
01120 
01121   if (PT == "Jacobi") {
01122     PrecType_ = IFPACK_JACOBI;
01123   }
01124   else if (PT == "Gauss-Seidel") {
01125     PrecType_ = IFPACK_GS;
01126   }
01127   else if (PT == "symmetric Gauss-Seidel") {
01128     PrecType_ = IFPACK_SGS;
01129   } else {
01130     cerr << "Option `relaxation: type' has an incorrect value ("
01131       << PT << ")'" << endl;
01132     cerr << "(file " << __FILE__ << ", line " << __LINE__ << ")" << endl;
01133     exit(EXIT_FAILURE);
01134   }
01135 
01136   NumSweeps_            = List.get("relaxation: sweeps", NumSweeps_);
01137   DampingFactor_        = List.get("relaxation: damping factor", 
01138                                    DampingFactor_);
01139   ZeroStartingSolution_ = List.get("relaxation: zero starting solution", 
01140                                    ZeroStartingSolution_);
01141   PartitionerType_      = List.get("partitioner: type", 
01142                                    PartitionerType_);
01143   NumLocalBlocks_       = List.get("partitioner: local parts", 
01144                                    NumLocalBlocks_);
01145   // only Jacobi can work with overlap among local domains,
01146   OverlapLevel_         = List.get("partitioner: overlap", 
01147                                    OverlapLevel_);
01148 
01149   // check parameters
01150   if (PrecType_ != IFPACK_JACOBI)
01151     OverlapLevel_ = 0;
01152   if (NumLocalBlocks_ < 0)
01153     NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
01154   // other checks are performed in Partitioner_
01155   
01156   // copy the list as each subblock's constructor will
01157   // require it later
01158   List_ = List;
01159 
01160   // set the label
01161   string PT2;
01162   if (PrecType_ == IFPACK_JACOBI)
01163     PT2 = "BJ";
01164   else if (PrecType_ == IFPACK_GS)
01165     PT2 = "BGS";
01166   else if (PrecType_ == IFPACK_SGS)
01167     PT2 = "BSGS";
01168   Label_ = "IFPACK (" + PT2 + ", sweeps=" 
01169     + Ifpack_toString(NumSweeps_) + ", damping="
01170     + Ifpack_toString(DampingFactor_) + ", blocks="
01171     + Ifpack_toString(NumLocalBlocks()) + ")";
01172 
01173   return(0);
01174 }
01175 
01176 //==============================================================================
01177 template<typename T>
01178 int Ifpack_BlockRelaxation<T>::Initialize()
01179 {
01180   IsInitialized_ = false;
01181   Time_.ResetStartTime();
01182 
01183   Graph_ = Teuchos::rcp( new Ifpack_Graph_Epetra_RowMatrix(Teuchos::rcp(&Matrix(),false)) );
01184   if (Graph_ == Teuchos::null) IFPACK_CHK_ERR(-5);
01185 
01186   if (PartitionerType_ == "linear")
01187     Partitioner_ = Teuchos::rcp( new Ifpack_LinearPartitioner(&*Graph_) );
01188   else if (PartitionerType_ == "greedy")
01189     Partitioner_ = Teuchos::rcp( new Ifpack_GreedyPartitioner(&*Graph_) );
01190   else if (PartitionerType_ == "metis")
01191     Partitioner_ = Teuchos::rcp( new Ifpack_METISPartitioner(&*Graph_) );
01192   else if (PartitionerType_ == "equation")
01193     Partitioner_ = Teuchos::rcp( new Ifpack_EquationPartitioner(&*Graph_) );
01194   else if (PartitionerType_ == "user")
01195     Partitioner_ = Teuchos::rcp( new Ifpack_UserPartitioner(&*Graph_) );
01196   else if (PartitionerType_ == "line")
01197     Partitioner_ = Teuchos::rcp( new Ifpack_LinePartitioner(&*Graph_) );
01198   else
01199     IFPACK_CHK_ERR(-2);
01200 
01201   if (Partitioner_ == Teuchos::null) IFPACK_CHK_ERR(-5);
01202 
01203   // need to partition the graph of A
01204   IFPACK_CHK_ERR(Partitioner_->SetParameters(List_));
01205   IFPACK_CHK_ERR(Partitioner_->Compute());
01206 
01207   // get actual number of partitions
01208   NumLocalBlocks_ = Partitioner_->NumLocalParts();
01209 
01210   // weight of each vertex
01211   W_ = Teuchos::rcp( new Epetra_Vector(Matrix().RowMatrixRowMap()) );
01212   W_->PutScalar(0.0);
01213 
01214   for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
01215 
01216     for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
01217       int LID = (*Partitioner_)(i,j);
01218       (*W_)[LID]++;
01219     }
01220   }
01221   W_->Reciprocal(*W_);
01222 
01223   // Update Label_ if line smoothing
01224   if (PartitionerType_ == "line") {
01225     // set the label
01226     string PT2;
01227     if (PrecType_ == IFPACK_JACOBI)
01228       PT2 = "BJ";
01229     else if (PrecType_ == IFPACK_GS)
01230       PT2 = "BGS";
01231     else if (PrecType_ == IFPACK_SGS)
01232       PT2 = "BSGS";
01233     Label_ = "IFPACK (" + PT2 + ", auto-line, sweeps=" 
01234       + Ifpack_toString(NumSweeps_) + ", damping="
01235       + Ifpack_toString(DampingFactor_) + ", blocks="
01236       + Ifpack_toString(NumLocalBlocks()) + ")";
01237   }
01238 
01239   InitializeTime_ += Time_.ElapsedTime();
01240   IsInitialized_ = true;
01241   ++NumInitialize_;
01242 
01243   return(0);
01244 }
01245 
01246 //==============================================================================
01247 #endif // IFPACK_BLOCKPRECONDITIONER_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends