IFPACK Development
Loading...
Searching...
No Matches
Ifpack_AdditiveSchwarz.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_ADDITIVESCHWARZ_H
45#define IFPACK_ADDITIVESCHWARZ_H
46
47#include "Ifpack_ConfigDefs.h"
48#include "Ifpack_Preconditioner.h"
49#include "Ifpack_ConfigDefs.h"
50#include "Ifpack_Preconditioner.h"
51#include "Ifpack_Reordering.h"
52#include "Ifpack_RCMReordering.h"
53#include "Ifpack_METISReordering.h"
54#include "Ifpack_LocalFilter.h"
55#include "Ifpack_NodeFilter.h"
56#include "Ifpack_SingletonFilter.h"
57#include "Ifpack_ReorderFilter.h"
58#include "Ifpack_Utils.h"
59#include "Ifpack_OverlappingRowMatrix.h"
60#include "Epetra_CombineMode.h"
61#include "Epetra_MultiVector.h"
62#include "Epetra_Map.h"
63#include "Epetra_Comm.h"
64#include "Epetra_Time.h"
65#include "Epetra_LinearProblem.h"
66#include "Epetra_RowMatrix.h"
67#include "Epetra_CrsMatrix.h"
68#include "Teuchos_ParameterList.hpp"
69#include "Teuchos_RefCountPtr.hpp"
70
71#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
72#include "Ifpack_SubdomainFilter.h"
73#include "EpetraExt_Reindex_CrsMatrix.h"
74#include "EpetraExt_Reindex_MultiVector.h"
75#endif
76#ifdef IFPACK_NODE_AWARE_CODE
77#include "EpetraExt_OperatorOut.h"
78#include "EpetraExt_RowMatrixOut.h"
79#include "EpetraExt_BlockMapOut.h"
80#endif
81
82#ifdef HAVE_IFPACK_AMESOS
83 #include "Ifpack_AMDReordering.h"
84#endif
85
86
88
141template<typename T>
143
144public:
145
147
160 int OverlapLevel_in = 0);
161
165
167
169
178 virtual int SetUseTranspose(bool UseTranspose_in);
180
182
184
194 virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
195
197
208 virtual int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
209
211 virtual double NormInf() const;
213
215
217 virtual const char * Label() const;
218
220 virtual bool UseTranspose() const;
221
223 virtual bool HasNormInf() const;
224
226 virtual const Epetra_Comm & Comm() const;
227
229 virtual const Epetra_Map & OperatorDomainMap() const;
230
232 virtual const Epetra_Map & OperatorRangeMap() const;
234
236 virtual bool IsInitialized() const
237 {
238 return(IsInitialized_);
239 }
240
242 virtual bool IsComputed() const
243 {
244 return(IsComputed_);
245 }
246
248
267 virtual int SetParameters(Teuchos::ParameterList& List);
268
269 // @}
270
271 // @{ Query methods
272
274 virtual int Initialize();
275
277 virtual int Compute();
278
280 virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
281 const int MaxIters = 1550,
282 const double Tol = 1e-9,
283 Epetra_RowMatrix* Matrix_in = 0);
284
286 virtual double Condest() const
287 {
288 return(Condest_);
289 }
290
292 virtual const Epetra_RowMatrix& Matrix() const
293 {
294 return(*Matrix_);
295 }
296
298 virtual bool IsOverlapping() const
299 {
300 return(IsOverlapping_);
301 }
302
304 virtual std::ostream& Print(std::ostream&) const;
305
306 virtual const T* Inverse() const
307 {
308 return(&*Inverse_);
309 }
310
312 virtual int NumInitialize() const
313 {
314 return(NumInitialize_);
315 }
316
318 virtual int NumCompute() const
319 {
320 return(NumCompute_);
321 }
322
324 virtual int NumApplyInverse() const
325 {
326 return(NumApplyInverse_);
327 }
328
330 virtual double InitializeTime() const
331 {
332 return(InitializeTime_);
333 }
334
336 virtual double ComputeTime() const
337 {
338 return(ComputeTime_);
339 }
340
342 virtual double ApplyInverseTime() const
343 {
344 return(ApplyInverseTime_);
345 }
346
348 virtual double InitializeFlops() const
349 {
350 return(InitializeFlops_);
351 }
352
353 virtual double ComputeFlops() const
354 {
355 return(ComputeFlops_);
356 }
357
358 virtual double ApplyInverseFlops() const
359 {
360 return(ApplyInverseFlops_);
361 }
362
364 virtual int OverlapLevel() const
365 {
366 return(OverlapLevel_);
367 }
368
370 virtual const Teuchos::ParameterList& List() const
371 {
372 return(List_);
373 }
374
375protected:
376
377 // @}
378
379 // @{ Internal merhods.
380
383 { }
384
386 int Setup();
387
388 // @}
389
390 // @{ Internal data.
391
393 Teuchos::RefCountPtr<const Epetra_RowMatrix> Matrix_;
395 Teuchos::RefCountPtr<Ifpack_OverlappingRowMatrix> OverlappingMatrix_;
397/*
398 //TODO if we choose to expose the node aware code, i.e., no ifdefs,
399 //TODO then we should switch to this definition.
400 Teuchos::RefCountPtr<Epetra_RowMatrix> LocalizedMatrix_;
401*/
402#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
404 Teuchos::RCP<Epetra_RowMatrix> LocalizedMatrix_;
406 Teuchos::RCP<Epetra_CrsMatrix> SubdomainCrsMatrix_;
408 Teuchos::RCP<Epetra_CrsMatrix> ReindexedCrsMatrix_;
409
410 // The following data members are needed when doing ApplyInverse
411 // with an Ifpack_SubdomainFilter local matrix
412 Teuchos::RCP<Epetra_Map> tempMap_;
413 Teuchos::RCP<Epetra_Map> tempDomainMap_;
414 Teuchos::RCP<Epetra_Map> tempRangeMap_;
415 Teuchos::RCP<EpetraExt::CrsMatrix_Reindex> SubdomainMatrixReindexer_;
416 Teuchos::RCP<EpetraExt::MultiVector_Reindex> DomainVectorReindexer_;
417 Teuchos::RCP<EpetraExt::MultiVector_Reindex> RangeVectorReindexer_;
418 mutable Teuchos::RCP<Epetra_MultiVector> tempX_;
419 mutable Teuchos::RCP<Epetra_MultiVector> tempY_;
420#else
421# ifdef IFPACK_NODE_AWARE_CODE
422 Teuchos::RefCountPtr<Ifpack_NodeFilter> LocalizedMatrix_;
423# else
424 Teuchos::RefCountPtr<Ifpack_LocalFilter> LocalizedMatrix_;
425# endif
426#endif
427#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
428 // Some data members used for determining the subdomain id (color)
430 int MpiRank_;
432 int NumMpiProcs_;
434 int NumMpiProcsPerSubdomain_;
436 int NumSubdomains_;
438 int SubdomainId_;
439#endif
441 std::string Label_;
453 Teuchos::ParameterList List_;
457 double Condest_;
463 std::string ReorderingType_;
465 Teuchos::RefCountPtr<Ifpack_Reordering> Reordering_;
467 Teuchos::RefCountPtr<Ifpack_ReorderFilter> ReorderedLocalizedMatrix_;
471 Teuchos::RefCountPtr<Ifpack_SingletonFilter> SingletonFilter_;
477 mutable int NumApplyInverse_;
483 mutable double ApplyInverseTime_;
489 mutable double ApplyInverseFlops_;
491 Teuchos::RefCountPtr<Epetra_Time> Time_;
493 Teuchos::RefCountPtr<T> Inverse_;
495#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
496 mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
497 mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
498#endif
499# ifdef IFPACK_NODE_AWARE_CODE
500 mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
501 mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
502#endif
503
504}; // class Ifpack_AdditiveSchwarz<T>
505
506//==============================================================================
507template<typename T>
510 int OverlapLevel_in) :
511#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
512 MpiRank_(0),
513 NumMpiProcs_(1),
514 NumMpiProcsPerSubdomain_(1),
515 NumSubdomains_(1),
516 SubdomainId_(0),
517#endif
518 IsInitialized_(false),
519 IsComputed_(false),
520 UseTranspose_(false),
521 IsOverlapping_(false),
522 OverlapLevel_(OverlapLevel_in),
523 CombineMode_(Zero),
524 Condest_(-1.0),
525 ComputeCondest_(true),
526 UseReordering_(false),
527 ReorderingType_("none"),
528 FilterSingletons_(false),
529 NumInitialize_(0),
530 NumCompute_(0),
531 NumApplyInverse_(0),
532 InitializeTime_(0.0),
533 ComputeTime_(0.0),
534 ApplyInverseTime_(0.0),
535 InitializeFlops_(0.0),
536 ComputeFlops_(0.0),
537 ApplyInverseFlops_(0.0)
538{
539 // Construct a reference-counted pointer with the input matrix, don't manage the memory.
540 Matrix_ = Teuchos::rcp( Matrix_in, false );
541
542 if (Matrix_->Comm().NumProc() == 1)
543 OverlapLevel_ = 0;
544
545 if ((OverlapLevel_ != 0) && (Matrix_->Comm().NumProc() > 1))
546 IsOverlapping_ = true;
547 // Sets parameters to default values
548 Teuchos::ParameterList List_in;
549 SetParameters(List_in);
550}
551
552#ifdef IFPACK_NODE_AWARE_CODE
553extern int ML_NODE_ID;
554#endif
555
556//==============================================================================
557template<typename T>
559{
560 using std::cerr;
561 using std::endl;
562
563 Epetra_RowMatrix* MatrixPtr;
564
565#ifndef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
566# ifdef IFPACK_NODE_AWARE_CODE
567/*
568 sleep(3);
569 if (Comm().MyPID() == 0) cout << "Printing out ovArowmap" << endl;
570 Comm().Barrier();
571
572 EpetraExt::BlockMapToMatrixMarketFile("ovArowmap",OverlappingMatrix_->RowMatrixRowMap());
573 if (Comm().MyPID() == 0) cout << "Printing out ovAcolmap" << endl;
574 Comm().Barrier();
575 EpetraExt::BlockMapToMatrixMarketFile("ovAcolmap",OverlappingMatrix_->RowMatrixColMap());
576 Comm().Barrier();
577*/
578/*
579 EpetraExt::RowMatrixToMatlabFile("ovA",*OverlappingMatrix_);
580 fprintf(stderr,"p %d n %d matrix file done\n",Comm().MyPID(),ML_NODE_ID);
581 Comm().Barrier();
582*/
583 int nodeID;
584 try{ nodeID = List_.get("ML node id",0);}
585 catch(...){fprintf(stderr,"%s","Ifpack_AdditiveSchwarz<T>::Setup(): no parameter \"ML node id\"\n\n");
586 std::cout << List_ << std::endl;}
587# endif
588#endif
589
590 try{
591 if (OverlappingMatrix_ != Teuchos::null)
592 {
593#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
594 if (NumMpiProcsPerSubdomain_ > 1) {
595 LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(OverlappingMatrix_, SubdomainId_));
596 } else {
597 LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(OverlappingMatrix_));
598 }
599#else
600# ifdef IFPACK_NODE_AWARE_CODE
601 Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
602 LocalizedMatrix_ = Teuchos::rcp(tt);
603 //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
604# else
605 LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
606# endif
607#endif
608 }
609 else
610 {
611#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
612
613 if (NumMpiProcsPerSubdomain_ > 1) {
614 LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(Matrix_, SubdomainId_));
615 } else {
616 LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(Matrix_));
617 }
618#else
619# ifdef IFPACK_NODE_AWARE_CODE
620 Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(Matrix_,nodeID); //FIXME
621 LocalizedMatrix_ = Teuchos::rcp(tt);
622 //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
623# else
624 LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
625# endif
626#endif
627 }
628 }
629 catch(...) {
630 fprintf(stderr,"%s","AdditiveSchwarz Setup: problem creating local filter matrix.\n");
631 }
632
633 if (LocalizedMatrix_ == Teuchos::null)
634 IFPACK_CHK_ERR(-5);
635
636 // users may want to skip singleton check
637 if (FilterSingletons_) {
638 SingletonFilter_ = Teuchos::rcp( new Ifpack_SingletonFilter(LocalizedMatrix_) );
639 MatrixPtr = &*SingletonFilter_;
640 }
641 else
642 MatrixPtr = &*LocalizedMatrix_;
643
644 if (UseReordering_) {
645
646 // create reordering and compute it
647 if (ReorderingType_ == "rcm")
648 Reordering_ = Teuchos::rcp( new Ifpack_RCMReordering() );
649 else if (ReorderingType_ == "metis")
650 Reordering_ = Teuchos::rcp( new Ifpack_METISReordering() );
651#ifdef HAVE_IFPACK_AMESOS
652 else if (ReorderingType_ == "amd" )
653 Reordering_ = Teuchos::rcp( new Ifpack_AMDReordering() );
654#endif
655 else {
656 cerr << "reordering type not correct (" << ReorderingType_ << ")" << endl;
657 exit(EXIT_FAILURE);
658 }
659 if (Reordering_ == Teuchos::null) IFPACK_CHK_ERR(-5);
660
661 IFPACK_CHK_ERR(Reordering_->SetParameters(List_));
662 IFPACK_CHK_ERR(Reordering_->Compute(*MatrixPtr));
663
664 // now create reordered localized matrix
665 ReorderedLocalizedMatrix_ =
666 Teuchos::rcp( new Ifpack_ReorderFilter(Teuchos::rcp( MatrixPtr, false ), Reordering_) );
667
668 if (ReorderedLocalizedMatrix_ == Teuchos::null) IFPACK_CHK_ERR(-5);
669
670 MatrixPtr = &*ReorderedLocalizedMatrix_;
671 }
672
673#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
674 // The subdomain matrix needs to be reindexed by Amesos so we need to make a CrsMatrix
675 // and then reindex it with EpetraExt.
676 // The reindexing is done here because this feature is only implemented in Amesos_Klu,
677 // and it is desired to have reindexing with other direct solvers in the Amesos package
678
679 SubdomainCrsMatrix_.reset(new Epetra_CrsMatrix(Copy, MatrixPtr->RowMatrixRowMap(), -1));
680 Teuchos::RCP<Epetra_Import> tempImporter = Teuchos::rcp(new Epetra_Import(SubdomainCrsMatrix_->Map(), MatrixPtr->Map()));
681 SubdomainCrsMatrix_->Import(*MatrixPtr, *tempImporter, Insert);
682 SubdomainCrsMatrix_->FillComplete();
683
684 if (NumMpiProcsPerSubdomain_ > 1) {
685 tempMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->RowMap().NumGlobalElements(),
686 SubdomainCrsMatrix_->RowMap().NumMyElements(),
687 0, SubdomainCrsMatrix_->Comm()));
688 tempRangeMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorRangeMap().NumGlobalElements(),
689 SubdomainCrsMatrix_->OperatorRangeMap().NumMyElements(),
690 0, SubdomainCrsMatrix_->Comm()));
691 tempDomainMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorDomainMap().NumGlobalElements(),
692 SubdomainCrsMatrix_->OperatorDomainMap().NumMyElements(),
693 0, SubdomainCrsMatrix_->Comm()));
694
695 SubdomainMatrixReindexer_.reset(new EpetraExt::CrsMatrix_Reindex(*tempMap_));
696 DomainVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempDomainMap_));
697 RangeVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempRangeMap_));
698
699 ReindexedCrsMatrix_.reset(&((*SubdomainMatrixReindexer_)(*SubdomainCrsMatrix_)), false);
700
701 MatrixPtr = &*ReindexedCrsMatrix_;
702
703 Inverse_ = Teuchos::rcp( new T(&*ReindexedCrsMatrix_) );
704 } else {
705 Inverse_ = Teuchos::rcp( new T(&*SubdomainCrsMatrix_) );
706 }
707#else
708 Inverse_ = Teuchos::rcp( new T(MatrixPtr) );
709#endif
710
711 if (Inverse_ == Teuchos::null)
712 IFPACK_CHK_ERR(-5);
713
714 return(0);
715}
716
717//==============================================================================
718template<typename T>
719int Ifpack_AdditiveSchwarz<T>::SetParameters(Teuchos::ParameterList& List_in)
720{
721#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
722 MpiRank_ = Matrix_->Comm().MyPID();
723 NumMpiProcs_ = Matrix_->Comm().NumProc();
724 NumMpiProcsPerSubdomain_ = List_in.get("subdomain: number-of-processors", 1);
725 NumSubdomains_ = NumMpiProcs_ / NumMpiProcsPerSubdomain_;
726 SubdomainId_ = MpiRank_ / NumMpiProcsPerSubdomain_;
727
728 if (NumSubdomains_ == 1) {
729 IsOverlapping_ = false;
730 }
731
732#endif
733 // compute the condition number each time Compute() is invoked.
734 ComputeCondest_ = List_in.get("schwarz: compute condest", ComputeCondest_);
735 // combine mode
736 if( Teuchos::ParameterEntry *combineModeEntry = List_in.getEntryPtr("schwarz: combine mode") )
737 {
738 if( typeid(std::string) == combineModeEntry->getAny().type() )
739 {
740 std::string mode = List_in.get("schwarz: combine mode", "Add");
741 if (mode == "Add")
742 CombineMode_ = Add;
743 else if (mode == "Zero")
744 CombineMode_ = Zero;
745 else if (mode == "Insert")
746 CombineMode_ = Insert;
747 else if (mode == "InsertAdd")
748 CombineMode_ = InsertAdd;
749 else if (mode == "Average")
750 CombineMode_ = Average;
751 else if (mode == "AbsMax")
752 CombineMode_ = AbsMax;
753 else
754 {
755 TEUCHOS_TEST_FOR_EXCEPTION(
756 true,std::logic_error
757 ,"Error, The (Epetra) combine mode of \""<<mode<<"\" is not valid! Only the values"
758 " \"Add\", \"Zero\", \"Insert\", \"InsertAdd\", \"Average\", and \"AbsMax\" are accepted!"
759 );
760 }
761 }
762 else if ( typeid(Epetra_CombineMode) == combineModeEntry->getAny().type() )
763 {
764 CombineMode_ = Teuchos::any_cast<Epetra_CombineMode>(combineModeEntry->getAny());
765 }
766 else
767 {
768 // Throw exception with good error message!
769 Teuchos::getParameter<std::string>(List_in,"schwarz: combine mode");
770 }
771 }
772 else
773 {
774 // Make the default be a std::string to be consistent with the valid parameters!
775 List_in.get("schwarz: combine mode","Zero");
776 }
777 // type of reordering
778 ReorderingType_ = List_in.get("schwarz: reordering type", ReorderingType_);
779 if (ReorderingType_ == "none")
780 UseReordering_ = false;
781 else
782 UseReordering_ = true;
783 // if true, filter singletons. NOTE: the filtered matrix can still have
784 // singletons! A simple example: upper triangular matrix, if I remove
785 // the lower node, I still get a matrix with a singleton! However, filter
786 // singletons should help for PDE problems with Dirichlet BCs.
787 FilterSingletons_ = List_in.get("schwarz: filter singletons", FilterSingletons_);
788
789 // This copy may be needed by Amesos or other preconditioners.
790 List_ = List_in;
791
792 return(0);
793}
794
795//==============================================================================
796template<typename T>
798{
799 IsInitialized_ = false;
800 IsComputed_ = false; // values required
801 Condest_ = -1.0; // zero-out condest
802
803 if (Time_ == Teuchos::null)
804 Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
805
806 Time_->ResetStartTime();
807
808 // compute the overlapping matrix if necessary
809 if (IsOverlapping_) {
810#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
811 if (NumMpiProcsPerSubdomain_ > 1) {
812 OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, SubdomainId_) );
813 } else {
814 OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_));
815 }
816#else
817# ifdef IFPACK_NODE_AWARE_CODE
818 int myNodeID;
819 try{ myNodeID = List_.get("ML node id",-1);}
820 catch(...){fprintf(stderr,"pid %d: no such entry (returned %d)\n",Comm().MyPID(),myNodeID);}
821/*
822 cout << "pid " << Comm().MyPID()
823 << ": calling Ifpack_OverlappingRowMatrix with myNodeID = "
824 << myNodeID << ", OverlapLevel_ = " << OverlapLevel_ << endl;
825*/
826 OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, myNodeID) );
827# else
828 OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_) );
829# endif
830#endif
831
832 if (OverlappingMatrix_ == Teuchos::null) {
833 IFPACK_CHK_ERR(-5);
834 }
835 }
836
837# ifdef IFPACK_NODE_AWARE_CODE
838/*
839 sleep(1);
840 Comm().Barrier();
841*/
842# endif
843
844 IFPACK_CHK_ERR(Setup());
845
846# ifdef IFPACK_NODE_AWARE_CODE
847/*
848 sleep(1);
849 Comm().Barrier();
850*/
851#endif
852
853 if (Inverse_ == Teuchos::null)
854 IFPACK_CHK_ERR(-5);
855
856 if (LocalizedMatrix_ == Teuchos::null)
857 IFPACK_CHK_ERR(-5);
858
859 IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose()));
860 IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
861 IFPACK_CHK_ERR(Inverse_->Initialize());
862
863 // Label is for Aztec-like solvers
864 Label_ = "Ifpack_AdditiveSchwarz, ";
865 if (UseTranspose())
866 Label_ += ", transp";
867 Label_ += ", ov = " + Ifpack_toString(OverlapLevel_)
868 + ", local solver = \n\t\t***** `" + std::string(Inverse_->Label()) + "'";
869
870 IsInitialized_ = true;
871 ++NumInitialize_;
872 InitializeTime_ += Time_->ElapsedTime();
873
874#ifdef IFPACK_FLOPCOUNTERS
875 // count flops by summing up all the InitializeFlops() in each
876 // Inverse. Each Inverse() can only give its flops -- it acts on one
877 // process only
878 double partial = Inverse_->InitializeFlops();
879 double total;
880 Comm().SumAll(&partial, &total, 1);
881 InitializeFlops_ += total;
882#endif
883
884 return(0);
885}
886
887//==============================================================================
888template<typename T>
890{
891 if (IsInitialized() == false)
892 IFPACK_CHK_ERR(Initialize());
893
894 Time_->ResetStartTime();
895 IsComputed_ = false;
896 Condest_ = -1.0;
897
898 IFPACK_CHK_ERR(Inverse_->Compute());
899
900 IsComputed_ = true; // need this here for Condest(Ifpack_Cheap)
901 ++NumCompute_;
902 ComputeTime_ += Time_->ElapsedTime();
903
904#ifdef IFPACK_FLOPCOUNTERS
905 // sum up flops
906 double partial = Inverse_->ComputeFlops();
907 double total;
908 Comm().SumAll(&partial, &total, 1);
909 ComputeFlops_ += total;
910#endif
911
912 // reset the Label
913 std::string R = "";
914 if (UseReordering_)
915 R = ReorderingType_ + " reord, ";
916
917 if (ComputeCondest_)
918 Condest(Ifpack_Cheap);
919
920 // add Condest() to label
921 Label_ = "Ifpack_AdditiveSchwarz, ov = " + Ifpack_toString(OverlapLevel_)
922 + ", local solver = \n\t\t***** `" + std::string(Inverse_->Label()) + "'"
923 + "\n\t\t***** " + R + "Condition number estimate = "
924 + Ifpack_toString(Condest_);
925
926 return(0);
927}
928
929//==============================================================================
930template<typename T>
932{
933 // store the flag -- it will be set in Initialize() if Inverse_ does not
934 // exist.
935 UseTranspose_ = UseTranspose_in;
936
937 // If Inverse_ exists, pass it right now.
938 if (Inverse_!=Teuchos::null)
939 IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
940 return(0);
941}
942
943//==============================================================================
944template<typename T>
947{
948 IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
949 return(0);
950}
951
952//==============================================================================
953template<typename T>
955{
956 return(-1.0);
957}
958
959//==============================================================================
960template<typename T>
962{
963 return(Label_.c_str());
964}
965
966//==============================================================================
967template<typename T>
969{
970 return(UseTranspose_);
971}
972
973//==============================================================================
974template<typename T>
976{
977 return(false);
978}
979
980//==============================================================================
981template<typename T>
983{
984 return(Matrix_->Comm());
985}
986
987//==============================================================================
988template<typename T>
990{
991 return(Matrix_->OperatorDomainMap());
992}
993
994//==============================================================================
995template<typename T>
997{
998 return(Matrix_->OperatorRangeMap());
999}
1000
1001//==============================================================================
1002template<typename T>
1005{
1006 // compute the preconditioner is not done by the user
1007 if (!IsComputed())
1008 IFPACK_CHK_ERR(-3);
1009
1010 int NumVectors = X.NumVectors();
1011
1012 if (NumVectors != Y.NumVectors())
1013 IFPACK_CHK_ERR(-2); // wrong input
1014
1015 Time_->ResetStartTime();
1016
1017 Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
1018 Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
1019 Teuchos::RefCountPtr<Epetra_MultiVector> Xtmp;
1020
1021 // for flop count, see bottom of this function
1022#ifdef IFPACK_FLOPCOUNTERS
1023 double pre_partial = Inverse_->ApplyInverseFlops();
1024 double pre_total;
1025 Comm().SumAll(&pre_partial, &pre_total, 1);
1026#endif
1027
1028 // process overlap, may need to create vectors and import data
1029 if (IsOverlapping()) {
1030#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
1031 if (OverlappingX == Teuchos::null) {
1032 OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), X.NumVectors()) );
1033 if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
1034 } else assert(OverlappingX->NumVectors() == X.NumVectors());
1035 if (OverlappingY == Teuchos::null) {
1036 OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), Y.NumVectors()) );
1037 if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
1038 } else assert(OverlappingY->NumVectors() == Y.NumVectors());
1039#else
1040# ifdef IFPACK_NODE_AWARE_CODE
1041 if (OverlappingX == Teuchos::null) {
1042 OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1043 X.NumVectors()) );
1044 if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
1045 } else assert(OverlappingX->NumVectors() == X.NumVectors());
1046 if (OverlappingY == Teuchos::null) {
1047 OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1048 Y.NumVectors()) );
1049 if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
1050 } else assert(OverlappingY->NumVectors() == Y.NumVectors());
1051#else
1052 OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1053 X.NumVectors()) );
1054 OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1055 Y.NumVectors()) );
1056 if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
1057# endif
1058#endif
1059 OverlappingY->PutScalar(0.0);
1060 OverlappingX->PutScalar(0.0);
1061 IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(X,*OverlappingX,Insert));
1062 // FIXME: this will not work with overlapping and non-zero starting
1063 // solutions. The same for other cases below.
1064 // IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(Y,*OverlappingY,Insert));
1065 }
1066 else {
1067 Xtmp = Teuchos::rcp( new Epetra_MultiVector(X) );
1068 OverlappingX = Xtmp;
1069 OverlappingY = Teuchos::rcp( &Y, false );
1070 }
1071
1072 if (FilterSingletons_) {
1073 // process singleton filter
1074 Epetra_MultiVector ReducedX(SingletonFilter_->Map(),NumVectors);
1075 Epetra_MultiVector ReducedY(SingletonFilter_->Map(),NumVectors);
1076 IFPACK_CHK_ERR(SingletonFilter_->SolveSingletons(*OverlappingX,*OverlappingY));
1077 IFPACK_CHK_ERR(SingletonFilter_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX));
1078
1079 // process reordering
1080 if (!UseReordering_) {
1081 IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReducedX,ReducedY));
1082 }
1083 else {
1084 Epetra_MultiVector ReorderedX(ReducedX);
1085 Epetra_MultiVector ReorderedY(ReducedY);
1086 IFPACK_CHK_ERR(Reordering_->P(ReducedX,ReorderedX));
1087 IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
1088 IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,ReducedY));
1089 }
1090
1091 // finish up with singletons
1092 IFPACK_CHK_ERR(SingletonFilter_->UpdateLHS(ReducedY,*OverlappingY));
1093 }
1094 else {
1095 // process reordering
1096 if (!UseReordering_) {
1097#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
1098 if (NumMpiProcsPerSubdomain_ > 1) {
1099 tempX_.reset(&((*RangeVectorReindexer_)(*OverlappingX)), false);
1100 tempY_.reset(&((*DomainVectorReindexer_)(*OverlappingY)), false);
1101 IFPACK_CHK_ERR(Inverse_->ApplyInverse(*tempX_,*tempY_));
1102 } else {
1103 IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX, *OverlappingY));
1104 }
1105#else
1106 IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
1107#endif
1108 }
1109 else {
1110 Epetra_MultiVector ReorderedX(*OverlappingX);
1111 Epetra_MultiVector ReorderedY(*OverlappingY);
1112 IFPACK_CHK_ERR(Reordering_->P(*OverlappingX,ReorderedX));
1113 IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
1114 IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,*OverlappingY));
1115 }
1116 }
1117
1118 if (IsOverlapping()) {
1119 IFPACK_CHK_ERR(OverlappingMatrix_->ExportMultiVector(*OverlappingY,Y,
1120 CombineMode_));
1121 }
1122
1123#ifdef IFPACK_FLOPCOUNTERS
1124 // add flops. Note the we only have to add the newly counted
1125 // flops -- and each Inverse returns the cumulative sum
1126 double partial = Inverse_->ApplyInverseFlops();
1127 double total;
1128 Comm().SumAll(&partial, &total, 1);
1129 ApplyInverseFlops_ += total - pre_total;
1130#endif
1131
1132 // FIXME: right now I am skipping the overlap and singletons
1133 ++NumApplyInverse_;
1134 ApplyInverseTime_ += Time_->ElapsedTime();
1135
1136 return(0);
1137
1138}
1139
1140//==============================================================================
1141template<typename T>
1143Print(std::ostream& os) const
1144{
1145 using std::endl;
1146
1147#ifdef IFPACK_FLOPCOUNTERS
1148 double IF = InitializeFlops();
1149 double CF = ComputeFlops();
1150 double AF = ApplyInverseFlops();
1151
1152 double IFT = 0.0, CFT = 0.0, AFT = 0.0;
1153 if (InitializeTime() != 0.0)
1154 IFT = IF / InitializeTime();
1155 if (ComputeTime() != 0.0)
1156 CFT = CF / ComputeTime();
1157 if (ApplyInverseTime() != 0.0)
1158 AFT = AF / ApplyInverseTime();
1159#endif
1160
1161 if (Matrix().Comm().MyPID())
1162 return(os);
1163
1164 os << endl;
1165 os << "================================================================================" << endl;
1166 os << "Ifpack_AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
1167 if (CombineMode_ == Insert)
1168 os << "Combine mode = Insert" << endl;
1169 else if (CombineMode_ == Add)
1170 os << "Combine mode = Add" << endl;
1171 else if (CombineMode_ == Zero)
1172 os << "Combine mode = Zero" << endl;
1173 else if (CombineMode_ == Average)
1174 os << "Combine mode = Average" << endl;
1175 else if (CombineMode_ == AbsMax)
1176 os << "Combine mode = AbsMax" << endl;
1177
1178 os << "Condition number estimate = " << Condest_ << endl;
1179 os << "Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
1180
1181#ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
1182 os << endl;
1183 os << "================================================================================" << endl;
1184 os << "Subcommunicator stats" << endl;
1185 os << "Number of MPI processes in simulation: " << NumMpiProcs_ << endl;
1186 os << "Number of subdomains: " << NumSubdomains_ << endl;
1187 os << "Number of MPI processes per subdomain: " << NumMpiProcsPerSubdomain_ << endl;
1188#endif
1189
1190 os << endl;
1191 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
1192 os << "----- ------- -------------- ------------ --------" << endl;
1193 os << "Initialize() " << std::setw(5) << NumInitialize()
1194 << " " << std::setw(15) << InitializeTime()
1195#ifdef IFPACK_FLOPCOUNTERS
1196 << " " << std::setw(15) << 1.0e-6 * IF
1197 << " " << std::setw(15) << 1.0e-6 * IFT
1198#endif
1199 << endl;
1200 os << "Compute() " << std::setw(5) << NumCompute()
1201 << " " << std::setw(15) << ComputeTime()
1202#ifdef IFPACK_FLOPCOUNTERS
1203 << " " << std::setw(15) << 1.0e-6 * CF
1204 << " " << std::setw(15) << 1.0e-6 * CFT
1205#endif
1206 << endl;
1207 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
1208 << " " << std::setw(15) << ApplyInverseTime()
1209#ifdef IFPACK_FLOPCOUNTERS
1210 << " " << std::setw(15) << 1.0e-6 * AF
1211 << " " << std::setw(15) << 1.0e-6 * AFT
1212#endif
1213 << endl;
1214 os << "================================================================================" << endl;
1215 os << endl;
1216
1217 return(os);
1218}
1219
1220#include "Ifpack_Condest.h"
1221//==============================================================================
1222template<typename T>
1224Condest(const Ifpack_CondestType CT, const int MaxIters,
1225 const double Tol, Epetra_RowMatrix* Matrix_in)
1226{
1227 if (!IsComputed()) // cannot compute right now
1228 return(-1.0);
1229
1230 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
1231
1232 return(Condest_);
1233}
1234
1235#endif // IFPACK_ADDITIVESCHWARZ_H
Epetra_CombineMode
Insert
Add
AbsMax
InsertAdd
Average
Zero
Copy
int NumVectors() const
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual const Epetra_BlockMap & Map() const=0
Ifpack_AMDReordering: approximate minimum degree reordering.
Ifpack_AdditiveSchwarz: a class to define Additive Schwarz preconditioners of Epetra_RowMatrix's.
int NumCompute_
Contains the number of successful call to Compute().
Ifpack_AdditiveSchwarz(Epetra_RowMatrix *Matrix_in, int OverlapLevel_in=0)
Ifpack_AdditiveSchwarz constructor with given Epetra_RowMatrix.
Teuchos::RefCountPtr< Ifpack_Reordering > Reordering_
Pointer to a reordering object.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
bool FilterSingletons_
Filter for singletons.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual const Teuchos::ParameterList & List() const
Returns a reference to the internally stored list.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
bool ComputeCondest_
If true, compute the condition number estimate each time Compute() is called.
double InitializeFlops_
Contains the number of flops for Initialize().
virtual const Epetra_RowMatrix & Matrix() const
Returns a refernence to the internally stored matrix.
bool IsComputed_
If true, the preconditioner has been successfully computed.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
Ifpack_AdditiveSchwarz(const Ifpack_AdditiveSchwarz &RHS)
Copy constructor (should never be used)
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
int NumInitialize_
Contains the number of successful calls to Initialize().
Epetra_CombineMode CombineMode_
Combine mode for off-process elements (only if overlap is used)
std::string Label_
Contains the label of this object.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets the parameters.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Teuchos::RefCountPtr< Ifpack_OverlappingRowMatrix > OverlappingMatrix_
Pointers to the overlapping matrix.
int OverlapLevel_
Level of overlap among the processors.
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
Teuchos::ParameterList List_
Stores a copy of the list given in SetParameters()
std::string ReorderingType_
Type of reordering of the local matrix.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
double ComputeFlops_
Contains the number of flops for Compute().
Teuchos::RefCountPtr< Ifpack_ReorderFilter > ReorderedLocalizedMatrix_
Pointer to the reorderd matrix.
virtual int Compute()
Computes the preconditioner.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
virtual int Initialize()
Initialized the preconditioner.
Teuchos::RefCountPtr< Epetra_Time > Time_
Object used for timing purposes.
bool UseTranspose_
If true, solve with the transpose (not supported by all solvers).
Teuchos::RefCountPtr< Ifpack_SingletonFilter > SingletonFilter_
filtering object.
virtual bool IsOverlapping() const
Returns true is an overlapping matrix is present.
Teuchos::RefCountPtr< T > Inverse_
Pointer to the local solver.
virtual int OverlapLevel() const
Returns the level of overlap.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual std::ostream & Print(std::ostream &) const
Prints major information about this preconditioner.
virtual const char * Label() const
Returns a character string describing the operator.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
double Condest_
Contains the estimated condition number.
double ComputeTime_
Contains the time for all successful calls to Compute().
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to X, returns the result in Y.
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.
Teuchos::RefCountPtr< Ifpack_LocalFilter > LocalizedMatrix_
Localized version of Matrix_ or OverlappingMatrix_.
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
virtual ~Ifpack_AdditiveSchwarz()
Destructor.
bool IsOverlapping_
If true, overlapping is used.
bool UseReordering_
If true, reorder the local matrix.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
int Setup()
Sets up the localized matrix and the singleton filter.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
virtual double Condest() const
Returns the estimated condition number, or -1.0 if not computed.
Ifpack_LocalFilter a class for light-weight extraction of the submatrix corresponding to local rows a...
Ifpack_METISReordering: A class to reorder a graph using METIS.
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
Ifpack_RCMReordering: reverse Cuthill-McKee reordering.
Ifpack_ReorderFilter: a class for light-weight reorder of local rows and columns of an Epetra_RowMatr...
Ifpack_SingletonFilter: Filter based on matrix entries.