IFPACK Development
Loading...
Searching...
No Matches
Ifpack_BlockRelaxation.h
1/*
2//@HEADER
3// ***********************************************************************
4//
5// Ifpack: Object-Oriented Algebraic Preconditioner Package
6// Copyright (2002) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44#ifndef IFPACK_BLOCKPRECONDITIONER_H
45#define IFPACK_BLOCKPRECONDITIONER_H
46
47#include "Ifpack_ConfigDefs.h"
48#include "Ifpack_Preconditioner.h"
49#include "Ifpack_Partitioner.h"
50#include "Ifpack_LinePartitioner.h"
51#include "Ifpack_LinearPartitioner.h"
52#include "Ifpack_GreedyPartitioner.h"
53#include "Ifpack_METISPartitioner.h"
54#include "Ifpack_EquationPartitioner.h"
55#include "Ifpack_UserPartitioner.h"
56#include "Ifpack_Graph_Epetra_RowMatrix.h"
57#include "Ifpack_DenseContainer.h"
58#include "Ifpack_Utils.h"
59#include "Teuchos_ParameterList.hpp"
60#include "Teuchos_RefCountPtr.hpp"
61
62#include "Epetra_Map.h"
63#include "Epetra_RowMatrix.h"
64#include "Epetra_MultiVector.h"
65#include "Epetra_Vector.h"
66#include "Epetra_Time.h"
67#include "Epetra_Import.h"
68
69static const int IFPACK_JACOBI = 0;
70static const int IFPACK_GS = 1;
71static const int IFPACK_SGS = 2;
72
73
75
138template<typename T>
140
141public:
142
144
151
152 virtual ~Ifpack_BlockRelaxation();
153
155
157
159
167 virtual int Apply(const Epetra_MultiVector& X,
168 Epetra_MultiVector& Y) const;
169
171
180 virtual int ApplyInverse(const Epetra_MultiVector& X,
181 Epetra_MultiVector& Y) const;
182
184 virtual double NormInf() const
185 {
186 return(-1.0);
187 }
189
191
192 virtual int SetUseTranspose(bool UseTranspose_in)
193 {
194 if (UseTranspose_in)
195 IFPACK_CHK_ERR(-98); // FIXME: can I work with the transpose?
196 return(0);
197 }
198
199 virtual const char* Label() const;
200
202 virtual bool UseTranspose() const
203 {
204 return(false);
205 }
206
208 virtual bool HasNormInf() const
209 {
210 return(false);
211 }
212
214 virtual const Epetra_Comm & Comm() const;
215
217 virtual const Epetra_Map & OperatorDomainMap() const;
218
220 virtual const Epetra_Map & OperatorRangeMap() const;
222
224 int NumLocalBlocks() const
225 {
226 return(NumLocalBlocks_);
227 }
228
230 virtual bool IsInitialized() const
231 {
232 return(IsInitialized_);
233 }
234
236 virtual bool IsComputed() const
237 {
238 return(IsComputed_);
239 }
240
242 virtual int SetParameters(Teuchos::ParameterList& List);
243
245 virtual int Initialize();
246
248 virtual int Compute();
249
250 virtual const Epetra_RowMatrix& Matrix() const
251 {
252 return(*Matrix_);
253 }
254
255 virtual double Condest(const Ifpack_CondestType /* CT */ = Ifpack_Cheap,
256 const int /* MaxIters */ = 1550,
257 const double /* Tol */ = 1e-9,
258 Epetra_RowMatrix* /* Matrix_in */ = 0)
259 {
260 return(-1.0);
261 }
262
263 virtual double Condest() const
264 {
265 return(-1.0);
266 }
267
268 std::ostream& Print(std::ostream& os) const;
269
271 virtual int NumInitialize() const
272 {
273 return(NumInitialize_);
274 }
275
277 virtual int NumCompute() const
278 {
279 return(NumCompute_);
280 }
281
283 virtual int NumApplyInverse() const
284 {
285 return(NumApplyInverse_);
286 }
287
289 virtual double InitializeTime() const
290 {
291 return(InitializeTime_);
292 }
293
295 virtual double ComputeTime() const
296 {
297 return(ComputeTime_);
298 }
299
301 virtual double ApplyInverseTime() const
302 {
303 return(ApplyInverseTime_);
304 }
305
307 virtual double InitializeFlops() const
308 {
309#ifdef IFPACK_FLOPCOUNTERS
310 if (Containers_.size() == 0)
311 return(0.0);
312
313 // the total number of flops is computed each time InitializeFlops() is
314 // called. This is becase I also have to add the contribution from each
315 // container.
316 double total = InitializeFlops_;
317 for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
318 if(Containers_[i]) total += Containers_[i]->InitializeFlops();
319 return(total);
320#else
321 return(0.0);
322#endif
323 }
324
325 virtual double ComputeFlops() const
326 {
327#ifdef IFPACK_FLOPCOUNTERS
328 if (Containers_.size() == 0)
329 return(0.0);
330
331 double total = ComputeFlops_;
332 for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
333 if(Containers_[i]) total += Containers_[i]->ComputeFlops();
334 return(total);
335#else
336 return(0.0);
337#endif
338 }
339
340 virtual double ApplyInverseFlops() const
341 {
342#ifdef IFPACK_FLOPCOUNTERS
343 if (Containers_.size() == 0)
344 return(0.0);
345
346 double total = ApplyInverseFlops_;
347 for (unsigned int i = 0 ; i < Containers_.size() ; ++i) {
348 if(Containers_[i]) total += Containers_[i]->ApplyInverseFlops();
349 }
350 return(total);
351#else
352 return(0.0);
353#endif
354 }
355
356private:
357
360
362 Ifpack_BlockRelaxation & operator=(const Ifpack_BlockRelaxation& rhs)
363 {
364 return(*this);
365 }
366
367 virtual int ApplyInverseJacobi(const Epetra_MultiVector& X,
368 Epetra_MultiVector& Y) const;
369
370 virtual int DoJacobi(const Epetra_MultiVector& X,
371 Epetra_MultiVector& Y) const;
372
373 virtual int ApplyInverseGS(const Epetra_MultiVector& X,
374 Epetra_MultiVector& Y) const;
375
376 virtual int DoGaussSeidel(Epetra_MultiVector& X,
377 Epetra_MultiVector& Y) const;
378
379 virtual int ApplyInverseSGS(const Epetra_MultiVector& X,
380 Epetra_MultiVector& Y) const;
381
382 virtual int DoSGS(const Epetra_MultiVector& X,
383 Epetra_MultiVector& Xtmp,
384 Epetra_MultiVector& Y) const;
385
386 int ExtractSubmatrices();
387
388 // @{ Initializations, timing and flops
389
391 bool IsInitialized_;
393 bool IsComputed_;
395 int NumInitialize_;
397 int NumCompute_;
399 mutable int NumApplyInverse_;
401 double InitializeTime_;
403 double ComputeTime_;
405 mutable double ApplyInverseTime_;
407 double InitializeFlops_;
409 double ComputeFlops_;
411 mutable double ApplyInverseFlops_;
412 // @}
413
414 // @{ Settings
416 int NumSweeps_;
418 double DampingFactor_;
420 int NumLocalBlocks_;
422 Teuchos::ParameterList List_;
423 // @}
424
425 // @{ Other data
428 Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_;
429 mutable std::vector<Teuchos::RefCountPtr<T> > Containers_;
430 Epetra_Vector Diagonal_ ;
431
433 Teuchos::RefCountPtr<Ifpack_Partitioner> Partitioner_;
434 std::string PartitionerType_;
435 int PrecType_;
437 std::string Label_;
439 bool ZeroStartingSolution_;
440 Teuchos::RefCountPtr<Ifpack_Graph> Graph_;
442 Teuchos::RefCountPtr<Epetra_Vector> W_;
443 // Level of overlap among blocks (for Jacobi only).
444 int OverlapLevel_;
445 mutable Epetra_Time Time_;
446 bool IsParallel_;
447 Teuchos::RefCountPtr<Epetra_Import> Importer_;
448 // @}
449
450}; // class Ifpack_BlockRelaxation
451
452//==============================================================================
453template<typename T>
456 IsInitialized_(false),
457 IsComputed_(false),
458 NumInitialize_(0),
459 NumCompute_(0),
460 NumApplyInverse_(0),
461 InitializeTime_(0.0),
462 ComputeTime_(0.0),
463 ApplyInverseTime_(0.0),
464 InitializeFlops_(0.0),
465 ComputeFlops_(0.0),
466 ApplyInverseFlops_(0.0),
467 NumSweeps_(1),
468 DampingFactor_(1.0),
469 NumLocalBlocks_(1),
470 Matrix_(Teuchos::rcp(Matrix_in,false)),
471 Diagonal_( Matrix_in->Map()),
472 PartitionerType_("greedy"),
473 PrecType_(IFPACK_JACOBI),
474 ZeroStartingSolution_(true),
475 OverlapLevel_(0),
476 Time_(Comm()),
477 IsParallel_(false)
478{
479 if (Matrix_in->Comm().NumProc() != 1)
480 IsParallel_ = true;
481}
482
483//==============================================================================
484template<typename T>
486{
487}
488
489//==============================================================================
490template<typename T>
491const char* Ifpack_BlockRelaxation<T>::Label() const
492{
493 return(Label_.c_str());
494}
495
496//==============================================================================
497template<typename T>
500{
501 int ierr = Matrix().Apply(X,Y);
502 IFPACK_RETURN(ierr);
503}
504
505//==============================================================================
506template<typename T>
508Comm() const
509{
510 return(Matrix().Comm());
511}
512
513//==============================================================================
514template<typename T>
516OperatorDomainMap() const
517{
518 return(Matrix().OperatorDomainMap());
519}
520
521//==============================================================================
522template<typename T>
524OperatorRangeMap() const
525{
526 return(Matrix().OperatorRangeMap());
527}
528
529//==============================================================================
530template<typename T>
532{
533
534 if (Partitioner_ == Teuchos::null)
535 IFPACK_CHK_ERR(-3);
536
537 NumLocalBlocks_ = Partitioner_->NumLocalParts();
538
539 Containers_.resize(NumLocalBlocks());
540
541 Diagonal_ = Epetra_Vector(Matrix_->Map());
542 Matrix_->ExtractDiagonalCopy(Diagonal_);
543
544 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
545
546 int rows = Partitioner_->NumRowsInPart(i);
547 // if rows == 1, then this is a singleton block, and should not be
548 // created. For now, allow creation, and just force the compute step below.
549
550 if( rows != 1 ) {
551 Containers_[i] = Teuchos::rcp( new T(rows) );
552
553 IFPACK_CHK_ERR(Containers_[i]->SetParameters(List_));
554 IFPACK_CHK_ERR(Containers_[i]->Initialize());
555 // flops in Initialize() will be computed on-the-fly in method InitializeFlops().
556
557 // set "global" ID of each partitioner row
558 for (int j = 0 ; j < rows ; ++j) {
559 int LRID = (*Partitioner_)(i,j);
560 Containers_[i]->ID(j) = LRID;
561 }
562
563 IFPACK_CHK_ERR(Containers_[i]->Compute(*Matrix_));
564 }
565 // otherwise leave Containers_[i] as null
566 }
567
568#ifdef SINGLETON_DEBUG
569 int issing = 0;
570
571 for (int i = 0 ; i < NumLocalBlocks() ; ++i)
572 issing += (int) ( Partitioner_->NumRowsInPart(i) == 1);
573 printf( " %d of %d containers are singleton \n",issing,NumLocalBlocks());
574#endif
575
576 return(0);
577}
578
579//==============================================================================
580template<typename T>
582{
583
584 if (!IsInitialized())
585 IFPACK_CHK_ERR(Initialize());
586
587 Time_.ResetStartTime();
588
589 IsComputed_ = false;
590
591 if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
592 IFPACK_CHK_ERR(-2); // only square matrices
593
594 IFPACK_CHK_ERR(ExtractSubmatrices());
595
596 if (IsParallel_ && PrecType_ != IFPACK_JACOBI) {
597 // not needed by Jacobi (done by matvec)
598 Importer_ = Teuchos::rcp( new Epetra_Import(Matrix().RowMatrixColMap(),
599 Matrix().RowMatrixRowMap()) );
600
601 if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
602 }
603 IsComputed_ = true;
604 ComputeTime_ += Time_.ElapsedTime();
605 ++NumCompute_;
606
607 return(0);
608
609}
610
611//==============================================================================
612template<typename T>
615{
616 if (!IsComputed())
617 IFPACK_CHK_ERR(-3);
618
619 if (X.NumVectors() != Y.NumVectors())
620 IFPACK_CHK_ERR(-2);
621
622 Time_.ResetStartTime();
623
624 // AztecOO gives X and Y pointing to the same memory location,
625 // need to create an auxiliary vector, Xcopy
626 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
627 if (X.Pointers()[0] == Y.Pointers()[0])
628 Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
629 else
630 Xcopy = Teuchos::rcp( &X, false );
631
632 switch (PrecType_) {
633 case IFPACK_JACOBI:
634 IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
635 break;
636 case IFPACK_GS:
637 IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
638 break;
639 case IFPACK_SGS:
640 IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
641 break;
642 }
643
644 ApplyInverseTime_ += Time_.ElapsedTime();
645 ++NumApplyInverse_;
646
647 return(0);
648}
649
650//==============================================================================
651// This method in general will not work with AztecOO if used
652// outside Ifpack_AdditiveSchwarz and OverlapLevel_ != 0
653//
654template<typename T>
657 Epetra_MultiVector& Y) const
658{
659
660 if (ZeroStartingSolution_)
661 Y.PutScalar(0.0);
662
663 // do not compute the residual in this case
664 if (NumSweeps_ == 1 && ZeroStartingSolution_) {
665 int ierr = DoJacobi(X,Y);
666 IFPACK_RETURN(ierr);
667 }
668
669 Epetra_MultiVector AX(Y);
670
671 for (int j = 0; j < NumSweeps_ ; j++) {
672 IFPACK_CHK_ERR(Apply(Y,AX));
673 ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalNonzeros64();
674 IFPACK_CHK_ERR(AX.Update(1.0,X,-1.0));
675 ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalRows64();
676 IFPACK_CHK_ERR(DoJacobi(AX,Y));
677 // flops counted in DoJacobi()
678 }
679
680
681 return(0);
682}
683
684//==============================================================================
685template<typename T>
688{
689 int NumVectors = X.NumVectors();
690
691 if (OverlapLevel_ == 0) {
692
693 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
694
695 int rows = Partitioner_->NumRowsInPart(i);
696 // may happen that a partition is empty
697 if (rows == 0)
698 continue;
699
700 if(rows != 1) {
701 int LID;
702
703 // extract RHS from X
704 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
705 LID = Containers_[i]->ID(j);
706 for (int k = 0 ; k < NumVectors ; ++k) {
707 Containers_[i]->RHS(j,k) = X[k][LID];
708 }
709 }
710
711 // apply the inverse of each block. NOTE: flops occurred
712 // in ApplyInverse() of each block are summed up in method
713 // ApplyInverseFlops().
714
715 IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
716
717 // copy back into solution vector Y
718 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
719 LID = Containers_[i]->ID(j);
720 for (int k = 0 ; k < NumVectors ; ++k) {
721 Y[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
722 }
723 }
724 } //end if(rows !=1)
725 else {
726 // rows == 1, this is a singleton. compute directly.
727 int LRID = (*Partitioner_)(i,0);
728 double b = X[0][LRID];
729 double a = Diagonal_[LRID];
730 Y[0][LRID] += DampingFactor_* b/a;
731 }
732 // NOTE: flops for ApplyInverse() of each block are summed up
733 // in method ApplyInverseFlops()
734#ifdef IFPACK_FLOPCOUNTERS
735 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
736#endif
737
738 }
739 }
740 else { // overlap test
741
742 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
743
744 int rows = Partitioner_->NumRowsInPart(i);
745
746 // may happen that a partition is empty
747 if (rows == 0)
748 continue;
749 if(rows != 1) {
750 int LID;
751
752 // extract RHS from X
753 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
754 LID = Containers_[i]->ID(j);
755 for (int k = 0 ; k < NumVectors ; ++k) {
756 Containers_[i]->RHS(j,k) = (*W_)[LID] * X[k][LID];
757 }
758 }
759
760 // apply the inverse of each block
761 // if(rows != 1)
762 IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
763
764 // copy back into solution vector Y
765 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
766 LID = Containers_[i]->ID(j);
767 for (int k = 0 ; k < NumVectors ; ++k) {
768 Y[k][LID] += DampingFactor_ * (*W_)[LID] * Containers_[i]->LHS(j,k);
769 }
770 }
771 } // end if(rows != 1)
772 else { // rows == 1, this is a singleton. compute directly.
773 int LRID = (*Partitioner_)(i,0);
774 double w = (*W_)[LRID];
775 double b = w * X[0][LRID];
776 double a = Diagonal_[LRID];
777
778 Y[0][LRID] += DampingFactor_ * w * b / a;
779 }
780 }
781 // NOTE: flops for ApplyInverse() of each block are summed up
782 // in method ApplyInverseFlops()
783 // NOTE: do not count for simplicity the flops due to overlapping rows
784#ifdef IFPACK_FLOPCOUNTERS
785 ApplyInverseFlops_ += NumVectors * 4 * Matrix_->NumGlobalRows();
786#endif
787 }
788
789 return(0);
790}
791
792//==============================================================================
793template<typename T>
796 Epetra_MultiVector& Y) const
797{
798
799 if (ZeroStartingSolution_)
800 Y.PutScalar(0.0);
801
802 Epetra_MultiVector Xcopy(X);
803 for (int j = 0; j < NumSweeps_ ; j++) {
804 IFPACK_CHK_ERR(DoGaussSeidel(Xcopy,Y));
805 if (j != NumSweeps_ - 1)
806 Xcopy = X;
807 }
808
809 return(0);
810
811}
812
813//==============================================================================
814template<typename T>
817{
818
819 // cycle over all local subdomains
820
821 int Length = Matrix().MaxNumEntries();
822 std::vector<int> Indices(Length);
823 std::vector<double> Values(Length);
824
825 int NumMyRows = Matrix().NumMyRows();
826 int NumVectors = X.NumVectors();
827
828 // an additonal vector is needed by parallel computations
829 // (note that applications through Ifpack_AdditiveSchwarz
830 // are always seen are serial)
831 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
832 if (IsParallel_)
833 Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
834 else
835 Y2 = Teuchos::rcp( &Y, false );
836
837 double** y_ptr;
838 double** y2_ptr;
839 Y.ExtractView(&y_ptr);
840 Y2->ExtractView(&y2_ptr);
841
842 // data exchange is here, once per sweep
843 if (IsParallel_)
844 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
845
846 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
847 int rows = Partitioner_->NumRowsInPart(i);
848
849 // may happen that a partition is empty, but if rows == 1, the container is null
850 if (rows!=1 && Containers_[i]->NumRows() == 0)
851 continue;
852
853 int LID;
854
855 // update from previous block
856
857 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
858 LID = (*Partitioner_)(i,j);
859
860 int NumEntries;
861 IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
862 &Values[0], &Indices[0]));
863
864 for (int k = 0 ; k < NumEntries ; ++k) {
865 int col = Indices[k];
866
867 for (int kk = 0 ; kk < NumVectors ; ++kk) {
868 X[kk][LID] -= Values[k] * y2_ptr[kk][col];
869 }
870 }
871 }
872
873 if(rows != 1) {
874 // solve with this block
875
876 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
877 LID = Containers_[i]->ID(j);
878 for (int k = 0 ; k < NumVectors ; ++k) {
879 Containers_[i]->RHS(j,k) = X[k][LID];
880 }
881 }
882
883 IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
884#ifdef IFPACK_FLOPCOUNTERS
885 ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
886#endif
887
888 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
889 LID = Containers_[i]->ID(j);
890 for (int k = 0 ; k < NumVectors ; ++k) {
891 double temp = DampingFactor_ * Containers_[i]->LHS(j,k);
892 y2_ptr[k][LID] += temp;
893 }
894 }
895 } // end if(rows != 1)
896 else {
897 int LRID = (*Partitioner_)(i,0);
898 double b = X[0][LRID];
899 double a = Diagonal_[LRID];
900 y2_ptr[0][LRID]+= DampingFactor_* b/a;
901 }
902 }
903 // operations for all getrow()'s
904 // NOTE: flops for ApplyInverse() of each block are summed up
905 // in method ApplyInverseFlops()
906#ifdef IFPACK_FLOPCOUNTERS
907 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
908 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
909#endif
910
911 // Attention: this is delicate... Not all combinations
912 // of Y2 and Y will always work (tough for ML it should be ok)
913 if (IsParallel_)
914 for (int m = 0 ; m < NumVectors ; ++m)
915 for (int i = 0 ; i < NumMyRows ; ++i)
916 y_ptr[m][i] = y2_ptr[m][i];
917
918 return(0);
919 }
920
921//==============================================================================
922template<typename T>
925{
926
927 if (ZeroStartingSolution_)
928 Y.PutScalar(0.0);
929
930 Epetra_MultiVector Xcopy(X);
931 for (int j = 0; j < NumSweeps_ ; j++) {
932 IFPACK_CHK_ERR(DoSGS(X,Xcopy,Y));
933 if (j != NumSweeps_ - 1)
934 Xcopy = X;
935 }
936 return(0);
937}
938
939//==============================================================================
940template<typename T>
943 Epetra_MultiVector& Y) const
944{
945
946 int NumMyRows = Matrix().NumMyRows();
947 int NumVectors = X.NumVectors();
948
949 int Length = Matrix().MaxNumEntries();
950 std::vector<int> Indices;
951 std::vector<double> Values;
952 Indices.resize(Length);
953 Values.resize(Length);
954
955 // an additonal vector is needed by parallel computations
956 // (note that applications through Ifpack_AdditiveSchwarz
957 // are always seen are serial)
958 Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
959 if (IsParallel_)
960 Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
961 else
962 Y2 = Teuchos::rcp( &Y, false );
963
964 double** y_ptr;
965 double** y2_ptr;
966 Y.ExtractView(&y_ptr);
967 Y2->ExtractView(&y2_ptr);
968
969 // data exchange is here, once per sweep
970 if (IsParallel_)
971 IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
972
973 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
974 int rows = Partitioner_->NumRowsInPart(i);
975 // may happen that a partition is empty
976 if (rows !=1 && Containers_[i]->NumRows() == 0)
977 continue;
978
979 int LID;
980
981 // update from previous block
982
983 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
984 LID = (*Partitioner_)(i,j);
985 int NumEntries;
986 IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
987 &Values[0], &Indices[0]));
988
989 for (int k = 0 ; k < NumEntries ; ++k) {
990 int col = Indices[k];
991
992 for (int kk = 0 ; kk < NumVectors ; ++kk) {
993 Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
994 }
995 }
996 }
997
998 // solve with this block
999
1000 if(rows != 1) {
1001 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1002 LID = Containers_[i]->ID(j);
1003 for (int k = 0 ; k < NumVectors ; ++k) {
1004 Containers_[i]->RHS(j,k) = Xcopy[k][LID];
1005 }
1006 }
1007
1008 IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
1009#ifdef IFPACK_FLOPCOUNTERS
1010 ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
1011#endif
1012
1013 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1014 LID = Containers_[i]->ID(j);
1015 for (int k = 0 ; k < NumVectors ; ++k) {
1016 y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
1017 }
1018 }
1019 }
1020 else {
1021 int LRID = (*Partitioner_)(i,0);
1022 double b = Xcopy[0][LRID];
1023 double a = Diagonal_[LRID];
1024 y2_ptr[0][LRID]+= DampingFactor_* b/a;
1025 }
1026 }
1027
1028#ifdef IFPACK_FLOPCOUNTERS
1029 // operations for all getrow()'s
1030 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
1031 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
1032#endif
1033
1034 Xcopy = X;
1035
1036 for (int i = NumLocalBlocks() - 1; i >=0 ; --i) {
1037 int rows = Partitioner_->NumRowsInPart(i);
1038 if (rows != 1 &&Containers_[i]->NumRows() == 0)
1039 continue;
1040
1041 int LID;
1042
1043 // update from previous block
1044
1045 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1046 LID = (*Partitioner_)(i,j);
1047
1048 int NumEntries;
1049 IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
1050 &Values[0], &Indices[0]));
1051
1052 for (int k = 0 ; k < NumEntries ; ++k) {
1053 int col = Indices[k];
1054
1055 for (int kk = 0 ; kk < NumVectors ; ++kk) {
1056 Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
1057 }
1058 }
1059 }
1060
1061 // solve with this block
1062 if(rows != 1) {
1063 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1064 LID = Containers_[i]->ID(j);
1065 for (int k = 0 ; k < NumVectors ; ++k) {
1066 Containers_[i]->RHS(j,k) = Xcopy[k][LID];
1067 }
1068 }
1069
1070 IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
1071#ifdef IFPACK_FLOPCOUNTERS
1072 ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
1073#endif
1074
1075 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1076 LID = Containers_[i]->ID(j);
1077 for (int k = 0 ; k < NumVectors ; ++k) {
1078 y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
1079 }
1080 }
1081 }
1082 else {
1083 int LRID = (*Partitioner_)(i,0);
1084 double b = Xcopy[0][LRID];
1085 double a = Diagonal_[LRID];
1086 y2_ptr[0][LRID]+= DampingFactor_* b/a;
1087 }
1088 }
1089
1090#ifdef IFPACK_FLOPCOUNTERS
1091 // operations for all getrow()'s
1092 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
1093 ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
1094#endif
1095
1096 // Attention: this is delicate... Not all combinations
1097 // of Y2 and Y will always work (tough for ML it should be ok)
1098 if (IsParallel_)
1099 for (int m = 0 ; m < NumVectors ; ++m)
1100 for (int i = 0 ; i < NumMyRows ; ++i)
1101 y_ptr[m][i] = y2_ptr[m][i];
1102
1103 return(0);
1104}
1105
1106//==============================================================================
1107template<typename T>
1108std::ostream& Ifpack_BlockRelaxation<T>::Print(std::ostream & os) const
1109{
1110 using std::endl;
1111
1112 std::string PT;
1113 if (PrecType_ == IFPACK_JACOBI)
1114 PT = "Jacobi";
1115 else if (PrecType_ == IFPACK_GS)
1116 PT = "Gauss-Seidel";
1117 else if (PrecType_ == IFPACK_SGS)
1118 PT = "symmetric Gauss-Seidel";
1119
1120 if (!Comm().MyPID()) {
1121 os << endl;
1122 os << "================================================================================" << endl;
1123 os << "Ifpack_BlockRelaxation, " << PT << endl;
1124 os << "Sweeps = " << NumSweeps_ << endl;
1125 os << "Damping factor = " << DampingFactor_;
1126 if (ZeroStartingSolution_)
1127 os << ", using zero starting solution" << endl;
1128 else
1129 os << ", using input starting solution" << endl;
1130 os << "Number of local blocks = " << Partitioner_->NumLocalParts() << endl;
1131 //os << "Condition number estimate = " << Condest_ << endl;
1132 os << "Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
1133 os << endl;
1134 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
1135 os << "----- ------- -------------- ------------ --------" << endl;
1136 os << "Initialize() " << std::setw(5) << NumInitialize()
1137 << " " << std::setw(15) << InitializeTime()
1138 << " " << std::setw(15) << 1.0e-6 * InitializeFlops();
1139 if (InitializeTime() != 0.0)
1140 os << " " << std::setw(15) << 1.0e-6 * InitializeFlops() / InitializeTime() << endl;
1141 else
1142 os << " " << std::setw(15) << 0.0 << endl;
1143 os << "Compute() " << std::setw(5) << NumCompute()
1144 << " " << std::setw(15) << ComputeTime()
1145 << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
1146 if (ComputeTime() != 0.0)
1147 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
1148 else
1149 os << " " << std::setw(15) << 0.0 << endl;
1150 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
1151 << " " << std::setw(15) << ApplyInverseTime()
1152 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
1153 if (ApplyInverseTime() != 0.0)
1154 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
1155 else
1156 os << " " << std::setw(15) << 0.0 << endl;
1157 os << "================================================================================" << endl;
1158 os << endl;
1159 }
1160
1161 return(os);
1162}
1163
1164//==============================================================================
1165template<typename T>
1166int Ifpack_BlockRelaxation<T>::SetParameters(Teuchos::ParameterList& List)
1167{
1168 using std::cerr;
1169 using std::endl;
1170
1171 std::string PT;
1172 if (PrecType_ == IFPACK_JACOBI)
1173 PT = "Jacobi";
1174 else if (PrecType_ == IFPACK_GS)
1175 PT = "Gauss-Seidel";
1176 else if (PrecType_ == IFPACK_SGS)
1177 PT = "symmetric Gauss-Seidel";
1178
1179 PT = List.get("relaxation: type", PT);
1180
1181 if (PT == "Jacobi") {
1182 PrecType_ = IFPACK_JACOBI;
1183 }
1184 else if (PT == "Gauss-Seidel") {
1185 PrecType_ = IFPACK_GS;
1186 }
1187 else if (PT == "symmetric Gauss-Seidel") {
1188 PrecType_ = IFPACK_SGS;
1189 } else {
1190 cerr << "Option `relaxation: type' has an incorrect value ("
1191 << PT << ")'" << endl;
1192 cerr << "(file " << __FILE__ << ", line " << __LINE__ << ")" << endl;
1193 exit(EXIT_FAILURE);
1194 }
1195
1196 NumSweeps_ = List.get("relaxation: sweeps", NumSweeps_);
1197 DampingFactor_ = List.get("relaxation: damping factor",
1198 DampingFactor_);
1199 ZeroStartingSolution_ = List.get("relaxation: zero starting solution",
1200 ZeroStartingSolution_);
1201 PartitionerType_ = List.get("partitioner: type",
1202 PartitionerType_);
1203 NumLocalBlocks_ = List.get("partitioner: local parts",
1204 NumLocalBlocks_);
1205 // only Jacobi can work with overlap among local domains,
1206 OverlapLevel_ = List.get("partitioner: overlap",
1207 OverlapLevel_);
1208
1209 // check parameters
1210 if (PrecType_ != IFPACK_JACOBI)
1211 OverlapLevel_ = 0;
1212 if (NumLocalBlocks_ < 0)
1213 NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
1214 // other checks are performed in Partitioner_
1215
1216 // copy the list as each subblock's constructor will
1217 // require it later
1218 List_ = List;
1219
1220 // set the label
1221 std::string PT2;
1222 if (PrecType_ == IFPACK_JACOBI)
1223 PT2 = "BJ";
1224 else if (PrecType_ == IFPACK_GS)
1225 PT2 = "BGS";
1226 else if (PrecType_ == IFPACK_SGS)
1227 PT2 = "BSGS";
1228 Label_ = "IFPACK (" + PT2 + ", sweeps="
1229 + Ifpack_toString(NumSweeps_) + ", damping="
1230 + Ifpack_toString(DampingFactor_) + ", blocks="
1231 + Ifpack_toString(NumLocalBlocks()) + ")";
1232
1233 return(0);
1234}
1235
1236//==============================================================================
1237template<typename T>
1239{
1240 IsInitialized_ = false;
1241 Time_.ResetStartTime();
1242
1243 Graph_ = Teuchos::rcp( new Ifpack_Graph_Epetra_RowMatrix(Teuchos::rcp(&Matrix(),false)) );
1244 if (Graph_ == Teuchos::null) IFPACK_CHK_ERR(-5);
1245
1246 if (PartitionerType_ == "linear")
1247 Partitioner_ = Teuchos::rcp( new Ifpack_LinearPartitioner(&*Graph_) );
1248 else if (PartitionerType_ == "greedy")
1249 Partitioner_ = Teuchos::rcp( new Ifpack_GreedyPartitioner(&*Graph_) );
1250 else if (PartitionerType_ == "metis")
1251 Partitioner_ = Teuchos::rcp( new Ifpack_METISPartitioner(&*Graph_) );
1252 else if (PartitionerType_ == "equation")
1253 Partitioner_ = Teuchos::rcp( new Ifpack_EquationPartitioner(&*Graph_) );
1254 else if (PartitionerType_ == "user")
1255 Partitioner_ = Teuchos::rcp( new Ifpack_UserPartitioner(&*Graph_) );
1256 else if (PartitionerType_ == "line")
1257 Partitioner_ = Teuchos::rcp( new Ifpack_LinePartitioner(&Matrix()) );
1258 else
1259 IFPACK_CHK_ERR(-2);
1260
1261 if (Partitioner_ == Teuchos::null) IFPACK_CHK_ERR(-5);
1262
1263 // need to partition the graph of A
1264 IFPACK_CHK_ERR(Partitioner_->SetParameters(List_));
1265 IFPACK_CHK_ERR(Partitioner_->Compute());
1266
1267 // get actual number of partitions
1268 NumLocalBlocks_ = Partitioner_->NumLocalParts();
1269
1270 // weight of each vertex
1271 W_ = Teuchos::rcp( new Epetra_Vector(Matrix().RowMatrixRowMap()) );
1272 W_->PutScalar(0.0);
1273
1274 for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
1275
1276 for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1277 int LID = (*Partitioner_)(i,j);
1278 (*W_)[LID]++;
1279 }
1280 }
1281 W_->Reciprocal(*W_);
1282
1283 // Update Label_ if line smoothing
1284 if (PartitionerType_ == "line") {
1285 // set the label
1286 std::string PT2;
1287 if (PrecType_ == IFPACK_JACOBI)
1288 PT2 = "BJ";
1289 else if (PrecType_ == IFPACK_GS)
1290 PT2 = "BGS";
1291 else if (PrecType_ == IFPACK_SGS)
1292 PT2 = "BSGS";
1293 Label_ = "IFPACK (" + PT2 + ", auto-line, sweeps="
1294 + Ifpack_toString(NumSweeps_) + ", damping="
1295 + Ifpack_toString(DampingFactor_) + ", blocks="
1296 + Ifpack_toString(NumLocalBlocks()) + ")";
1297 }
1298
1299 InitializeTime_ += Time_.ElapsedTime();
1300 IsInitialized_ = true;
1301 ++NumInitialize_;
1302
1303 return(0);
1304}
1305
1306//==============================================================================
1307#endif // IFPACK_BLOCKPRECONDITIONER_H
virtual int NumProc() const=0
int NumVectors() const
int ExtractView(double **A, int *MyLDA) const
double ** Pointers() const
int PutScalar(double ScalarConstant)
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Comm & Comm() const=0
virtual int NumMyRows() const=0
virtual int MaxNumEntries() const=0
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
Ifpack_BlockRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_BlockRelaxation constructor with given Epetra_RowMatrix.
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual double ComputeTime() const
Returns the time spent in Compute().
std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
virtual int Initialize()
Initializes the preconditioner.
virtual double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the block Jacobi preconditioner to X, returns the result in Y.
virtual double Condest(const Ifpack_CondestType=Ifpack_Cheap, const int=1550, const double=1e-9, Epetra_RowMatrix *=0)
Computes the condition number estimate, returns its value.
int NumLocalBlocks() const
Returns the number local blocks.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully computed.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int Compute()
Computes the preconditioner.
Ifpack_EquationPartitioner: A class to decompose an Ifpack_Graph so that each block will contain all ...
Ifpack_Graph_Epetra_RowMatrix: a class to define Ifpack_Graph as a light-weight conversion of Epetra_...
Ifpack_GreedyPartitioner: A class to decompose Ifpack_Graph's using a simple greedy algorithm.
Ifpack_LinearPartitioner: A class to define linear partitions.
Ifpack_METISPartitioner: A class to decompose Ifpack_Graph's using METIS.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
Ifpack_UserPartitioner: A class to define linear partitions.