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);
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
00267
00268
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
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
00358 int NumSweeps_;
00360 double DampingFactor_;
00362 int NumLocalBlocks_;
00364 Teuchos::ParameterList List_;
00365
00366
00367
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
00384 int OverlapLevel_;
00385 mutable Epetra_Time Time_;
00386 bool IsParallel_;
00387 Teuchos::RefCountPtr<Epetra_Import> Importer_;
00388
00389
00390 };
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
00485
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
00493
00494
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
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);
00522
00523 IFPACK_CHK_ERR(ExtractSubmatrices());
00524
00525 if (IsParallel_ && PrecType_ != IFPACK_JACOBI) {
00526
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
00554
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
00581
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
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
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
00624 if (Containers_[i]->NumRows() == 0)
00625 continue;
00626
00627 int LID;
00628
00629
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
00638
00639
00640 IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00641
00642
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
00652
00653 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00654
00655 }
00656 else {
00657
00658 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00659
00660
00661 if (Containers_[i]->NumRows() == 0)
00662 continue;
00663
00664 int LID;
00665
00666
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
00675 IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
00676
00677
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
00687
00688
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
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
00732
00733
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
00746 if (IsParallel_)
00747 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00748
00749 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00750
00751
00752 if (Containers_[i]->NumRows() == 0)
00753 continue;
00754
00755 int LID;
00756
00757
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
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
00797
00798
00799 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00800 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00801
00802
00803
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
00847
00848
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
00861 if (IsParallel_)
00862 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
00863
00864 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
00865
00866
00867 if (Containers_[i]->NumRows() == 0)
00868 continue;
00869
00870 int LID;
00871
00872
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
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
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
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
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
00962 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
00963 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
00964
00965
00966
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
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
01072 OverlapLevel_ = List.get("partitioner: overlap",
01073 OverlapLevel_);
01074
01075
01076 if (PrecType_ != IFPACK_JACOBI)
01077 OverlapLevel_ = 0;
01078 if (NumLocalBlocks_ < 0)
01079 NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
01080
01081
01082
01083
01084 List_ = List;
01085
01086
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
01128 IFPACK_CHK_ERR(Partitioner_->SetParameters(List_));
01129 IFPACK_CHK_ERR(Partitioner_->Compute());
01130
01131
01132 NumLocalBlocks_ = Partitioner_->NumLocalParts();
01133
01134
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