Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_AdditiveSchwarz_def.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 1
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef IFPACK2_ADDITIVESCHWARZ_DEF_HPP
00030 #define IFPACK2_ADDITIVESCHWARZ_DEF_HPP
00031 
00032 #include "Ifpack2_AdditiveSchwarz_decl.hpp"
00033 
00034 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
00035 #include "Xpetra_RowMatrix.hpp"
00036 #include "Xpetra_TpetraRowMatrix.hpp"
00037 #include "Zoltan2_XpetraRowMatrixInput.hpp"
00038 #include "Zoltan2_OrderingProblem.hpp"
00039 #endif
00040 
00041 #include "Ifpack2_Condest.hpp"
00042 #include "Ifpack2_OverlappingRowMatrix_def.hpp"
00043 #include "Ifpack2_LocalFilter_def.hpp"
00044 #include "Ifpack2_ReorderFilter_def.hpp"
00045 #include "Ifpack2_SingletonFilter_def.hpp"
00046 #ifdef HAVE_MPI
00047 #include "Teuchos_DefaultMpiComm.hpp"
00048 #endif
00049 
00050 namespace Ifpack2 {
00051 
00052 //==============================================================================
00053 template<class MatrixType,class LocalInverseType>
00054 AdditiveSchwarz<MatrixType,LocalInverseType>::AdditiveSchwarz(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > & Matrix_in, int OverlapLevel_in):
00055   Matrix_(Matrix_in),
00056   IsInitialized_(false),
00057   IsComputed_(false),
00058   IsOverlapping_(false),
00059   OverlapLevel_(OverlapLevel_in),
00060   CombineMode_(Tpetra::ADD),
00061   Condest_(-1.0),
00062   ComputeCondest_(true),
00063   UseReordering_(false),
00064   UseSubdomain_(false),
00065   FilterSingletons_(false),
00066   NumInitialize_(0),
00067   NumCompute_(0),
00068   NumApply_(0),
00069   InitializeTime_(0.0),
00070   ComputeTime_(0.0),
00071   ApplyTime_(0.0)
00072 {
00073   
00074   // Reset Overlap parameters
00075   if (Matrix_->getComm()->getSize() == 1)
00076     OverlapLevel_ = 0;
00077   
00078   if ((OverlapLevel_ != 0) && (Matrix_->getComm()->getSize() > 1))
00079     IsOverlapping_ = true;
00080 
00081   // Sets parameters to default values
00082   Teuchos::ParameterList List_in;
00083   setParameters(List_in);
00084 }
00085 
00086 //==============================================================================
00087 // Destructor
00088 template<class MatrixType,class LocalInverseType>
00089 AdditiveSchwarz<MatrixType,LocalInverseType>::~AdditiveSchwarz()
00090 {
00091 }
00092 
00093 //==============================================================================
00094 // Returns the Map associated with the domain of this operator, which must be compatible with X.getMap().
00095 template<class MatrixType,class LocalInverseType>
00096 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > > & AdditiveSchwarz<MatrixType,LocalInverseType>::getDomainMap() const 
00097 { 
00098   return Matrix_->getDomainMap();
00099 }
00100   
00101 //==============================================================================
00102 // Returns the Map associated with the range of this operator, which must be compatible with Y.getMap().
00103 template<class MatrixType,class LocalInverseType>
00104 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > & AdditiveSchwarz<MatrixType,LocalInverseType>::getRangeMap() const 
00105 {
00106   return Matrix_->getRangeMap();
00107 }
00108 
00109 //==============================================================================
00110 template<class MatrixType,class LocalInverseType>  
00111 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> > AdditiveSchwarz<MatrixType,LocalInverseType>::getMatrix() const
00112 {
00113   return Matrix_;
00114 }
00115   
00116 //==============================================================================
00117 // Applies the effect of the preconditione.
00118 template<class MatrixType,class LocalInverseType>
00119 void AdditiveSchwarz<MatrixType,LocalInverseType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
00120           Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
00121           Teuchos::ETransp mode,
00122           Scalar alpha,
00123           Scalar beta) const
00124 {
00125   typedef typename Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVectorType;
00126 
00127   TEUCHOS_TEST_FOR_EXCEPTION(IsComputed_ == false, std::runtime_error,
00128      "Ifpack2::AdditiveSchwarz::apply ERROR: isComputed() must be true prior to calling apply.");
00129 
00130   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00131      "Ifpack2::AdditiveSchwarz::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
00132 
00133   size_t NumVectors = X.getNumVectors();
00134 
00135   Time_->start();
00136 
00137   Teuchos::RCP<MultiVectorType> OverlappingX,OverlappingY,Xtmp;
00138   
00139   if(IsOverlapping_){
00140     // Setup if we're overlapping
00141     OverlappingX = Teuchos::rcp( new MultiVectorType(OverlappingMatrix_->getRowMap(), X.getNumVectors()) );
00142     OverlappingY = Teuchos::rcp( new MultiVectorType(OverlappingMatrix_->getRowMap(), X.getNumVectors()) );
00143     OverlappingY->putScalar(0.0);
00144     OverlappingX->putScalar(0.0);    
00145     OverlappingMatrix_->importMultiVector(X,*OverlappingX,Tpetra::INSERT);
00146     // FIXME from Ifpack1: Will not work with non-zero starting solutions.
00147   }
00148   else{
00149     Xtmp=Teuchos::rcp(new MultiVectorType(X));
00150     OverlappingX=Xtmp;
00151     OverlappingY=Teuchos::rcp(&Y,false);          
00152   }
00153 
00154   if (FilterSingletons_) {
00155     // process singleton filter
00156     MultiVectorType ReducedX(SingletonMatrix_->getRowMap(),NumVectors);
00157     MultiVectorType ReducedY(SingletonMatrix_->getRowMap(),NumVectors);
00158     SingletonMatrix_->SolveSingletons(*OverlappingX,*OverlappingY);
00159     SingletonMatrix_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX);
00160 
00161     // process reordering
00162     if (!UseReordering_) {
00163       Inverse_->apply(ReducedX,ReducedY);
00164     }
00165     else {
00166       MultiVectorType ReorderedX(ReducedX);
00167       MultiVectorType ReorderedY(ReducedY);
00168       ReorderedLocalizedMatrix_->permuteOriginalToReordered(ReducedX,ReorderedX);
00169       Inverse_->apply(ReorderedX,ReorderedY);
00170       ReorderedLocalizedMatrix_->permuteReorderedToOriginal(ReorderedY,ReducedY);
00171     }
00172 
00173     // finish up with singletons
00174     SingletonMatrix_->UpdateLHS(ReducedY,*OverlappingY);
00175   }
00176   else {
00177 
00178     // process reordering
00179     if (!UseReordering_) {
00180       Inverse_->apply(*OverlappingX,*OverlappingY);
00181     }
00182     else {
00183       MultiVectorType ReorderedX(*OverlappingX);
00184       MultiVectorType ReorderedY(*OverlappingY);
00185       ReorderedLocalizedMatrix_->permuteOriginalToReordered(*OverlappingX,ReorderedX);
00186       Inverse_->apply(ReorderedX,ReorderedY);
00187       ReorderedLocalizedMatrix_->permuteReorderedToOriginal(ReorderedY,*OverlappingY);
00188     }
00189   }
00190 
00191   if(IsOverlapping_)
00192     OverlappingMatrix_->exportMultiVector(*OverlappingY,Y,CombineMode_);
00193 
00194   NumApply_++;
00195   ApplyTime_ += Time_->stop();
00196 }
00197 
00198 //==============================================================================
00199 // Sets all parameters for the preconditioner.
00200 template<class MatrixType,class LocalInverseType>
00201 void AdditiveSchwarz<MatrixType,LocalInverseType>::setParameters(const Teuchos::ParameterList& List_in)
00202 {
00203   List_=List_in;
00204 
00205   // compute the condition number each time compute() is invoked.
00206   ComputeCondest_ = List_.get("schwarz: compute condest", ComputeCondest_);
00207   // combine mode
00208   if( Teuchos::ParameterEntry *combineModeEntry = List_.getEntryPtr("schwarz: combine mode") ) {
00209     if( typeid(std::string) == combineModeEntry->getAny().type() ) {
00210       std::string mode = List_.get("schwarz: combine mode", "Add");
00211       if (mode == "Add")
00212         CombineMode_ = Tpetra::ADD;
00213       else if (mode == "Insert")
00214         CombineMode_ = Tpetra::INSERT;
00215       else if (mode == "Replace")
00216         CombineMode_ = Tpetra::REPLACE;
00217       else if (mode == "AbsMax")
00218         CombineMode_ = Tpetra::ABSMAX;
00219       else {  
00220         TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error
00221            ,"Error, The Tpetra combine mode of \""<<mode<<"\" is not valid!  Only the values"
00222            " \"Add\", \"Insert\", \"Replace\", and \"AbsMax\" are accepted!"
00223            );
00224       }
00225     }
00226     else if ( typeid(Tpetra::CombineMode) == combineModeEntry->getAny().type() ) {
00227       CombineMode_ = Teuchos::any_cast<Tpetra::CombineMode>(combineModeEntry->getAny());
00228     }
00229     else {
00230       // Throw exception with good error message!
00231       Teuchos::getParameter<std::string>(List_,"schwarz: combine mode");
00232     }
00233   }
00234   else {
00235     // Make the default be a string to be consistent with the valid parameters!
00236     List_.get("schwarz: combine mode","Add");
00237   }
00238   
00239   // Will we be doing reordering?
00240   // Note: Unlike Ifpack we'll use a "schwarz: reordering list" to give to Zoltan2...
00241   UseReordering_ = List_.get("schwarz: use reordering",false);
00242 
00243 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
00244   // If we don't have Zoltan2, we just turn the reordering off completely...
00245   UseReordering_=false;
00246 #endif
00247 
00248   // Subdomain check
00249   if (List_.isParameter("schwarz: subdomain id") && List_.get("schwarz: subdomain id",-1) >  0)
00250     UseSubdomain_=true;
00251 
00252   // if true, filter singletons. NOTE: the filtered matrix can still have
00253   // singletons! A simple example: upper triangular matrix, if I remove
00254   // the lower node, I still get a matrix with a singleton! However, filter
00255   // singletons should help for PDE problems with Dirichlet BCs.
00256   FilterSingletons_ = List_.get("schwarz: filter singletons", FilterSingletons_);
00257 }
00258 
00259   
00260 //==============================================================================
00261 // Computes all (graph-related) data necessary to initialize the preconditioner.
00262 template<class MatrixType,class LocalInverseType>
00263 void AdditiveSchwarz<MatrixType,LocalInverseType>::initialize()
00264 {
00265   IsInitialized_ = false;
00266   IsComputed_ = false; // values required
00267   Condest_ = -1.0; // zero-out condest
00268 
00269   if (Time_ == Teuchos::null)
00270     Time_ = Teuchos::rcp( new Teuchos::Time("Ifpack2::AdditiveSchwarz"));
00271 
00272   Time_->start();
00273   
00274   // compute the overlapping matrix if necessary
00275   if (IsOverlapping_) {
00276     if(UseSubdomain_) {
00277       int sid = List_.get("subdomain id",-1);
00278       OverlappingMatrix_ = Teuchos::rcp( new Ifpack2::OverlappingRowMatrix<LocalMatrixType>(Matrix_, OverlapLevel_, sid) );
00279     }
00280     else
00281       OverlappingMatrix_ = Teuchos::rcp( new Ifpack2::OverlappingRowMatrix<LocalMatrixType>(Matrix_, OverlapLevel_) );
00282 
00283     TEUCHOS_TEST_FOR_EXCEPTION(OverlappingMatrix_ == Teuchos::null, std::runtime_error,
00284      "Ifpack2::AdditiveSchwarz::initialize ERROR: OverlappingRowMatrix constructor failed.");
00285   }
00286 
00287   // Setup
00288   setup();
00289 
00290   // Sanity Checks
00291   TEUCHOS_TEST_FOR_EXCEPTION(Inverse_ == Teuchos::null, std::runtime_error,
00292      "Ifpack2::AdditiveSchwarz::initialize ERROR: Setup() failed.");
00293  
00294   // Initialize inverse
00295   Inverse_->setParameters(List_);
00296   Inverse_->initialize();
00297 
00298   IsInitialized_ = true;
00299   NumInitialize_++;
00300   InitializeTime_ += Time_->stop();
00301 }
00302   
00303 //==============================================================================
00304 // Returns true if the  preconditioner has been successfully initialized, false otherwise.
00305 template<class MatrixType,class LocalInverseType>
00306 bool AdditiveSchwarz<MatrixType,LocalInverseType>::isInitialized() const
00307 {
00308   return IsInitialized_;
00309 }
00310   
00311 //==============================================================================
00312 // Computes all (coefficient) data necessary to apply the preconditioner.
00313 template<class MatrixType,class LocalInverseType>
00314 void AdditiveSchwarz<MatrixType,LocalInverseType>::compute()
00315 {
00316   if (!IsInitialized_) initialize();
00317 
00318   Time_->start();
00319   IsComputed_ = false;
00320   Condest_ = -1.0;
00321 
00322   Inverse_->compute();
00323   
00324   IsComputed_ = true; 
00325   ++NumCompute_;
00326   ComputeTime_ += Time_->stop();
00327 
00328 }
00329   
00330 //==============================================================================
00331 // Returns true if the  preconditioner has been successfully computed, false otherwise.
00332 template<class MatrixType,class LocalInverseType>
00333 bool AdditiveSchwarz<MatrixType,LocalInverseType>::isComputed() const
00334 {
00335   return IsComputed_;
00336 }
00337   
00338 //==============================================================================
00339 // Computes the condition number estimate and returns its value.
00340 template<class MatrixType,class LocalInverseType>
00341 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType 
00342 AdditiveSchwarz<MatrixType,LocalInverseType>::computeCondEst(CondestType CT,
00343                    LocalOrdinal MaxIters,
00344                    magnitudeType Tol,
00345                    const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &Matrix_in)
00346 {
00347   
00348   // If we haven't computed, we can't do a condest
00349   if (!isComputed()) return (-Teuchos::ScalarTraits<magnitudeType>::one());
00350 
00351   Condest_ = Ifpack2::Condest(*this, CT, MaxIters, Tol, Matrix_in);
00352   return(Condest_);
00353 }
00354 
00355 
00356 //==============================================================================
00357 // Returns the computed condition number estimate, or -1.0 if not computed.
00358 template<class MatrixType,class LocalInverseType>
00359 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType 
00360 AdditiveSchwarz<MatrixType,LocalInverseType>::getCondEst() const
00361 {
00362   return Condest_;
00363 }
00364   
00365 //==============================================================================
00366 // Returns the number of calls to initialize().
00367 template<class MatrixType,class LocalInverseType>  
00368 int AdditiveSchwarz<MatrixType,LocalInverseType>::getNumInitialize() const 
00369 {
00370   return NumInitialize_;
00371 }
00372   
00373 //==============================================================================
00374 // Returns the number of calls to compute().
00375 template<class MatrixType,class LocalInverseType>
00376 int AdditiveSchwarz<MatrixType,LocalInverseType>::getNumCompute() const
00377 {
00378   return NumCompute_;
00379 }
00380   
00381 //==============================================================================
00382 // Returns the number of calls to apply().
00383 template<class MatrixType,class LocalInverseType>
00384 int AdditiveSchwarz<MatrixType,LocalInverseType>::getNumApply() const 
00385 {
00386   return NumApply_;
00387 }
00388   
00389 //==============================================================================
00390 // Returns the time spent in initialize().
00391 template<class MatrixType,class LocalInverseType>
00392 double AdditiveSchwarz<MatrixType,LocalInverseType>::getInitializeTime() const
00393 {
00394   return InitializeTime_;
00395 }
00396   
00397 //==============================================================================
00398 // Returns the time spent in compute().
00399 template<class MatrixType,class LocalInverseType>
00400 double AdditiveSchwarz<MatrixType,LocalInverseType>::getComputeTime() const
00401 {
00402   return ComputeTime_;
00403 }
00404 
00405 //==============================================================================
00406 // Returns the time spent in Apply().
00407 template<class MatrixType,class LocalInverseType>
00408 double AdditiveSchwarz<MatrixType,LocalInverseType>::getApplyTime() const 
00409 {
00410   return ApplyTime_;
00411 }
00412 //==============================================================================
00413 /* Return a simple one-line description of this object. */
00414 template<class MatrixType,class LocalInverseType>
00415 std::string AdditiveSchwarz<MatrixType,LocalInverseType>::description() const
00416 {
00417   std::ostringstream oss;
00418   oss << Teuchos::Describable::description();
00419   if (isInitialized()) {
00420     if (isComputed()) {
00421       oss << "{status = initialized, computed";
00422     }
00423     else {
00424       oss << "{status = initialized, not computed";
00425     }
00426   }
00427   else {
00428     oss << "{status = not initialized, not computed";
00429   }
00430   oss<<"overlap level ="<<OverlapLevel_;
00431   oss << "}";
00432   return oss.str();
00433 }
00434 
00435 //==============================================================================
00436 /* Print the object with some verbosity level to an FancyOStream object. */
00437 template<class MatrixType,class LocalInverseType>
00438 void AdditiveSchwarz<MatrixType,LocalInverseType>::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const
00439 {
00440   using std::endl;
00441 
00442   os << endl;
00443   os << "================================================================================" << endl;
00444   os << "Ifpack2::AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
00445   if (CombineMode_ == Tpetra::INSERT)
00446     os << "Combine mode                          = Insert" << endl;
00447   else if (CombineMode_ == Tpetra::ADD)
00448     os << "Combine mode                          = Add" << endl;
00449   else if (CombineMode_ == Tpetra::REPLACE)
00450     os << "Combine mode                          = Replace" << endl;
00451   else if (CombineMode_ == Tpetra::ABSMAX)
00452     os << "Combine mode                          = AbsMax" << endl;
00453 
00454   os << "Condition number estimate             = " << Condest_ << endl;
00455   os << "Global number of rows                 = " << Matrix_->getGlobalNumRows() << endl;
00456 
00457   os << endl;
00458   os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00459   os << "-----           -------   --------------       ------------     --------" << endl;
00460   os << "Initialize()    "   << std::setw(5) << getNumInitialize()
00461      << "  " << std::setw(15) << getInitializeTime() << endl;
00462   os << "Compute()       "   << std::setw(5) << getNumCompute() 
00463      << "  " << std::setw(15) << getComputeTime() << endl;
00464   os << "Apply()         "   << std::setw(5) << getNumApply() 
00465      << "  " << std::setw(15) << getApplyTime() <<endl;
00466   os << "================================================================================" << endl;
00467   os << endl;
00468 }
00469 
00470 //==============================================================================
00471 // Prints basic information on iostream. This function is used by operator<<.
00472 template<class MatrixType,class LocalInverseType>
00473 std::ostream& AdditiveSchwarz<MatrixType,LocalInverseType>::print(std::ostream& os) const
00474 {
00475   Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
00476   fos.setOutputToRootOnly(0);
00477   describe(fos);
00478   return(os);
00479 }
00480 
00481 //==============================================================================
00482 // Returns the level of overlap.
00483 template<class MatrixType,class LocalInverseType>
00484 int AdditiveSchwarz<MatrixType,LocalInverseType>::getOverlapLevel() const
00485 {
00486   return OverlapLevel_;
00487 }
00488 
00489 //==============================================================================
00490 // Sets up the localized matrix and the singleton filter.
00491 template<class MatrixType,class LocalInverseType>
00492 void AdditiveSchwarz<MatrixType,LocalInverseType>::setup() 
00493 {  
00494 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
00495   typedef typename Xpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpetraMatrixType;
00496   typedef typename Xpetra::TpetraRowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpetraTpetraMatrixType;
00497 #endif
00498 
00499   Teuchos::RCP<LocalMatrixType> ActiveMatrix;
00500 
00501   // Create Localized Matrix
00502   if (OverlappingMatrix_ != Teuchos::null) {
00503     if (UseSubdomain_) {
00504       //      int sid = List_.get("subdomain id",-1);
00505       throw std::runtime_error("Ifpack2::AdditiveSchwarz subdomain code not yet supported.");
00506       //Ifpack2_NodeFilter *tt = new Ifpack2_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
00507       //LocalizedMatrix_ = Teuchos::rcp(tt);
00508     }
00509     else
00510       LocalizedMatrix_ = Teuchos::rcp( new Ifpack2::LocalFilter<LocalMatrixType>(OverlappingMatrix_) );
00511   }
00512   else {
00513     if (UseSubdomain_) {
00514       //      int sid = List_.get("subdomain id",-1);
00515       throw std::runtime_error("Ifpack2::AdditiveSchwarz subdomain code not yet supported.");
00516     }
00517     else
00518       LocalizedMatrix_ = Teuchos::rcp( new Ifpack2::LocalFilter<LocalMatrixType>(Matrix_) );
00519   }
00520 
00521   // Mark localized matrix as active
00522   ActiveMatrix = LocalizedMatrix_;
00523   TEUCHOS_TEST_FOR_EXCEPTION(LocalizedMatrix_ == Teuchos::null, std::runtime_error,
00524            "Ifpack2::AdditiveSchwarz::Setup() ERROR: LocalFilter failed.");
00525 
00526   // Singleton Filtering
00527   if (FilterSingletons_) {
00528     SingletonMatrix_ = Teuchos::rcp( new Ifpack2::SingletonFilter<LocalMatrixType>(LocalizedMatrix_) );
00529     ActiveMatrix = SingletonMatrix_;
00530   }
00531 
00532 
00533   // Do reordering
00534   if (UseReordering_) {
00535 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
00536     // Unlike Ifpack, Zoltan2 does all the dirty work here.
00537     Teuchos::ParameterList zlist = List_.sublist("schwarz: reordering list");
00538     XpetraTpetraMatrixType XpetraMatrix(ActiveMatrix);
00539     Zoltan2::XpetraRowMatrixInput<XpetraMatrixType> Zoltan2Matrix(Teuchos::rcp<XpetraMatrixType>(&XpetraMatrix,false));
00540 
00541 #ifdef HAVE_MPI
00542     // Grab the MPI Communicator and build the ordering problem with that
00543     MPI_Comm mycomm;
00544 
00545     Teuchos::RCP<const Teuchos::MpiComm<int> > mpicomm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(ActiveMatrix->getComm());
00546     if(mpicomm == Teuchos::null) mycomm = MPI_COMM_SELF;
00547     else mycomm = *(mpicomm->getRawMpiComm());
00548     Zoltan2::OrderingProblem<Zoltan2::XpetraRowMatrixInput<XpetraMatrixType> > MyOrderingProblem(&Zoltan2Matrix, &zlist,mycomm);    
00549 #else
00550     Zoltan2::OrderingProblem<Zoltan2::XpetraRowMatrixInput<XpetraMatrixType> > MyOrderingProblem(&Zoltan2Matrix, &zlist);    
00551 #endif
00552     MyOrderingProblem.solve();
00553 
00554     // Now create the reordered matrix & mark it as active
00555     ReorderedLocalizedMatrix_ =  Teuchos::rcp(new Ifpack2::ReorderFilter<LocalMatrixType>(ActiveMatrix,Teuchos::rcp( new Zoltan2::OrderingSolution<GlobalOrdinal,LocalOrdinal>(*MyOrderingProblem.getSolution()))));
00556 
00557     ActiveMatrix = ReorderedLocalizedMatrix_;
00558 #else
00559     throw std::runtime_error("Ifpack2::AdditiveSchwarz::Setup() You need to compile with Zoltan2 to support reordering.");
00560 #endif
00561   }
00562 
00563   // Build the inverse
00564   Inverse_ = Teuchos::rcp(new LocalInverseType(ActiveMatrix));
00565   TEUCHOS_TEST_FOR_EXCEPTION(Inverse_ == Teuchos::null, std::runtime_error,
00566            "Ifpack2::AdditiveSchwarz::Setup() ERROR: Inverse constructor failed.");
00567 }
00568 
00569 
00570 }// namespace Ifpack2
00571 
00572 #endif // IFPACK2_ADDITIVESCHWARZ_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends