IFPACK  Development
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 
141 template<typename T>
143 
144 public:
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 
375 protected:
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
403  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
440  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_;
481  double ComputeTime_;
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 //==============================================================================
507 template<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
553 extern int ML_NODE_ID;
554 #endif
555 
556 //==============================================================================
557 template<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 //==============================================================================
718 template<typename T>
719 int 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 //==============================================================================
796 template<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 //==============================================================================
888 template<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 //==============================================================================
930 template<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 //==============================================================================
944 template<typename T>
947 {
948  IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
949  return(0);
950 }
951 
952 //==============================================================================
953 template<typename T>
955 {
956  return(-1.0);
957 }
958 
959 //==============================================================================
960 template<typename T>
962 {
963  return(Label_.c_str());
964 }
965 
966 //==============================================================================
967 template<typename T>
969 {
970  return(UseTranspose_);
971 }
972 
973 //==============================================================================
974 template<typename T>
976 {
977  return(false);
978 }
979 
980 //==============================================================================
981 template<typename T>
983 {
984  return(Matrix_->Comm());
985 }
986 
987 //==============================================================================
988 template<typename T>
990 {
991  return(Matrix_->OperatorDomainMap());
992 }
993 
994 //==============================================================================
995 template<typename T>
997 {
998  return(Matrix_->OperatorRangeMap());
999 }
1000 
1001 //==============================================================================
1002 template<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 //==============================================================================
1141 template<typename T>
1142 std::ostream& Ifpack_AdditiveSchwarz<T>::
1143 Print(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 //==============================================================================
1222 template<typename T>
1224 Condest(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
int Setup()
Sets up the localized matrix and the singleton filter.
virtual double Condest() const
Returns the estimated condition number, or -1.0 if not computed.
int OverlapLevel_
Level of overlap among the processors.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
std::string ReorderingType_
Type of reordering of the local matrix.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
virtual const Epetra_Map & RowMatrixRowMap() const=0
double ComputeFlops_
Contains the number of flops for Compute().
Copy
virtual int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
Ifpack_AMDReordering: approximate minimum degree reordering.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual std::ostream & Print(std::ostream &) const
Prints major information about this preconditioner.
virtual int OverlapLevel() const
Returns the level of overlap.
bool UseTranspose_
If true, solve with the transpose (not supported by all solvers).
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
Zero
bool ComputeCondest_
If true, compute the condition number estimate each time Compute() is called.
Insert
double Condest_
Contains the estimated condition number.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Ifpack_METISReordering: A class to reorder a graph using METIS.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
Teuchos::RefCountPtr< Ifpack_OverlappingRowMatrix > OverlappingMatrix_
Pointers to the overlapping matrix.
AbsMax
virtual int NumCompute() const
Returns the number of calls to Compute().
Average
Teuchos::RefCountPtr< Ifpack_Reordering > Reordering_
Pointer to a reordering object.
virtual const Epetra_RowMatrix & Matrix() const
Returns a refernence to the internally stored matrix.
Teuchos::RefCountPtr< Ifpack_ReorderFilter > ReorderedLocalizedMatrix_
Pointer to the reorderd matrix.
virtual const Teuchos::ParameterList & List() const
Returns a reference to the internally stored list.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual double InitializeTime() const
Returns the time spent in Initialize().
std::string Label_
Contains the label of this object.
virtual bool IsOverlapping() const
Returns true is an overlapping matrix is present.
bool IsOverlapping_
If true, overlapping is used.
virtual int Initialize()
Initialized the preconditioner.
Ifpack_AdditiveSchwarz(const Ifpack_AdditiveSchwarz &RHS)
Copy constructor (should never be used)
double ComputeTime_
Contains the time for all successful calls to Compute().
Ifpack_RCMReordering: reverse Cuthill-McKee reordering.
Ifpack_AdditiveSchwarz: a class to define Additive Schwarz preconditioners of Epetra_RowMatrix&#39;s.
int NumVectors() const
virtual const char * Label() const
Returns a character string describing the operator.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual const Epetra_BlockMap & Map() const=0
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Ifpack_AdditiveSchwarz(Epetra_RowMatrix *Matrix_in, int OverlapLevel_in=0)
Ifpack_AdditiveSchwarz constructor with given Epetra_RowMatrix.
double InitializeFlops_
Contains the number of flops for Initialize().
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
Ifpack_SingletonFilter: Filter based on matrix entries.
virtual ~Ifpack_AdditiveSchwarz()
Destructor.
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
Ifpack_LocalFilter a class for light-weight extraction of the submatrix corresponding to local rows a...
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
InsertAdd
virtual int NumInitialize() const
Returns the number of calls to Initialize().
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual int SetParameters(Teuchos::ParameterList &List)
Sets the parameters.
bool UseReordering_
If true, reorder the local matrix.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
virtual int Compute()
Computes the preconditioner.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
Teuchos::RefCountPtr< Ifpack_SingletonFilter > SingletonFilter_
filtering object.
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
bool IsComputed_
If true, the preconditioner has been successfully computed.
Epetra_CombineMode
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Teuchos::RefCountPtr< T > Inverse_
Pointer to the local solver.
int NumCompute_
Contains the number of successful call to Compute().
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to X, returns the result in Y.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Ifpack_ReorderFilter: a class for light-weight reorder of local rows and columns of an Epetra_RowMatr...
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
Add
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
Teuchos::ParameterList List_
Stores a copy of the list given in SetParameters()
Epetra_CombineMode CombineMode_
Combine mode for off-process elements (only if overlap is used)
bool FilterSingletons_
Filter for singletons.
Teuchos::RefCountPtr< Epetra_Time > Time_
Object used for timing purposes.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Teuchos::RefCountPtr< Ifpack_LocalFilter > LocalizedMatrix_
Localized version of Matrix_ or OverlappingMatrix_.