IFPACK  Development
 All Classes Files Functions Variables Enumerations Friends Pages
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 
159  Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
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  string Label_;
453  Teuchos::ParameterList List_;
455  Epetra_CombineMode CombineMode_;
457  double Condest_;
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>
509 Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
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 
561  Epetra_RowMatrix* MatrixPtr;
562 
563 #ifndef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
564 # ifdef IFPACK_NODE_AWARE_CODE
565 /*
566  sleep(3);
567  if (Comm().MyPID() == 0) cout << "Printing out ovArowmap" << endl;
568  Comm().Barrier();
569 
570  EpetraExt::BlockMapToMatrixMarketFile("ovArowmap",OverlappingMatrix_->RowMatrixRowMap());
571  if (Comm().MyPID() == 0) cout << "Printing out ovAcolmap" << endl;
572  Comm().Barrier();
573  EpetraExt::BlockMapToMatrixMarketFile("ovAcolmap",OverlappingMatrix_->RowMatrixColMap());
574  Comm().Barrier();
575 */
576 /*
577  EpetraExt::RowMatrixToMatlabFile("ovA",*OverlappingMatrix_);
578  fprintf(stderr,"p %d n %d matrix file done\n",Comm().MyPID(),ML_NODE_ID);
579  Comm().Barrier();
580 */
581  int nodeID;
582  try{ nodeID = List_.get("ML node id",0);}
583  catch(...){fprintf(stderr,"%s","Ifpack_AdditiveSchwarz<T>::Setup(): no parameter \"ML node id\"\n\n");
584  cout << List_ << endl;}
585 # endif
586 #endif
587 
588  try{
589  if (OverlappingMatrix_ != Teuchos::null)
590  {
591 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
592  if (NumMpiProcsPerSubdomain_ > 1) {
593  LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(OverlappingMatrix_, SubdomainId_));
594  } else {
595  LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(OverlappingMatrix_));
596  }
597 #else
598 # ifdef IFPACK_NODE_AWARE_CODE
599  Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
600  LocalizedMatrix_ = Teuchos::rcp(tt);
601  //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
602 # else
603  LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
604 # endif
605 #endif
606  }
607  else
608  {
609 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
610 
611  if (NumMpiProcsPerSubdomain_ > 1) {
612  LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(Matrix_, SubdomainId_));
613  } else {
614  LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(Matrix_));
615  }
616 #else
617 # ifdef IFPACK_NODE_AWARE_CODE
618  Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(Matrix_,nodeID); //FIXME
619  LocalizedMatrix_ = Teuchos::rcp(tt);
620  //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
621 # else
622  LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
623 # endif
624 #endif
625  }
626  }
627  catch(...) {
628  fprintf(stderr,"%s","AdditiveSchwarz Setup: problem creating local filter matrix.\n");
629  }
630 
631  if (LocalizedMatrix_ == Teuchos::null)
632  IFPACK_CHK_ERR(-5);
633 
634  // users may want to skip singleton check
635  if (FilterSingletons_) {
636  SingletonFilter_ = Teuchos::rcp( new Ifpack_SingletonFilter(LocalizedMatrix_) );
637  MatrixPtr = &*SingletonFilter_;
638  }
639  else
640  MatrixPtr = &*LocalizedMatrix_;
641 
642  if (UseReordering_) {
643 
644  // create reordering and compute it
645  if (ReorderingType_ == "rcm")
646  Reordering_ = Teuchos::rcp( new Ifpack_RCMReordering() );
647  else if (ReorderingType_ == "metis")
648  Reordering_ = Teuchos::rcp( new Ifpack_METISReordering() );
649 #ifdef HAVE_IFPACK_AMESOS
650  else if (ReorderingType_ == "amd" )
651  Reordering_ = Teuchos::rcp( new Ifpack_AMDReordering() );
652 #endif
653  else {
654  cerr << "reordering type not correct (" << ReorderingType_ << ")" << endl;
655  exit(EXIT_FAILURE);
656  }
657  if (Reordering_ == Teuchos::null) IFPACK_CHK_ERR(-5);
658 
659  IFPACK_CHK_ERR(Reordering_->SetParameters(List_));
660  IFPACK_CHK_ERR(Reordering_->Compute(*MatrixPtr));
661 
662  // now create reordered localized matrix
663  ReorderedLocalizedMatrix_ =
664  Teuchos::rcp( new Ifpack_ReorderFilter(Teuchos::rcp( MatrixPtr, false ), Reordering_) );
665 
666  if (ReorderedLocalizedMatrix_ == Teuchos::null) IFPACK_CHK_ERR(-5);
667 
668  MatrixPtr = &*ReorderedLocalizedMatrix_;
669  }
670 
671 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
672  // The subdomain matrix needs to be reindexed by Amesos so we need to make a CrsMatrix
673  // and then reindex it with EpetraExt.
674  // The reindexing is done here because this feature is only implemented in Amesos_Klu,
675  // and it is desired to have reindexing with other direct solvers in the Amesos package
676 
677  SubdomainCrsMatrix_.reset(new Epetra_CrsMatrix(Copy, MatrixPtr->RowMatrixRowMap(), -1));
678  Teuchos::RCP<Epetra_Import> tempImporter = Teuchos::rcp(new Epetra_Import(SubdomainCrsMatrix_->Map(), MatrixPtr->Map()));
679  SubdomainCrsMatrix_->Import(*MatrixPtr, *tempImporter, Insert);
680  SubdomainCrsMatrix_->FillComplete();
681 
682  if (NumMpiProcsPerSubdomain_ > 1) {
683  tempMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->RowMap().NumGlobalElements(),
684  SubdomainCrsMatrix_->RowMap().NumMyElements(),
685  0, SubdomainCrsMatrix_->Comm()));
686  tempRangeMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorRangeMap().NumGlobalElements(),
687  SubdomainCrsMatrix_->OperatorRangeMap().NumMyElements(),
688  0, SubdomainCrsMatrix_->Comm()));
689  tempDomainMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorDomainMap().NumGlobalElements(),
690  SubdomainCrsMatrix_->OperatorDomainMap().NumMyElements(),
691  0, SubdomainCrsMatrix_->Comm()));
692 
693  SubdomainMatrixReindexer_.reset(new EpetraExt::CrsMatrix_Reindex(*tempMap_));
694  DomainVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempDomainMap_));
695  RangeVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempRangeMap_));
696 
697  ReindexedCrsMatrix_.reset(&((*SubdomainMatrixReindexer_)(*SubdomainCrsMatrix_)), false);
698 
699  MatrixPtr = &*ReindexedCrsMatrix_;
700 
701  Inverse_ = Teuchos::rcp( new T(&*ReindexedCrsMatrix_) );
702  } else {
703  Inverse_ = Teuchos::rcp( new T(&*SubdomainCrsMatrix_) );
704  }
705 #else
706  Inverse_ = Teuchos::rcp( new T(MatrixPtr) );
707 #endif
708 
709  if (Inverse_ == Teuchos::null)
710  IFPACK_CHK_ERR(-5);
711 
712  return(0);
713 }
714 
715 //==============================================================================
716 template<typename T>
717 int Ifpack_AdditiveSchwarz<T>::SetParameters(Teuchos::ParameterList& List_in)
718 {
719 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
720  MpiRank_ = Matrix_->Comm().MyPID();
721  NumMpiProcs_ = Matrix_->Comm().NumProc();
722  NumMpiProcsPerSubdomain_ = List_in.get("subdomain: number-of-processors", 1);
723  NumSubdomains_ = NumMpiProcs_ / NumMpiProcsPerSubdomain_;
724  SubdomainId_ = MpiRank_ / NumMpiProcsPerSubdomain_;
725 
726  if (NumSubdomains_ == 1) {
727  IsOverlapping_ = false;
728  }
729 
730 #endif
731  // compute the condition number each time Compute() is invoked.
732  ComputeCondest_ = List_in.get("schwarz: compute condest", ComputeCondest_);
733  // combine mode
734  if( Teuchos::ParameterEntry *combineModeEntry = List_in.getEntryPtr("schwarz: combine mode") )
735  {
736  if( typeid(std::string) == combineModeEntry->getAny().type() )
737  {
738  std::string mode = List_in.get("schwarz: combine mode", "Add");
739  if (mode == "Add")
740  CombineMode_ = Add;
741  else if (mode == "Zero")
742  CombineMode_ = Zero;
743  else if (mode == "Insert")
744  CombineMode_ = Insert;
745  else if (mode == "InsertAdd")
746  CombineMode_ = InsertAdd;
747  else if (mode == "Average")
748  CombineMode_ = Average;
749  else if (mode == "AbsMax")
750  CombineMode_ = AbsMax;
751  else
752  {
753  TEUCHOS_TEST_FOR_EXCEPTION(
754  true,std::logic_error
755  ,"Error, The (Epetra) combine mode of \""<<mode<<"\" is not valid! Only the values"
756  " \"Add\", \"Zero\", \"Insert\", \"InsertAdd\", \"Average\", and \"AbsMax\" are accepted!"
757  );
758  }
759  }
760  else if ( typeid(Epetra_CombineMode) == combineModeEntry->getAny().type() )
761  {
762  CombineMode_ = Teuchos::any_cast<Epetra_CombineMode>(combineModeEntry->getAny());
763  }
764  else
765  {
766  // Throw exception with good error message!
767  Teuchos::getParameter<std::string>(List_in,"schwarz: combine mode");
768  }
769  }
770  else
771  {
772  // Make the default be a string to be consistent with the valid parameters!
773  List_in.get("schwarz: combine mode","Zero");
774  }
775  // type of reordering
776  ReorderingType_ = List_in.get("schwarz: reordering type", ReorderingType_);
777  if (ReorderingType_ == "none")
778  UseReordering_ = false;
779  else
780  UseReordering_ = true;
781  // if true, filter singletons. NOTE: the filtered matrix can still have
782  // singletons! A simple example: upper triangular matrix, if I remove
783  // the lower node, I still get a matrix with a singleton! However, filter
784  // singletons should help for PDE problems with Dirichlet BCs.
785  FilterSingletons_ = List_in.get("schwarz: filter singletons", FilterSingletons_);
786 
787  // This copy may be needed by Amesos or other preconditioners.
788  List_ = List_in;
789 
790  return(0);
791 }
792 
793 //==============================================================================
794 template<typename T>
796 {
797  IsInitialized_ = false;
798  IsComputed_ = false; // values required
799  Condest_ = -1.0; // zero-out condest
800 
801  if (Time_ == Teuchos::null)
802  Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
803 
804  Time_->ResetStartTime();
805 
806  // compute the overlapping matrix if necessary
807  if (IsOverlapping_) {
808 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
809  if (NumMpiProcsPerSubdomain_ > 1) {
810  OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, SubdomainId_) );
811  } else {
812  OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_));
813  }
814 #else
815 # ifdef IFPACK_NODE_AWARE_CODE
816  int myNodeID;
817  try{ myNodeID = List_.get("ML node id",-1);}
818  catch(...){fprintf(stderr,"pid %d: no such entry (returned %d)\n",Comm().MyPID(),myNodeID);}
819 /*
820  cout << "pid " << Comm().MyPID()
821  << ": calling Ifpack_OverlappingRowMatrix with myNodeID = "
822  << myNodeID << ", OverlapLevel_ = " << OverlapLevel_ << endl;
823 */
824  OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, myNodeID) );
825 # else
826  OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_) );
827 # endif
828 #endif
829 
830  if (OverlappingMatrix_ == Teuchos::null) {
831  IFPACK_CHK_ERR(-5);
832  }
833  }
834 
835 # ifdef IFPACK_NODE_AWARE_CODE
836 /*
837  sleep(1);
838  Comm().Barrier();
839 */
840 # endif
841 
842  IFPACK_CHK_ERR(Setup());
843 
844 # ifdef IFPACK_NODE_AWARE_CODE
845 /*
846  sleep(1);
847  Comm().Barrier();
848 */
849 #endif
850 
851  if (Inverse_ == Teuchos::null)
852  IFPACK_CHK_ERR(-5);
853 
854  if (LocalizedMatrix_ == Teuchos::null)
855  IFPACK_CHK_ERR(-5);
856 
857  IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose()));
858  IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
859  IFPACK_CHK_ERR(Inverse_->Initialize());
860 
861  // Label is for Aztec-like solvers
862  Label_ = "Ifpack_AdditiveSchwarz, ";
863  if (UseTranspose())
864  Label_ += ", transp";
865  Label_ += ", ov = " + Ifpack_toString(OverlapLevel_)
866  + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'";
867 
868  IsInitialized_ = true;
869  ++NumInitialize_;
870  InitializeTime_ += Time_->ElapsedTime();
871 
872 #ifdef IFPACK_FLOPCOUNTERS
873  // count flops by summing up all the InitializeFlops() in each
874  // Inverse. Each Inverse() can only give its flops -- it acts on one
875  // process only
876  double partial = Inverse_->InitializeFlops();
877  double total;
878  Comm().SumAll(&partial, &total, 1);
879  InitializeFlops_ += total;
880 #endif
881 
882  return(0);
883 }
884 
885 //==============================================================================
886 template<typename T>
888 {
889  if (IsInitialized() == false)
890  IFPACK_CHK_ERR(Initialize());
891 
892  Time_->ResetStartTime();
893  IsComputed_ = false;
894  Condest_ = -1.0;
895 
896  IFPACK_CHK_ERR(Inverse_->Compute());
897 
898  IsComputed_ = true; // need this here for Condest(Ifpack_Cheap)
899  ++NumCompute_;
900  ComputeTime_ += Time_->ElapsedTime();
901 
902 #ifdef IFPACK_FLOPCOUNTERS
903  // sum up flops
904  double partial = Inverse_->ComputeFlops();
905  double total;
906  Comm().SumAll(&partial, &total, 1);
907  ComputeFlops_ += total;
908 #endif
909 
910  // reset the Label
911  string R = "";
912  if (UseReordering_)
913  R = ReorderingType_ + " reord, ";
914 
915  if (ComputeCondest_)
916  Condest(Ifpack_Cheap);
917 
918  // add Condest() to label
919  Label_ = "Ifpack_AdditiveSchwarz, ov = " + Ifpack_toString(OverlapLevel_)
920  + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'"
921  + "\n\t\t***** " + R + "Condition number estimate = "
922  + Ifpack_toString(Condest_);
923 
924  return(0);
925 }
926 
927 //==============================================================================
928 template<typename T>
930 {
931  // store the flag -- it will be set in Initialize() if Inverse_ does not
932  // exist.
933  UseTranspose_ = UseTranspose_in;
934 
935  // If Inverse_ exists, pass it right now.
936  if (Inverse_!=Teuchos::null)
937  IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
938  return(0);
939 }
940 
941 //==============================================================================
942 template<typename T>
944 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
945 {
946  IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
947  return(0);
948 }
949 
950 //==============================================================================
951 template<typename T>
953 {
954  return(-1.0);
955 }
956 
957 //==============================================================================
958 template<typename T>
960 {
961  return(Label_.c_str());
962 }
963 
964 //==============================================================================
965 template<typename T>
967 {
968  return(UseTranspose_);
969 }
970 
971 //==============================================================================
972 template<typename T>
974 {
975  return(false);
976 }
977 
978 //==============================================================================
979 template<typename T>
980 const Epetra_Comm & Ifpack_AdditiveSchwarz<T>::Comm() const
981 {
982  return(Matrix_->Comm());
983 }
984 
985 //==============================================================================
986 template<typename T>
988 {
989  return(Matrix_->OperatorDomainMap());
990 }
991 
992 //==============================================================================
993 template<typename T>
995 {
996  return(Matrix_->OperatorRangeMap());
997 }
998 
999 //==============================================================================
1000 template<typename T>
1002 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
1003 {
1004  // compute the preconditioner is not done by the user
1005  if (!IsComputed())
1006  IFPACK_CHK_ERR(-3);
1007 
1008  int NumVectors = X.NumVectors();
1009 
1010  if (NumVectors != Y.NumVectors())
1011  IFPACK_CHK_ERR(-2); // wrong input
1012 
1013  Time_->ResetStartTime();
1014 
1015  Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
1016  Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
1017  Teuchos::RefCountPtr<Epetra_MultiVector> Xtmp;
1018 
1019  // for flop count, see bottom of this function
1020 #ifdef IFPACK_FLOPCOUNTERS
1021  double pre_partial = Inverse_->ApplyInverseFlops();
1022  double pre_total;
1023  Comm().SumAll(&pre_partial, &pre_total, 1);
1024 #endif
1025 
1026  // process overlap, may need to create vectors and import data
1027  if (IsOverlapping()) {
1028 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
1029  if (OverlappingX == Teuchos::null) {
1030  OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), X.NumVectors()) );
1031  if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
1032  } else assert(OverlappingX->NumVectors() == X.NumVectors());
1033  if (OverlappingY == Teuchos::null) {
1034  OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), Y.NumVectors()) );
1035  if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
1036  } else assert(OverlappingY->NumVectors() == Y.NumVectors());
1037 #else
1038 # ifdef IFPACK_NODE_AWARE_CODE
1039  if (OverlappingX == Teuchos::null) {
1040  OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1041  X.NumVectors()) );
1042  if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
1043  } else assert(OverlappingX->NumVectors() == X.NumVectors());
1044  if (OverlappingY == Teuchos::null) {
1045  OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1046  Y.NumVectors()) );
1047  if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
1048  } else assert(OverlappingY->NumVectors() == Y.NumVectors());
1049 #else
1050  OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1051  X.NumVectors()) );
1052  OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
1053  Y.NumVectors()) );
1054  if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
1055 # endif
1056 #endif
1057  OverlappingY->PutScalar(0.0);
1058  OverlappingX->PutScalar(0.0);
1059  IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(X,*OverlappingX,Insert));
1060  // FIXME: this will not work with overlapping and non-zero starting
1061  // solutions. The same for other cases below.
1062  // IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(Y,*OverlappingY,Insert));
1063  }
1064  else {
1065  Xtmp = Teuchos::rcp( new Epetra_MultiVector(X) );
1066  OverlappingX = Xtmp;
1067  OverlappingY = Teuchos::rcp( &Y, false );
1068  }
1069 
1070  if (FilterSingletons_) {
1071  // process singleton filter
1072  Epetra_MultiVector ReducedX(SingletonFilter_->Map(),NumVectors);
1073  Epetra_MultiVector ReducedY(SingletonFilter_->Map(),NumVectors);
1074  IFPACK_CHK_ERR(SingletonFilter_->SolveSingletons(*OverlappingX,*OverlappingY));
1075  IFPACK_CHK_ERR(SingletonFilter_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX));
1076 
1077  // process reordering
1078  if (!UseReordering_) {
1079  IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReducedX,ReducedY));
1080  }
1081  else {
1082  Epetra_MultiVector ReorderedX(ReducedX);
1083  Epetra_MultiVector ReorderedY(ReducedY);
1084  IFPACK_CHK_ERR(Reordering_->P(ReducedX,ReorderedX));
1085  IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
1086  IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,ReducedY));
1087  }
1088 
1089  // finish up with singletons
1090  IFPACK_CHK_ERR(SingletonFilter_->UpdateLHS(ReducedY,*OverlappingY));
1091  }
1092  else {
1093  // process reordering
1094  if (!UseReordering_) {
1095 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
1096  if (NumMpiProcsPerSubdomain_ > 1) {
1097  tempX_.reset(&((*RangeVectorReindexer_)(*OverlappingX)), false);
1098  tempY_.reset(&((*DomainVectorReindexer_)(*OverlappingY)), false);
1099  IFPACK_CHK_ERR(Inverse_->ApplyInverse(*tempX_,*tempY_));
1100  } else {
1101  IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX, *OverlappingY));
1102  }
1103 #else
1104  IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
1105 #endif
1106  }
1107  else {
1108  Epetra_MultiVector ReorderedX(*OverlappingX);
1109  Epetra_MultiVector ReorderedY(*OverlappingY);
1110  IFPACK_CHK_ERR(Reordering_->P(*OverlappingX,ReorderedX));
1111  IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
1112  IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,*OverlappingY));
1113  }
1114  }
1115 
1116  if (IsOverlapping()) {
1117  IFPACK_CHK_ERR(OverlappingMatrix_->ExportMultiVector(*OverlappingY,Y,
1118  CombineMode_));
1119  }
1120 
1121 #ifdef IFPACK_FLOPCOUNTERS
1122  // add flops. Note the we only have to add the newly counted
1123  // flops -- and each Inverse returns the cumulative sum
1124  double partial = Inverse_->ApplyInverseFlops();
1125  double total;
1126  Comm().SumAll(&partial, &total, 1);
1127  ApplyInverseFlops_ += total - pre_total;
1128 #endif
1129 
1130  // FIXME: right now I am skipping the overlap and singletons
1131  ++NumApplyInverse_;
1132  ApplyInverseTime_ += Time_->ElapsedTime();
1133 
1134  return(0);
1135 
1136 }
1137 
1138 //==============================================================================
1139 template<typename T>
1140 std::ostream& Ifpack_AdditiveSchwarz<T>::
1141 Print(std::ostream& os) const
1142 {
1143 #ifdef IFPACK_FLOPCOUNTERS
1144  double IF = InitializeFlops();
1145  double CF = ComputeFlops();
1146  double AF = ApplyInverseFlops();
1147 
1148  double IFT = 0.0, CFT = 0.0, AFT = 0.0;
1149  if (InitializeTime() != 0.0)
1150  IFT = IF / InitializeTime();
1151  if (ComputeTime() != 0.0)
1152  CFT = CF / ComputeTime();
1153  if (ApplyInverseTime() != 0.0)
1154  AFT = AF / ApplyInverseTime();
1155 #endif
1156 
1157  if (Matrix().Comm().MyPID())
1158  return(os);
1159 
1160  os << endl;
1161  os << "================================================================================" << endl;
1162  os << "Ifpack_AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
1163  if (CombineMode_ == Insert)
1164  os << "Combine mode = Insert" << endl;
1165  else if (CombineMode_ == Add)
1166  os << "Combine mode = Add" << endl;
1167  else if (CombineMode_ == Zero)
1168  os << "Combine mode = Zero" << endl;
1169  else if (CombineMode_ == Average)
1170  os << "Combine mode = Average" << endl;
1171  else if (CombineMode_ == AbsMax)
1172  os << "Combine mode = AbsMax" << endl;
1173 
1174  os << "Condition number estimate = " << Condest_ << endl;
1175  os << "Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
1176 
1177 #ifdef HAVE_IFPACK_PARALLEL_SUBDOMAIN_SOLVERS
1178  os << endl;
1179  os << "================================================================================" << endl;
1180  os << "Subcommunicator stats" << endl;
1181  os << "Number of MPI processes in simulation: " << NumMpiProcs_ << endl;
1182  os << "Number of subdomains: " << NumSubdomains_ << endl;
1183  os << "Number of MPI processes per subdomain: " << NumMpiProcsPerSubdomain_ << endl;
1184 #endif
1185 
1186  os << endl;
1187  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
1188  os << "----- ------- -------------- ------------ --------" << endl;
1189  os << "Initialize() " << std::setw(5) << NumInitialize()
1190  << " " << std::setw(15) << InitializeTime()
1191 #ifdef IFPACK_FLOPCOUNTERS
1192  << " " << std::setw(15) << 1.0e-6 * IF
1193  << " " << std::setw(15) << 1.0e-6 * IFT
1194 #endif
1195  << endl;
1196  os << "Compute() " << std::setw(5) << NumCompute()
1197  << " " << std::setw(15) << ComputeTime()
1198 #ifdef IFPACK_FLOPCOUNTERS
1199  << " " << std::setw(15) << 1.0e-6 * CF
1200  << " " << std::setw(15) << 1.0e-6 * CFT
1201 #endif
1202  << endl;
1203  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
1204  << " " << std::setw(15) << ApplyInverseTime()
1205 #ifdef IFPACK_FLOPCOUNTERS
1206  << " " << std::setw(15) << 1.0e-6 * AF
1207  << " " << std::setw(15) << 1.0e-6 * AFT
1208 #endif
1209  << endl;
1210  os << "================================================================================" << endl;
1211  os << endl;
1212 
1213  return(os);
1214 }
1215 
1216 #include "Ifpack_Condest.h"
1217 //==============================================================================
1218 template<typename T>
1220 Condest(const Ifpack_CondestType CT, const int MaxIters,
1221  const double Tol, Epetra_RowMatrix* Matrix_in)
1222 {
1223  if (!IsComputed()) // cannot compute right now
1224  return(-1.0);
1225 
1226  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
1227 
1228  return(Condest_);
1229 }
1230 
1231 #endif // IFPACK_ADDITIVESCHWARZ_H
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.
int OverlapLevel_
Level of overlap among the processors.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
string ReorderingType_
Type of reordering of the local matrix.
double ComputeFlops_
Contains the number of flops for Compute().
virtual int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Ifpack_AMDReordering: approximate minimum degree reordering.
bool UseTranspose_
If true, solve with the transpose (not supported by all solvers).
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
bool ComputeCondest_
If true, compute the condition number estimate each time Compute() is called.
double Condest_
Contains the estimated condition number.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
Ifpack_METISReordering: A class to reorder a graph using METIS.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
string Ifpack_toString(const int &x)
Convertes an integer to string.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
Teuchos::RefCountPtr< Ifpack_OverlappingRowMatrix > OverlappingMatrix_
Pointers to the overlapping matrix.
virtual int OverlapLevel() const
Returns the level of overlap.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
virtual bool IsOverlapping() const
Returns true is an overlapping matrix is present.
virtual double ComputeTime() const
Returns the time spent in Compute().
Teuchos::RefCountPtr< Ifpack_Reordering > Reordering_
Pointer to a reordering object.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Teuchos::RefCountPtr< Ifpack_ReorderFilter > ReorderedLocalizedMatrix_
Pointer to the reorderd matrix.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
bool IsOverlapping_
If true, overlapping is used.
virtual std::ostream & Print(std::ostream &) const
Prints major information about this preconditioner.
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's.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
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 HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
Ifpack_SingletonFilter: Filter based on matrix entries.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to X, returns the result in Y.
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.
string Label_
Contains the label of this object.
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 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 InitializeFlops() const
Returns the number of flops in the initialization phase.
bool IsComputed_
If true, the preconditioner has been successfully computed.
virtual const Teuchos::ParameterList & List() const
Returns a reference to the internally stored list.
virtual const char * Label() const
Returns a character string describing the operator.
Teuchos::RefCountPtr< T > Inverse_
Pointer to the local solver.
int NumCompute_
Contains the number of successful call to Compute().
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
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 double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
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)
virtual const Epetra_RowMatrix & Matrix() const
Returns a refernence to the internally stored matrix.
bool FilterSingletons_
Filter for singletons.
virtual double Condest() const
Returns the estimated condition number, or -1.0 if not computed.
virtual int NumCompute() const
Returns the number of calls to Compute().
Teuchos::RefCountPtr< Epetra_Time > Time_
Object used for timing purposes.
Teuchos::RefCountPtr< Ifpack_LocalFilter > LocalizedMatrix_
Localized version of Matrix_ or OverlappingMatrix_.