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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 */
00042 
00043 #ifndef IFPACK2_ADDITIVESCHWARZ_DEF_HPP
00044 #define IFPACK2_ADDITIVESCHWARZ_DEF_HPP
00045 
00046 #include "Ifpack2_AdditiveSchwarz_decl.hpp"
00047 
00048 // AdditiveSchwarz uses OneLevelFactory to create a default inner
00049 // preconditioner.
00050 //
00051 // FIXME (mfh 13 Dec 2013) For some inexplicable reason, I have to
00052 // include the _decl and _def headers separately here; including just
00053 // Ifpack2_Details_OneLevelFactory.hpp doesn't work.  It probably has
00054 // something to do with ETI, but I don't fully understand what.
00055 #include "Ifpack2_Details_OneLevelFactory_decl.hpp"
00056 #include "Ifpack2_Details_OneLevelFactory_def.hpp"
00057 
00058 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
00059 #include "Xpetra_RowMatrix.hpp"
00060 #include "Xpetra_TpetraRowMatrix.hpp"
00061 #include "Zoltan2_XpetraRowMatrixAdapter.hpp"
00062 #include "Zoltan2_OrderingProblem.hpp"
00063 #include "Zoltan2_OrderingSolution.hpp"
00064 #endif
00065 
00066 #if defined(HAVE_IFPACK2_EXPERIMENTAL) && defined(HAVE_IFPACK2_SUPPORTGRAPH)
00067 #include "Ifpack2_SupportGraph_decl.hpp"
00068 #endif
00069 
00070 #include "Ifpack2_Condest.hpp"
00071 #include "Ifpack2_Details_CanChangeMatrix.hpp"
00072 #include "Ifpack2_LocalFilter_def.hpp"
00073 #include "Ifpack2_OverlappingRowMatrix_def.hpp"
00074 #include "Ifpack2_Parameters.hpp"
00075 #include "Ifpack2_ReorderFilter_def.hpp"
00076 #include "Ifpack2_SingletonFilter_def.hpp"
00077 
00078 #ifdef HAVE_MPI
00079 #include "Teuchos_DefaultMpiComm.hpp"
00080 #endif
00081 
00082 #include <locale> // std::toupper
00083 
00084 namespace Ifpack2 {
00085 
00086 namespace Details {
00087 
00103 template<class PrecType>
00104 class OneLevelPreconditionerNamer {
00105 public:
00107   static std::string name () {
00108     // The default implementation returns an invalid preconditioner
00109     // name.  This ensures that AdditiveSchwarz won't try to create a
00110     // preconditioner it doesn't know how to create.  This is better
00111     // than providing a valid default that is a different class than
00112     // the user expects.
00113     return "INVALID";
00114   }
00115 };
00116 
00117 //
00118 // Partial specialization for Ifpack2::Preconditioner.
00119 // It picks a reasonable default subdomain solver.
00120 //
00121 
00122 template<class S, class LO, class GO, class NT>
00123 class OneLevelPreconditionerNamer< ::Ifpack2::Preconditioner<S, LO, GO, NT> > {
00124 public:
00125   static std::string name () {
00126     // The default inner preconditioner is "ILUT", for backwards
00127     // compatibility with the original AdditiveSchwarz implementation.
00128     return "ILUT";
00129   }
00130 };
00131 
00132 //
00133 // Partial specializations for each single-level preconditioner.
00134 //
00135 
00136 template<class MatrixType>
00137 class OneLevelPreconditionerNamer< ::Ifpack2::Chebyshev<MatrixType> > {
00138 public:
00139   static std::string name () {
00140     return "CHEBYSHEV";
00141   }
00142 };
00143 
00144 template<class MatrixType>
00145 class OneLevelPreconditionerNamer< ::Ifpack2::Details::DenseSolver<MatrixType> > {
00146 public:
00147   static std::string name () {
00148     return "DENSE";
00149   }
00150 };
00151 
00152 #ifdef HAVE_IFPACK2_AMESOS2
00153 template<class MatrixType>
00154 class OneLevelPreconditionerNamer< ::Ifpack2::Details::Amesos2Wrapper<MatrixType> > {
00155 public:
00156   static std::string name () {
00157     return "AMESOS2";
00158   }
00159 };
00160 #endif // HAVE_IFPACK2_AMESOS2
00161 
00162 template<class MatrixType>
00163 class OneLevelPreconditionerNamer< ::Ifpack2::Diagonal<MatrixType> > {
00164 public:
00165   static std::string name () {
00166     return "DIAGONAL";
00167   }
00168 };
00169 
00170 template<class MatrixType>
00171 class OneLevelPreconditionerNamer< ::Ifpack2::ILUT<MatrixType> > {
00172 public:
00173   static std::string name () {
00174     return "ILUT";
00175   }
00176 };
00177 
00178 template<class MatrixType>
00179 class OneLevelPreconditionerNamer< ::Ifpack2::Relaxation<MatrixType> > {
00180 public:
00181   static std::string name () {
00182     return "RELAXATION";
00183   }
00184 };
00185 
00186 template<class MatrixType>
00187 class OneLevelPreconditionerNamer< ::Ifpack2::RILUK<MatrixType> > {
00188 public:
00189   static std::string name () {
00190     return "RILUK";
00191   }
00192 };
00193 
00194 template<class MatrixType>
00195 class OneLevelPreconditionerNamer< ::Ifpack2::Krylov<MatrixType> > {
00196 public:
00197   static std::string name () {
00198     return "KRYLOV";
00199   }
00200 };
00201 
00202 template<class MatrixType>
00203 class OneLevelPreconditionerNamer< ::Ifpack2::IdentitySolver<MatrixType> > {
00204 public:
00205   static std::string name () {
00206     return "IDENTITY";
00207   }
00208 };
00209 
00210 } // namespace Details
00211 
00212 
00213 template<class MatrixType, class LocalInverseType>
00214 bool
00215 AdditiveSchwarz<MatrixType, LocalInverseType>::hasInnerPrecName () const
00216 {
00217   const char* options[4] = {
00218     "inner preconditioner name",
00219     "subdomain solver name",
00220     "schwarz: inner preconditioner name",
00221     "schwarz: subdomain solver name"
00222   };
00223   const int numOptions = 4;
00224   bool match = false;
00225   for (int k = 0; k < numOptions && ! match; ++k) {
00226     if (List_.isParameter (options[k])) {
00227       return true;
00228     }
00229   }
00230   return false;
00231 }
00232 
00233 
00234 template<class MatrixType, class LocalInverseType>
00235 void
00236 AdditiveSchwarz<MatrixType, LocalInverseType>::removeInnerPrecName ()
00237 {
00238   const char* options[4] = {
00239     "inner preconditioner name",
00240     "subdomain solver name",
00241     "schwarz: inner preconditioner name",
00242     "schwarz: subdomain solver name"
00243   };
00244   const int numOptions = 4;
00245   for (int k = 0; k < numOptions; ++k) {
00246     List_.remove (options[k], false);
00247   }
00248 }
00249 
00250 
00251 
00252 
00253 
00254 template<class MatrixType, class LocalInverseType>
00255 std::string
00256 AdditiveSchwarz<MatrixType, LocalInverseType>::innerPrecName () const
00257 {
00258   const char* options[4] = {
00259     "inner preconditioner name",
00260     "subdomain solver name",
00261     "schwarz: inner preconditioner name",
00262     "schwarz: subdomain solver name"
00263   };
00264   const int numOptions = 4;
00265   std::string newName;
00266   bool match = false;
00267 
00268   // As soon as one parameter option matches, ignore all others.
00269   for (int k = 0; k < numOptions && ! match; ++k) {
00270     if (List_.isParameter (options[k])) {
00271       // try-catch block protects against incorrect type errors.
00272       //
00273       // FIXME (mfh 04 Jan 2013) We should instead catch and report
00274       // type errors.
00275       try {
00276         newName = List_.get<std::string> (options[k]);
00277         match = true;
00278       } catch (...) {}
00279     }
00280   }
00281   return match ? newName : defaultInnerPrecName ();
00282 }
00283 
00284 
00285 template<class MatrixType, class LocalInverseType>
00286 void
00287 AdditiveSchwarz<MatrixType, LocalInverseType>::removeInnerPrecParams ()
00288 {
00289   const char* options[4] = {
00290     "inner preconditioner parameters",
00291     "subdomain solver parameters",
00292     "schwarz: inner preconditioner parameters",
00293     "schwarz: subdomain solver parameters"
00294   };
00295   const int numOptions = 4;
00296 
00297   // As soon as one parameter option matches, ignore all others.
00298   for (int k = 0; k < numOptions; ++k) {
00299     List_.remove (options[k], false);
00300   }
00301 }
00302 
00303 
00304 template<class MatrixType, class LocalInverseType>
00305 std::pair<Teuchos::ParameterList, bool>
00306 AdditiveSchwarz<MatrixType, LocalInverseType>::innerPrecParams () const
00307 {
00308   const char* options[4] = {
00309     "inner preconditioner parameters",
00310     "subdomain solver parameters",
00311     "schwarz: inner preconditioner parameters",
00312     "schwarz: subdomain solver parameters"
00313   };
00314   const int numOptions = 4;
00315   Teuchos::ParameterList params;
00316 
00317   // As soon as one parameter option matches, ignore all others.
00318   bool match = false;
00319   for (int k = 0; k < numOptions && ! match; ++k) {
00320     if (List_.isSublist (options[k])) {
00321       params = List_.sublist (options[k]);
00322       match = true;
00323     }
00324   }
00325   // Default is an empty list of parameters.
00326   return std::make_pair (params, match);
00327 }
00328 
00329 
00330 template<class MatrixType, class LocalInverseType>
00331 std::string
00332 AdditiveSchwarz<MatrixType, LocalInverseType>::defaultInnerPrecName ()
00333 {
00334   // FIXME (mfh 14 Dec 2013) We want to get rid of the
00335   // LocalInverseType template parameter.  Soon, we will add an "inner
00336   // preconditioner" string parameter to the input ParameterList.  For
00337   // now, we map statically from LocalInverseType to its string name,
00338   // and use the string name to create the inner preconditioner.
00339   return Details::OneLevelPreconditionerNamer<LocalInverseType>::name ();
00340 }
00341 
00342 
00343 template<class MatrixType, class LocalInverseType>
00344 AdditiveSchwarz<MatrixType, LocalInverseType>::
00345 AdditiveSchwarz (const Teuchos::RCP<const row_matrix_type>& A) :
00346   Matrix_ (A),
00347   IsInitialized_ (false),
00348   IsComputed_ (false),
00349   IsOverlapping_ (false),
00350   OverlapLevel_ (0),
00351   CombineMode_ (Tpetra::ADD),
00352   Condest_ (-Teuchos::ScalarTraits<magnitude_type>::one ()),
00353   ComputeCondest_ (true),
00354   UseReordering_ (false),
00355   ReorderingAlgorithm_ ("none"),
00356   UseSubdomain_ (false),
00357   FilterSingletons_ (false),
00358   NumInitialize_ (0),
00359   NumCompute_ (0),
00360   NumApply_ (0),
00361   InitializeTime_ (0.0),
00362   ComputeTime_ (0.0),
00363   ApplyTime_ (0.0)
00364 {
00365   Teuchos::ParameterList plist;
00366   setParameters (plist); // Set parameters to default values
00367 }
00368 
00369 template<class MatrixType, class LocalInverseType>
00370 AdditiveSchwarz<MatrixType, LocalInverseType>::
00371 AdditiveSchwarz (const Teuchos::RCP<const row_matrix_type>& A,
00372                  const int overlapLevel) :
00373   Matrix_ (A),
00374   IsInitialized_ (false),
00375   IsComputed_ (false),
00376   IsOverlapping_ (false),
00377   OverlapLevel_ (overlapLevel),
00378   CombineMode_ (Tpetra::ADD),
00379   Condest_ (-Teuchos::ScalarTraits<magnitude_type>::one ()),
00380   ComputeCondest_ (true),
00381   UseReordering_ (false),
00382   ReorderingAlgorithm_ ("none"),
00383   UseSubdomain_ (false),
00384   FilterSingletons_ (false),
00385   NumInitialize_ (0),
00386   NumCompute_ (0),
00387   NumApply_ (0),
00388   InitializeTime_ (0.0),
00389   ComputeTime_ (0.0),
00390   ApplyTime_ (0.0)
00391 {
00392   Teuchos::ParameterList plist;
00393   setParameters (plist); // Set parameters to default values
00394 }
00395 
00396 
00397 template<class MatrixType,class LocalInverseType>
00398 AdditiveSchwarz<MatrixType,LocalInverseType>::~AdditiveSchwarz () {}
00399 
00400 
00401 template<class MatrixType,class LocalInverseType>
00402 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type > >
00403 AdditiveSchwarz<MatrixType,LocalInverseType>::getDomainMap() const
00404 {
00405   TEUCHOS_TEST_FOR_EXCEPTION(
00406     Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
00407     "getDomainMap: The matrix to precondition is null.  You must either pass "
00408     "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
00409     "input, before you may call this method.");
00410   return Matrix_->getDomainMap ();
00411 }
00412 
00413 
00414 template<class MatrixType,class LocalInverseType>
00415 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
00416 AdditiveSchwarz<MatrixType,LocalInverseType>::getRangeMap () const
00417 {
00418   TEUCHOS_TEST_FOR_EXCEPTION(
00419     Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
00420     "getRangeMap: The matrix to precondition is null.  You must either pass "
00421     "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
00422     "input, before you may call this method.");
00423   return Matrix_->getRangeMap ();
00424 }
00425 
00426 
00427 template<class MatrixType,class LocalInverseType>
00428 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
00429 {
00430   return Matrix_;
00431 }
00432 
00433 
00434 template<class MatrixType,class LocalInverseType>
00435 void
00436 AdditiveSchwarz<MatrixType,LocalInverseType>::
00437 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
00438        Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
00439        Teuchos::ETransp mode,
00440        scalar_type alpha,
00441        scalar_type beta) const
00442 {
00443   using Teuchos::Time;
00444   using Teuchos::TimeMonitor;
00445   using Teuchos::RCP;
00446   using Teuchos::rcp;
00447 
00448   const std::string timerName ("Ifpack2::AdditiveSchwarz::apply");
00449   RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
00450   if (timer.is_null ()) {
00451     timer = TimeMonitor::getNewCounter (timerName);
00452   }
00453 
00454   { // Start timing here.
00455     TimeMonitor timeMon (*timer);
00456 
00457     const scalar_type ZERO = Teuchos::ScalarTraits<scalar_type>::zero ();
00458     const scalar_type ONE = Teuchos::ScalarTraits<scalar_type>::one ();
00459 
00460     TEUCHOS_TEST_FOR_EXCEPTION(
00461       ! IsComputed_, std::runtime_error,
00462       "Ifpack2::AdditiveSchwarz::apply: "
00463       "isComputed() must be true before you may call apply().");
00464     TEUCHOS_TEST_FOR_EXCEPTION(
00465       Matrix_.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::apply: "
00466       "The input matrix A is null, but the preconditioner says that it has "
00467       "been computed (isComputed() is true).  This should never happen, since "
00468       "setMatrix() should always mark the preconditioner as not computed if "
00469       "its argument is null.  "
00470       "Please report this bug to the Ifpack2 developers.");
00471     TEUCHOS_TEST_FOR_EXCEPTION(
00472       Inverse_.is_null (), std::runtime_error,
00473       "Ifpack2::AdditiveSchwarz::apply: The subdomain solver is null.  "
00474       "This can only happen if you called setInnerPreconditioner() with a null "
00475       "input, after calling initialize() or compute().  If you choose to call "
00476       "setInnerPreconditioner() with a null input, you must then call it with "
00477       "a nonnull input before you may call initialize() or compute().");
00478     TEUCHOS_TEST_FOR_EXCEPTION(
00479       X.getNumVectors() != Y.getNumVectors(), std::invalid_argument,
00480       "Ifpack2::AdditiveSchwarz::apply: "
00481       "X and Y must have the same number of columns.  X has "
00482       << X.getNumVectors() << " columns, but Y has " << Y.getNumVectors() << ".");
00483     TEUCHOS_TEST_FOR_EXCEPTION(
00484       alpha != ONE, std::logic_error,
00485       "Ifpack2::AdditiveSchwarz::apply: Not implemented for alpha != 1.");
00486     TEUCHOS_TEST_FOR_EXCEPTION(
00487       beta != ZERO, std::logic_error,
00488       "Ifpack2::AdditiveSchwarz::apply: Not implemented for beta != 0.");
00489 
00490     const size_t numVectors = X.getNumVectors ();
00491 
00492     RCP<MV> OverlappingX,OverlappingY;
00493 
00494     if (IsOverlapping_) {
00495       TEUCHOS_TEST_FOR_EXCEPTION(
00496         OverlappingMatrix_.is_null (), std::logic_error,
00497         "Ifpack2::AdditiveSchwarz::apply: The overlapping matrix is null.  "
00498         "This should never happen if IsOverlapping_ is true.  "
00499         "Please report this bug to the Ifpack2 developers.");
00500 
00501       // Setup if we're overlapping
00502       //
00503       // MV's constructor fills with zeros.
00504       OverlappingX = rcp (new MV (OverlappingMatrix_->getRowMap (), numVectors));
00505       OverlappingY = rcp (new MV (OverlappingMatrix_->getRowMap (), numVectors));
00506       OverlappingMatrix_->importMultiVector (X, *OverlappingX, Tpetra::INSERT);
00507       // FIXME from Ifpack1: Will not work with non-zero starting solutions.
00508     }
00509     else {
00510       TEUCHOS_TEST_FOR_EXCEPTION(
00511         localMap_.is_null (), std::logic_error,
00512         "Ifpack2::AdditiveSchwarz::apply: localMap_ is null.");
00513 
00514       // MV's constructor fills with zeros.
00515       //
00516       // localMap_ has the same number of indices on each process that
00517       // Matrix_->getRowMap() does on that process.  Thus, we can do
00518       // the Import step without creating a new MV, just by viewing
00519       // OverlappingX using Matrix_->getRowMap ().
00520       OverlappingX = rcp (new MV (localMap_, numVectors));
00521       OverlappingY = rcp (new MV (localMap_, numVectors));
00522 
00523       RCP<MV> globalOverlappingX =
00524         OverlappingX->offsetViewNonConst (Matrix_->getRowMap (), 0);
00525 
00526       // Create Import object on demand, if necessary.
00527       if (DistributedImporter_.is_null ()) {
00528         // FIXME (mfh 15 Apr 2014) Why can't we just ask the Matrix
00529         // for its Import object?  Of course a general RowMatrix might
00530         // not necessarily have one.
00531         DistributedImporter_ =
00532           rcp (new import_type (Matrix_->getRowMap (), Matrix_->getDomainMap ()));
00533       }
00534       globalOverlappingX->doImport (X, *DistributedImporter_, Tpetra::INSERT);
00535     }
00536 
00537     if (FilterSingletons_) {
00538       // process singleton filter
00539       MV ReducedX (SingletonMatrix_->getRowMap (), numVectors);
00540       MV ReducedY (SingletonMatrix_->getRowMap (), numVectors);
00541       SingletonMatrix_->SolveSingletons (*OverlappingX, *OverlappingY);
00542       SingletonMatrix_->CreateReducedRHS (*OverlappingY, *OverlappingX, ReducedX);
00543 
00544       // process reordering
00545       if (! UseReordering_) {
00546         Inverse_->apply (ReducedX, ReducedY);
00547       }
00548       else {
00549         MV ReorderedX (ReducedX);
00550         MV ReorderedY (ReducedY);
00551         ReorderedLocalizedMatrix_->permuteOriginalToReordered (ReducedX, ReorderedX);
00552         Inverse_->apply (ReorderedX, ReorderedY);
00553         ReorderedLocalizedMatrix_->permuteReorderedToOriginal (ReorderedY, ReducedY);
00554       }
00555 
00556       // finish up with singletons
00557       SingletonMatrix_->UpdateLHS (ReducedY, *OverlappingY);
00558     }
00559     else {
00560 
00561       // process reordering
00562       if (! UseReordering_) {
00563         Inverse_->apply (*OverlappingX, *OverlappingY);
00564       }
00565       else {
00566         MV ReorderedX = createCopy(*OverlappingX);
00567         MV ReorderedY = createCopy(*OverlappingY);
00568         ReorderedLocalizedMatrix_->permuteOriginalToReordered (*OverlappingX, ReorderedX);
00569         Inverse_->apply (ReorderedX, ReorderedY);
00570         ReorderedLocalizedMatrix_->permuteReorderedToOriginal (ReorderedY, *OverlappingY);
00571       }
00572     }
00573 
00574     if (IsOverlapping_) {
00575       OverlappingMatrix_->exportMultiVector (*OverlappingY, Y, CombineMode_);
00576     }
00577     else {
00578       // mfh 16 Apr 2014: Make a view of Y with the same Map as
00579       // OverlappingY, so that we can copy OverlappingY into Y.  This
00580       // replaces code that iterates over all entries of OverlappingY,
00581       // copying them one at a time into Y.  That code assumed that
00582       // the rows of Y and the rows of OverlappingY have the same
00583       // global indices in the same order; see Bug 5992.
00584       RCP<MV> Y_view = Y.offsetViewNonConst (OverlappingY->getMap (), 0);
00585       Tpetra::deep_copy (*Y_view, *OverlappingY);
00586     }
00587   } // Stop timing here.
00588 
00589   ++NumApply_;
00590 
00591   // timer->totalElapsedTime() returns the total time over all timer
00592   // calls.  Thus, we use = instead of +=.
00593   ApplyTime_ = timer->totalElapsedTime ();
00594 }
00595 
00596 
00597 template<class MatrixType,class LocalInverseType>
00598 void AdditiveSchwarz<MatrixType,LocalInverseType>::
00599 setParameters (const Teuchos::ParameterList& plist)
00600 {
00601   // mfh 18 Nov 2013: Ifpack2's setParameters() method passes in the
00602   // input list as const.  This means that we have to copy it before
00603   // validation or passing into setParameterList().
00604   List_ = plist;
00605   this->setParameterList (Teuchos::rcpFromRef (List_));
00606 }
00607 
00608 
00609 
00610 template<class MatrixType,class LocalInverseType>
00611 void AdditiveSchwarz<MatrixType,LocalInverseType>::
00612 setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00613 {
00614   using Tpetra::CombineMode;
00615   using Teuchos::getIntegralValue;
00616   using Teuchos::ParameterEntry;
00617   using Teuchos::ParameterEntryValidator;
00618   using Teuchos::RCP;
00619   using Teuchos::rcp_dynamic_cast;
00620   using Teuchos::StringToIntegralParameterEntryValidator;
00621 
00622   if (plist.is_null ()) {
00623     // Assume that the user meant to set default parameters by passing
00624     // in an empty list.
00625     this->setParameterList (Teuchos::parameterList ());
00626   }
00627   // At this point, plist should be nonnull.
00628   TEUCHOS_TEST_FOR_EXCEPTION(
00629     plist.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
00630     "setParameterList: plist is null.  This should never happen, since the "
00631     "method should have replaced a null input list with a nonnull empty list "
00632     "by this point.  Please report this bug to the Ifpack2 developers.");
00633 
00634   // try {
00635   //   List_.validateParameters (* getValidParameters ());
00636   // }
00637   // catch (std::exception& e) {
00638   //   std::cerr << "Ifpack2::AdditiveSchwarz::setParameterList: Validation failed with the following error message: " << e.what () << std::endl;
00639   //   throw e;
00640   // }
00641 
00642   // mfh 18 Nov 2013: Supplying the current value as the default value
00643   // when calling ParameterList::get() ensures "delta" behavior when
00644   // users pass in new parameters: any unspecified parameters in the
00645   // new list retain their values in the old list.  This preserves
00646   // backwards compatiblity with this class' previous behavior.  Note
00647   // that validateParametersAndSetDefaults() would have different
00648   // behavior: any parameters not in the new list would get default
00649   // values, which could be different than their values in the
00650   // original list.
00651 
00652   ComputeCondest_ = plist->get ("schwarz: compute condest", ComputeCondest_);
00653 
00654   bool gotCombineMode = false;
00655   try {
00656     CombineMode_ = getIntegralValue<Tpetra::CombineMode> (List_, "schwarz: combine mode");
00657     gotCombineMode = true;
00658   }
00659   catch (Teuchos::Exceptions::InvalidParameterName&) {
00660     // The caller didn't provide that parameter.  Just keep the
00661     // existing value of CombineMode_.
00662     gotCombineMode = true;
00663   }
00664   catch (Teuchos::Exceptions::InvalidParameterType&) {
00665     // The user perhaps supplied it as an Tpetra::CombineMode enum
00666     // value.  Let's try again (below).  If it doesn't succeed, we
00667     // know that the type is wrong, so we can let it throw whatever
00668     // exception it would throw.
00669   }
00670   // Try to get the combine mode as an integer.
00671   if (! gotCombineMode) {
00672     try {
00673       CombineMode_ = plist->get ("schwarz: combine mode", CombineMode_);
00674       gotCombineMode = true;
00675     }
00676     catch (Teuchos::Exceptions::InvalidParameterType&) {}
00677   }
00678   // Try to get the combine mode as a string.  If this works, use the
00679   // validator to convert to int.  This is painful, but necessary in
00680   // order to do validation, since the input list doesn't come with a
00681   // validator.
00682   if (! gotCombineMode) {
00683     const ParameterEntry& validEntry =
00684       getValidParameters ()->getEntry ("schwarz: combine mode");
00685     RCP<const ParameterEntryValidator> v = validEntry.validator ();
00686     typedef StringToIntegralParameterEntryValidator<CombineMode> vs2e_type;
00687     RCP<const vs2e_type> vs2e = rcp_dynamic_cast<const vs2e_type> (v, true);
00688 
00689     const ParameterEntry& inputEntry = plist->getEntry ("schwarz: combine mode");
00690     CombineMode_ = vs2e->getIntegralValue (inputEntry, "schwarz: combine mode");
00691     gotCombineMode = true;
00692   }
00693   (void) gotCombineMode; // forestall "set but not used" compiler warning
00694 
00695   OverlapLevel_ = plist->get ("schwarz: overlap level", OverlapLevel_);
00696 
00697   // We set IsOverlapping_ in initialize(), once we know that Matrix_ is nonnull.
00698 
00699   // Will we be doing reordering?  Unlike Ifpack, we'll use a
00700   // "schwarz: reordering list" to give to Zoltan2.
00701   UseReordering_ = plist->get ("schwarz: use reordering", UseReordering_);
00702 
00703 #if !defined(HAVE_IFPACK2_XPETRA) || !defined(HAVE_IFPACK2_ZOLTAN2)
00704   TEUCHOS_TEST_FOR_EXCEPTION(
00705     UseReordering_, std::invalid_argument, "Ifpack2::AdditiveSchwarz::"
00706     "setParameters: You specified \"schwarz: use reordering\" = true.  "
00707     "This is only valid when Trilinos was built with Ifpack2, Xpetra, and "
00708     "Zoltan2 enabled.  Either Xpetra or Zoltan2 was not enabled in your build "
00709     "of Trilinos.");
00710 #endif
00711 
00712   // FIXME (mfh 18 Nov 2013) Now would be a good time to validate the
00713   // "schwarz: reordering list" parameter list.  Currently, that list
00714   // gets extracted in setup().
00715 
00716   // Subdomain check
00717   if (plist->isParameter ("schwarz: subdomain id") &&
00718       plist->get ("schwarz: subdomain id", -1) > 0) {
00719     UseSubdomain_ = true;
00720   }
00721   TEUCHOS_TEST_FOR_EXCEPTION(
00722     UseSubdomain_, std::logic_error, "Ifpack2::AdditiveSchwarz::"
00723     "setParameters: You specified the \"schwarz: subdomain id\" parameter, "
00724     "with a value other than -1.  This parameter is not yet supported.");
00725 
00726   // if true, filter singletons. NOTE: the filtered matrix can still have
00727   // singletons! A simple example: upper triangular matrix, if I remove
00728   // the lower node, I still get a matrix with a singleton! However, filter
00729   // singletons should help for PDE problems with Dirichlet BCs.
00730   FilterSingletons_ = plist->get ("schwarz: filter singletons", FilterSingletons_);
00731 
00732   // If the inner solver doesn't exist yet, don't create it.
00733   // initialize() creates it.
00734   //
00735   // If the inner solver _does_ exist, there are three cases,
00736   // depending on what the user put in the input ParameterList.
00737   //
00738   //   1. The user did /not/ provide a parameter specifying the inner
00739   //      solver's type, nor did the user specify a sublist of
00740   //      parameters for the inner solver
00741   //   2. The user did /not/ provide a parameter specifying the inner
00742   //      solver's type, but /did/ specify a sublist of parameters for
00743   //      the inner solver
00744   //   3. The user provided a parameter specifying the inner solver's
00745   //      type (it does not matter in this case whether the user gave
00746   //      a sublist of parameters for the inner solver)
00747   //
00748   // AdditiveSchwarz has "delta" (relative) semantics for setting
00749   // parameters.  This means that if the user did not specify the
00750   // inner solver's type, we presume that the type has not changed.
00751   // Thus, if the inner solver exists, we don't need to recreate it.
00752   //
00753   // In Case 3, if the user bothered to specify the inner solver's
00754   // type, then we must assume it may differ than the current inner
00755   // solver's type.  Thus, we have to recreate the inner solver.  We
00756   // achieve this here by assigning null to Inverse_; initialize()
00757   // will recreate the solver when it is needed.  Our assumption here
00758   // is necessary because Ifpack2::Preconditioner does not have a
00759   // method for querying a preconditioner's "type" (i.e., name) as a
00760   // string.  Remember that the user may have previously set an
00761   // arbitrary inner solver by calling setInnerPreconditioner().
00762   //
00763   // See note at the end of setInnerPreconditioner().
00764 
00765   if (! Inverse_.is_null ()) {
00766     // "CUSTOM" explicitly indicates that the user called or plans to
00767     // call setInnerPreconditioner.
00768     if (hasInnerPrecName () && innerPrecName () != "CUSTOM") {
00769       // Wipe out the current inner solver.  initialize() will
00770       // recreate it with the correct type.
00771       Inverse_ = Teuchos::null;
00772     }
00773     else {
00774       // Extract and apply the sublist of parameters to give to the
00775       // inner solver, if there is such a sublist of parameters.
00776       std::pair<Teuchos::ParameterList, bool> result = innerPrecParams ();
00777       if (result.second) {
00778         Inverse_->setParameters (result.first);
00779       }
00780     }
00781   }
00782 }
00783 
00784 
00785 
00786 template<class MatrixType,class LocalInverseType>
00787 Teuchos::RCP<const Teuchos::ParameterList>
00788 AdditiveSchwarz<MatrixType,LocalInverseType>::
00789 getValidParameters () const
00790 {
00791   using Teuchos::ParameterList;
00792   using Teuchos::parameterList;
00793   using Teuchos::RCP;
00794   using Teuchos::rcp_const_cast;
00795 
00796   if (validParams_.is_null ()) {
00797     const int overlapLevel = 0;
00798     const bool useReordering = false;
00799     const bool computeCondest = false;
00800     const bool filterSingletons = false;
00801     ParameterList reorderingSublist;
00802     reorderingSublist.set ("order_method", std::string ("rcm"));
00803 
00804     RCP<ParameterList> plist = parameterList ("Ifpack2::AdditiveSchwarz");
00805 
00806     Tpetra::setCombineModeParameter (*plist, "schwarz: combine mode");
00807     plist->set ("schwarz: overlap level", overlapLevel);
00808     plist->set ("schwarz: use reordering", useReordering);
00809     plist->set ("schwarz: reordering list", reorderingSublist);
00810     plist->set ("schwarz: compute condest", computeCondest);
00811     plist->set ("schwarz: filter singletons", filterSingletons);
00812 
00813     // FIXME (mfh 18 Nov 2013) Get valid parameters from inner solver.
00814     //
00815     // FIXME (mfh 18 Nov 2013) Get valid parameters from Zoltan2, if
00816     // Zoltan2 was enabled in the build.
00817 
00818     validParams_ = rcp_const_cast<const ParameterList> (plist);
00819   }
00820   return validParams_;
00821 }
00822 
00823 
00824 template<class MatrixType,class LocalInverseType>
00825 void AdditiveSchwarz<MatrixType,LocalInverseType>::initialize ()
00826 {
00827   using Tpetra::global_size_t;
00828   using Teuchos::RCP;
00829   using Teuchos::rcp;
00830   using Teuchos::SerialComm;
00831   using Teuchos::Time;
00832   using Teuchos::TimeMonitor;
00833 
00834   const std::string timerName ("Ifpack2::AdditiveSchwarz::initialize");
00835   RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
00836   if (timer.is_null ()) {
00837     timer = TimeMonitor::getNewCounter (timerName);
00838   }
00839 
00840   { // Start timing here.
00841     TimeMonitor timeMon (*timer);
00842 
00843     TEUCHOS_TEST_FOR_EXCEPTION(
00844       Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
00845       "initialize: The matrix to precondition is null.  You must either pass "
00846       "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
00847       "input, before you may call this method.");
00848 
00849     IsInitialized_ = false;
00850     IsComputed_ = false;
00851     Condest_ = -Teuchos::ScalarTraits<magnitude_type>::one ();
00852 
00853     RCP<const Teuchos::Comm<int> > comm = Matrix_->getComm ();
00854     RCP<const map_type> rowMap = Matrix_->getRowMap ();
00855     RCP<node_type> node = Matrix_->getNode ();
00856     const global_size_t INVALID =
00857       Teuchos::OrdinalTraits<global_size_t>::invalid ();
00858 
00859     // If there's only one process in the matrix's communicator,
00860     // then there's no need to compute overlap.
00861     if (comm->getSize () == 1) {
00862       OverlapLevel_ = 0;
00863       IsOverlapping_ = false;
00864     } else if (OverlapLevel_ != 0) {
00865       IsOverlapping_ = true;
00866     }
00867 
00868     if (OverlapLevel_ == 0) {
00869       const global_ordinal_type indexBase = rowMap->getIndexBase ();
00870       RCP<const SerialComm<int> > localComm (new SerialComm<int> ());
00871       // FIXME (mfh 15 Apr 2014) What if indexBase isn't the least
00872       // global index in the list of GIDs on this process?
00873       localMap_ =
00874         rcp (new map_type (INVALID, rowMap->getNodeNumElements (),
00875                            indexBase, localComm, node));
00876     }
00877 
00878     // compute the overlapping matrix if necessary
00879     if (IsOverlapping_) {
00880       if (UseSubdomain_) {
00881         const int sid = List_.get ("subdomain id", -1);
00882         OverlappingMatrix_ = rcp (new OverlappingRowMatrix<row_matrix_type> (Matrix_, OverlapLevel_, sid));
00883       } else {
00884         OverlappingMatrix_ = rcp (new OverlappingRowMatrix<row_matrix_type> (Matrix_, OverlapLevel_));
00885       }
00886     }
00887 
00888     setup (); // This does a lot of the initialization work.
00889 
00890     if (! Inverse_.is_null ()) {
00891       Inverse_->initialize (); // Initialize subdomain solver.
00892     }
00893   } // Stop timing here.
00894 
00895   IsInitialized_ = true;
00896   ++NumInitialize_;
00897 
00898   // timer->totalElapsedTime() returns the total time over all timer
00899   // calls.  Thus, we use = instead of +=.
00900   InitializeTime_ = timer->totalElapsedTime ();
00901 }
00902 
00903 
00904 template<class MatrixType,class LocalInverseType>
00905 bool AdditiveSchwarz<MatrixType,LocalInverseType>::isInitialized () const
00906 {
00907   return IsInitialized_;
00908 }
00909 
00910 
00911 template<class MatrixType,class LocalInverseType>
00912 void AdditiveSchwarz<MatrixType,LocalInverseType>::compute ()
00913 {
00914   using Teuchos::RCP;
00915   using Teuchos::Time;
00916   using Teuchos::TimeMonitor;
00917 
00918   if (! IsInitialized_) {
00919     initialize ();
00920   }
00921 
00922   TEUCHOS_TEST_FOR_EXCEPTION(
00923     ! isInitialized (), std::logic_error, "Ifpack2::AdditiveSchwarz::compute: "
00924     "The preconditioner is not yet initialized, "
00925     "even though initialize() supposedly has been called.  "
00926     "This should never happen.  "
00927     "Please report this bug to the Ifpack2 developers.");
00928 
00929   TEUCHOS_TEST_FOR_EXCEPTION(
00930     Inverse_.is_null (), std::runtime_error,
00931     "Ifpack2::AdditiveSchwarz::compute: The subdomain solver is null.  "
00932     "This can only happen if you called setInnerPreconditioner() with a null "
00933     "input, after calling initialize() or compute().  If you choose to call "
00934     "setInnerPreconditioner() with a null input, you must then call it with a "
00935     "nonnull input before you may call initialize() or compute().");
00936 
00937   const std::string timerName ("Ifpack2::AdditiveSchwarz::compute");
00938   RCP<Time> timer = TimeMonitor::lookupCounter (timerName);
00939   if (timer.is_null ()) {
00940     timer = TimeMonitor::getNewCounter (timerName);
00941   }
00942 
00943   { // Start timing here.
00944     TimeMonitor timeMon (*timer);
00945 
00946     IsComputed_ = false;
00947     Condest_ = -Teuchos::ScalarTraits<magnitude_type>::one ();
00948     Inverse_->compute ();
00949   } // Stop timing here.
00950 
00951   IsComputed_ = true;
00952   ++NumCompute_;
00953 
00954   // timer->totalElapsedTime() returns the total time over all timer
00955   // calls.  Thus, we use = instead of +=.
00956   ComputeTime_ = timer->totalElapsedTime ();
00957 }
00958 
00959 //==============================================================================
00960 // Returns true if the  preconditioner has been successfully computed, false otherwise.
00961 template<class MatrixType,class LocalInverseType>
00962 bool AdditiveSchwarz<MatrixType,LocalInverseType>::isComputed() const
00963 {
00964   return IsComputed_;
00965 }
00966 
00967 
00968 template<class MatrixType,class LocalInverseType>
00969 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00970 AdditiveSchwarz<MatrixType,LocalInverseType>::
00971 computeCondEst (CondestType CT,
00972                 local_ordinal_type MaxIters,
00973                 magnitude_type Tol,
00974                 const Teuchos::Ptr<const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > &Matrix_in)
00975 {
00976   // The preconditioner must have been computed in order to estimate
00977   // its condition number.
00978   if (! isComputed ()) {
00979     return -Teuchos::ScalarTraits<magnitude_type>::one ();
00980   }
00981 
00982   Condest_ = Ifpack2::Condest (*this, CT, MaxIters, Tol, Matrix_in);
00983   return Condest_;
00984 }
00985 
00986 
00987 template<class MatrixType,class LocalInverseType>
00988 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00989 AdditiveSchwarz<MatrixType,LocalInverseType>::getCondEst() const
00990 {
00991   return Condest_;
00992 }
00993 
00994 
00995 template<class MatrixType,class LocalInverseType>
00996 int AdditiveSchwarz<MatrixType,LocalInverseType>::getNumInitialize() const
00997 {
00998   return NumInitialize_;
00999 }
01000 
01001 
01002 template<class MatrixType,class LocalInverseType>
01003 int AdditiveSchwarz<MatrixType,LocalInverseType>::getNumCompute() const
01004 {
01005   return NumCompute_;
01006 }
01007 
01008 
01009 template<class MatrixType,class LocalInverseType>
01010 int AdditiveSchwarz<MatrixType,LocalInverseType>::getNumApply() const
01011 {
01012   return NumApply_;
01013 }
01014 
01015 
01016 template<class MatrixType,class LocalInverseType>
01017 double AdditiveSchwarz<MatrixType,LocalInverseType>::getInitializeTime() const
01018 {
01019   return InitializeTime_;
01020 }
01021 
01022 
01023 template<class MatrixType,class LocalInverseType>
01024 double AdditiveSchwarz<MatrixType,LocalInverseType>::getComputeTime() const
01025 {
01026   return ComputeTime_;
01027 }
01028 
01029 
01030 template<class MatrixType,class LocalInverseType>
01031 double AdditiveSchwarz<MatrixType,LocalInverseType>::getApplyTime() const
01032 {
01033   return ApplyTime_;
01034 }
01035 
01036 
01037 template<class MatrixType,class LocalInverseType>
01038 std::string AdditiveSchwarz<MatrixType,LocalInverseType>::description () const
01039 {
01040   std::ostringstream out;
01041 
01042   out << "\"Ifpack2::AdditiveSchwarz\": {";
01043   if (this->getObjectLabel () != "") {
01044     out << "Label: \"" << this->getObjectLabel () << "\"";
01045   }
01046   out << ", Initialized: " << (isInitialized () ? "true" : "false")
01047       << ", Computed: " << (isComputed () ? "true" : "false")
01048       << ", Overlap level: " << OverlapLevel_
01049       << ", Subdomain reordering: \"" << ReorderingAlgorithm_ << "\"";
01050 
01051   if (Matrix_.is_null ()) {
01052     out << ", Matrix: null";
01053   }
01054   else {
01055     out << ", Matrix: not null"
01056         << ", Global matrix dimensions: ["
01057         << Matrix_->getGlobalNumRows () << ", "
01058         << Matrix_->getGlobalNumCols () << "]";
01059   }
01060 
01061   // It's useful to print this instance's overlap level.  If you want
01062   // more info than that, call describe() instead.
01063   out << ", Overlap level: " << OverlapLevel_;
01064 
01065   out << "}";
01066   return out.str ();
01067 }
01068 
01069 
01070 template<class MatrixType,class LocalInverseType>
01071 void
01072 AdditiveSchwarz<MatrixType,LocalInverseType>::
01073 describe (Teuchos::FancyOStream& out,
01074           const Teuchos::EVerbosityLevel verbLevel) const
01075 {
01076   using Teuchos::OSTab;
01077   using Teuchos::TypeNameTraits;
01078   using std::endl;
01079 
01080   const int myRank = Matrix_->getComm ()->getRank ();
01081   const int numProcs = Matrix_->getComm ()->getSize ();
01082   const Teuchos::EVerbosityLevel vl =
01083     (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
01084 
01085   if (vl > Teuchos::VERB_NONE) {
01086     // describe() starts with a tab, by convention.
01087     OSTab tab0 (out);
01088     if (myRank == 0) {
01089       out << "\"Ifpack2::AdditiveSchwarz\":";
01090     }
01091     OSTab tab1 (out);
01092     if (myRank == 0) {
01093       out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
01094       out << "LocalInverseType: " << TypeNameTraits<LocalInverseType>::name () << endl;
01095       if (this->getObjectLabel () != "") {
01096         out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
01097       }
01098 
01099       out << "Overlap level: " << OverlapLevel_ << endl
01100           << "Combine mode: \"";
01101       if (CombineMode_ == Tpetra::INSERT) {
01102         out << "INSERT";
01103       } else if (CombineMode_ == Tpetra::ADD) {
01104         out << "ADD";
01105       } else if (CombineMode_ == Tpetra::REPLACE) {
01106         out << "REPLACE";
01107       } else if (CombineMode_ == Tpetra::ABSMAX) {
01108         out << "ABSMAX";
01109       } else if (CombineMode_ == Tpetra::ZERO) {
01110         out << "ZERO";
01111       }
01112       out << "\"" << endl
01113           << "Subdomain reordering: \"" << ReorderingAlgorithm_ << "\"" << endl;
01114     }
01115 
01116     if (Matrix_.is_null ()) {
01117       if (myRank == 0) {
01118         out << "Matrix: null" << endl;
01119       }
01120     }
01121     else {
01122       if (myRank == 0) {
01123         out << "Matrix:" << endl;
01124         std::flush (out);
01125       }
01126       Matrix_->getComm ()->barrier (); // wait for output to finish
01127       Matrix_->describe (out, Teuchos::VERB_LOW);
01128     }
01129 
01130     if (myRank == 0) {
01131       out << "Number of initialize calls: " << getNumInitialize () << endl
01132           << "Number of compute calls: " << getNumCompute () << endl
01133           << "Number of apply calls: " << getNumApply () << endl
01134           << "Total time in seconds for initialize: " << getInitializeTime () << endl
01135           << "Total time in seconds for compute: " << getComputeTime () << endl
01136           << "Total time in seconds for apply: " << getApplyTime () << endl;
01137     }
01138 
01139     if (Inverse_.is_null ()) {
01140       if (myRank == 0) {
01141         out << "Subdomain solver: null" << endl;
01142       }
01143     }
01144     else {
01145       if (vl < Teuchos::VERB_EXTREME) {
01146         if (myRank == 0) {
01147           out << "Subdomain solver: not null" << endl;
01148         }
01149       }
01150       else { // vl >= Teuchos::VERB_EXTREME
01151         for (int p = 0; p < numProcs; ++p) {
01152           if (p == myRank) {
01153             out << "Subdomain solver on Process " << myRank << ":";
01154             if (Inverse_.is_null ()) {
01155               out << "null" << endl;
01156             } else {
01157               out << endl;
01158               Inverse_->describe (out, vl);
01159             }
01160           }
01161           Matrix_->getComm ()->barrier ();
01162           Matrix_->getComm ()->barrier ();
01163           Matrix_->getComm ()->barrier (); // wait for output to finish
01164         }
01165       }
01166     }
01167 
01168     Matrix_->getComm ()->barrier (); // wait for output to finish
01169   }
01170 }
01171 
01172 
01173 template<class MatrixType,class LocalInverseType>
01174 std::ostream& AdditiveSchwarz<MatrixType,LocalInverseType>::print(std::ostream& os) const
01175 {
01176   Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
01177   fos.setOutputToRootOnly(0);
01178   describe(fos);
01179   return(os);
01180 }
01181 
01182 
01183 template<class MatrixType,class LocalInverseType>
01184 int AdditiveSchwarz<MatrixType,LocalInverseType>::getOverlapLevel() const
01185 {
01186   return OverlapLevel_;
01187 }
01188 
01189 
01190 template<class MatrixType,class LocalInverseType>
01191 void AdditiveSchwarz<MatrixType,LocalInverseType>::setup ()
01192 {
01193 #ifdef HAVE_MPI
01194   using Teuchos::MpiComm;
01195 #endif // HAVE_MPI
01196   using Teuchos::ArrayRCP;
01197   using Teuchos::RCP;
01198   using Teuchos::rcp;
01199   using Teuchos::rcp_dynamic_cast;
01200   using Teuchos::rcpFromRef;
01201 
01202 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
01203   typedef Xpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> XpetraMatrixType;
01204   typedef Xpetra::TpetraRowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> XpetraTpetraMatrixType;
01205 #endif
01206 
01207   TEUCHOS_TEST_FOR_EXCEPTION(
01208     Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
01209     "initialize: The matrix to precondition is null.  You must either pass "
01210     "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
01211     "input, before you may call this method.");
01212 
01213   // Localized version of Matrix_ or OverlappingMatrix_.
01214   RCP<row_matrix_type> LocalizedMatrix;
01215 
01216   // The "most current local matrix."  At the end of this method, this
01217   // will be handed off to the inner solver.
01218   RCP<row_matrix_type> ActiveMatrix;
01219 
01220   // Create localized matrix.
01221   if (! OverlappingMatrix_.is_null ()) {
01222     if (UseSubdomain_) {
01223       //      int sid = List_.get("subdomain id",-1);
01224       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
01225         "Ifpack2::AdditiveSchwarz::setup: subdomain code not yet supported.");
01226       //
01227       // FIXME (mfh 18 Nov 2013) btw what's the difference between
01228       // Ifpack_NodeFilter and Ifpack_LocalFilter?  The former's
01229       // documentation sounds a lot like what Ifpack2::LocalFilter
01230       // does.
01231       //
01232       //Ifpack2_NodeFilter *tt = new Ifpack2_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
01233       //LocalizedMatrix = Teuchos::rcp(tt);
01234     }
01235     else
01236       LocalizedMatrix = rcp (new LocalFilter<row_matrix_type> (OverlappingMatrix_));
01237   }
01238   else {
01239     if (UseSubdomain_) {
01240       //      int sid = List_.get("subdomain id",-1);
01241       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
01242         "Ifpack2::AdditiveSchwarz::setup: subdomain code not yet supported.");
01243     }
01244     else {
01245       LocalizedMatrix = rcp (new LocalFilter<row_matrix_type> (Matrix_));
01246     }
01247   }
01248 
01249   // Sanity check; I don't trust the logic above to have created LocalizedMatrix.
01250   TEUCHOS_TEST_FOR_EXCEPTION(
01251     LocalizedMatrix.is_null (), std::logic_error,
01252     "Ifpack2::AdditiveSchwarz::setup: LocalizedMatrix is null, after the code "
01253     "that claimed to have created it.  This should never be the case.  Please "
01254     "report this bug to the Ifpack2 developers.");
01255 
01256   // Mark localized matrix as active
01257   ActiveMatrix = LocalizedMatrix;
01258 
01259   // Singleton Filtering
01260   if (FilterSingletons_) {
01261     SingletonMatrix_ = rcp (new SingletonFilter<row_matrix_type> (LocalizedMatrix));
01262     ActiveMatrix = SingletonMatrix_;
01263   }
01264 
01265   // Do reordering
01266   if (UseReordering_) {
01267 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
01268     // Unlike Ifpack, Zoltan2 does all the dirty work here.
01269     Teuchos::ParameterList zlist = List_.sublist ("schwarz: reordering list");
01270 
01271     // FIXME (mfh 18 Nov 2013) Shouldn't this come from the Zoltan2 sublist?
01272     ReorderingAlgorithm_ = List_.get<std::string> ("order_method", "rcm");
01273     XpetraTpetraMatrixType XpetraMatrix (ActiveMatrix);
01274     typedef Zoltan2::XpetraRowMatrixAdapter<XpetraMatrixType> z2_adapter_type;
01275     z2_adapter_type Zoltan2Matrix (rcpFromRef (XpetraMatrix));
01276     typedef Zoltan2::OrderingProblem<z2_adapter_type> ordering_problem_type;
01277 #ifdef HAVE_MPI
01278     // Grab the MPI Communicator and build the ordering problem with that
01279     MPI_Comm myRawComm;
01280 
01281     RCP<const MpiComm<int> > mpicomm =
01282       rcp_dynamic_cast<const MpiComm<int> > (ActiveMatrix->getComm ());
01283     if (mpicomm == Teuchos::null) {
01284       myRawComm = MPI_COMM_SELF;
01285     } else {
01286       myRawComm = * (mpicomm->getRawMpiComm ());
01287     }
01288     ordering_problem_type MyOrderingProblem (&Zoltan2Matrix, &zlist, myRawComm);
01289 #else
01290     ordering_problem_type MyOrderingProblem (&Zoltan2Matrix, &zlist);
01291 #endif
01292     MyOrderingProblem.solve ();
01293 
01294     // Now create the reordered matrix & mark it as active
01295     {
01296       typedef ReorderFilter<row_matrix_type> reorder_filter_type;
01297       typedef Zoltan2::OrderingSolution<global_ordinal_type,
01298         local_ordinal_type> ordering_solution_type;
01299 
01300       ordering_solution_type sol (*MyOrderingProblem.getSolution ());
01301 
01302       // perm[i] gives the where OLD index i shows up in the NEW
01303       // ordering.  revperm[i] gives the where NEW index i shows
01304       // up in the OLD ordering.  Note that perm is actually the
01305       // "inverse permutation," in Zoltan2 terms.
01306       ArrayRCP<local_ordinal_type> perm = sol.getPermutationRCPConst (true);
01307       ArrayRCP<local_ordinal_type> revperm = sol.getPermutationRCPConst ();
01308 
01309       ReorderedLocalizedMatrix_ =
01310         rcp (new reorder_filter_type (ActiveMatrix, perm, revperm));
01311 
01312       ActiveMatrix = ReorderedLocalizedMatrix_;
01313     }
01314 #else
01315     // This is a logic_error, not a runtime_error, because
01316     // setParameters() should have excluded this case already.
01317     TEUCHOS_TEST_FOR_EXCEPTION(
01318       true, std::logic_error, "Ifpack2::AdditiveSchwarz::setup: "
01319       "The Zoltan2 and Xpetra packages must be enabled in order "
01320       "to support reordering.");
01321 #endif
01322   }
01323 
01324   innerMatrix_ = ActiveMatrix;
01325 
01326   TEUCHOS_TEST_FOR_EXCEPTION(
01327     innerMatrix_.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
01328     "setup: Inner matrix is null right before constructing inner solver.  "
01329     "Please report this bug to the Ifpack2 developers.");
01330 
01331   // Construct the inner solver if necessary.
01332   if (Inverse_.is_null ()) {
01333     const std::string innerName = innerPrecName ();
01334     TEUCHOS_TEST_FOR_EXCEPTION(
01335       innerName == "INVALID", std::logic_error,
01336       "Ifpack2::AdditiveSchwarz::initialize: AdditiveSchwarz doesn't "
01337       "know how to create an instance of your LocalInverseType \""
01338       << Teuchos::TypeNameTraits<LocalInverseType>::name () << "\".  If "
01339       "LocalInverseType is a single-level preconditioner (does not take an "
01340       "inner solver), then you can fix this in one of two ways.  Either (a) "
01341       "create the LocalInverseType instance yourself and give it to "
01342       "AdditiveSchwarz by calling setInnerPreconditioner(), before calling "
01343       "initialize(), or (b) teach Details::OneLevelFactory how to create an "
01344       "inner preconditioner of that type.  "
01345       "Please talk to the Ifpack2 developers for details.");
01346 
01347     TEUCHOS_TEST_FOR_EXCEPTION(
01348       innerName == "CUSTOM", std::runtime_error, "Ifpack2::AdditiveSchwarz::"
01349       "initialize: If the \"inner preconditioner name\" parameter (or any "
01350       "alias thereof) has the value \"CUSTOM\", then you must first call "
01351       "setInnerPreconditioner with a nonnull inner preconditioner input before "
01352       "you may call initialize().");
01353 
01354     Details::OneLevelFactory<MatrixType> factory;
01355     RCP<prec_type> innerPrec = factory.create (innerName, innerMatrix_);
01356     TEUCHOS_TEST_FOR_EXCEPTION(
01357       innerPrec.is_null (), std::logic_error,
01358       "Ifpack2::AdditiveSchwarz::setup: Failed to create inner preconditioner "
01359       "with name \"" << innerName << "\".");
01360 
01361     // Extract and apply the sublist of parameters to give to the
01362     // inner solver, if there is such a sublist of parameters.
01363     std::pair<Teuchos::ParameterList, bool> result = innerPrecParams ();
01364     if (result.second) {
01365       innerPrec->setParameters (result.first);
01366     }
01367     Inverse_ = innerPrec; // "Commit" the inner solver.
01368   }
01369   else if (Inverse_->getMatrix ().getRawPtr () != innerMatrix_.getRawPtr ()) {
01370     // The new inner matrix is different from the inner
01371     // preconditioner's current matrix, so give the inner
01372     // preconditioner the new inner matrix.  First make sure that the
01373     // inner solver knows how to have its matrix changed.
01374     typedef Details::CanChangeMatrix<row_matrix_type> can_change_type;
01375     can_change_type* innerSolver =
01376       dynamic_cast<can_change_type*> (Inverse_.getRawPtr ());
01377     TEUCHOS_TEST_FOR_EXCEPTION(
01378       innerSolver == NULL, std::invalid_argument, "Ifpack2::AdditiveSchwarz::"
01379       "setup: The current inner preconditioner does not implement the "
01380       "setMatrix() feature.  Only preconditioners that inherit from "
01381       "Ifpack2::Details::CanChangeMatrix implement this feature.");
01382 
01383     // Give the new inner matrix to the inner preconditioner.
01384     innerSolver->setMatrix (innerMatrix_);
01385   }
01386   TEUCHOS_TEST_FOR_EXCEPTION(
01387     Inverse_.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
01388     "setup: Inverse_ is null right after we were supposed to have created it."
01389     "  Please report this bug to the Ifpack2 developers.");
01390 
01391   // We don't have to call setInnerPreconditioner() here, because we
01392   // had the inner matrix (innerMatrix_) before creation of the inner
01393   // preconditioner.  Calling setInnerPreconditioner here would be
01394   // legal, but it would require an unnecessary reset of the inner
01395   // preconditioner (i.e., calling initialize() and compute() again).
01396 }
01397 
01398 
01399 template<class MatrixType, class LocalInverseType>
01400 void AdditiveSchwarz<MatrixType, LocalInverseType>::
01401 setInnerPreconditioner (const Teuchos::RCP<Preconditioner<scalar_type,
01402                                                           local_ordinal_type,
01403                                                           global_ordinal_type,
01404                                                           node_type> >& innerPrec)
01405 {
01406   if (! innerPrec.is_null ()) {
01407     // Make sure that the new inner solver knows how to have its matrix changed.
01408     typedef Details::CanChangeMatrix<row_matrix_type> can_change_type;
01409     can_change_type* innerSolver = dynamic_cast<can_change_type*> (&*innerPrec);
01410     TEUCHOS_TEST_FOR_EXCEPTION(
01411       innerSolver == NULL, std::invalid_argument, "Ifpack2::AdditiveSchwarz::"
01412       "setInnerPreconditioner: The input preconditioner does not implement the "
01413       "setMatrix() feature.  Only input preconditioners that inherit from "
01414       "Ifpack2::Details::CanChangeMatrix implement this feature.");
01415 
01416     // If users provide an inner solver, we assume that
01417     // AdditiveSchwarz's current inner solver parameters no longer
01418     // apply.  (In fact, we will remove those parameters from
01419     // AdditiveSchwarz's current list below.)  Thus, we do /not/ apply
01420     // the current sublist of inner solver parameters to the input
01421     // inner solver.
01422 
01423     // mfh 03 Jan 2014: Thanks to Paul Tsuji for pointing out that it's
01424     // perfectly legal for innerMatrix_ to be null here.  This can
01425     // happen if initialize() has not been called yet.  For example,
01426     // when Factory creates an AdditiveSchwarz instance, it calls
01427     // setInnerPreconditioner() without first calling initialize().
01428 
01429     // Give the local matrix to the new inner solver.
01430     innerSolver->setMatrix (innerMatrix_);
01431 
01432     // If the user previously specified a parameter for the inner
01433     // preconditioner's type, then clear out that parameter and its
01434     // associated sublist.  Replace the inner preconditioner's type with
01435     // "CUSTOM", to make it obvious that AdditiveSchwarz's ParameterList
01436     // does not necessarily describe the current inner preconditioner.
01437     // We have to remove all allowed aliases of "inner preconditioner
01438     // name" before we may set it to "CUSTOM".  Users may also set this
01439     // parameter to "CUSTOM" themselves, but this is not required.
01440     removeInnerPrecName ();
01441     removeInnerPrecParams ();
01442     List_.set ("inner preconditioner name", "CUSTOM");
01443 
01444     // Bring the new inner solver's current status (initialized or
01445     // computed) in line with AdditiveSchwarz's current status.
01446     if (isInitialized ()) {
01447       innerPrec->initialize ();
01448     }
01449     if (isComputed ()) {
01450       innerPrec->compute ();
01451     }
01452   }
01453 
01454   // If the new inner solver is null, we don't change the initialized
01455   // or computed status of AdditiveSchwarz.  That way, AdditiveSchwarz
01456   // won't have to recompute innerMatrix_ if the inner solver changes.
01457   // This does introduce a new error condition in compute() and
01458   // apply(), but that's OK.
01459 
01460   // Set the new inner solver.
01461   Inverse_ = innerPrec;
01462 }
01463 
01464 template<class MatrixType, class LocalInverseType>
01465 void AdditiveSchwarz<MatrixType, LocalInverseType>::
01466 setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
01467 {
01468   // Don't set the matrix unless it is different from the current one.
01469   if (A.getRawPtr () != Matrix_.getRawPtr ()) {
01470     IsInitialized_ = false;
01471     IsComputed_ = false;
01472 
01473     // Reset all the state computed in initialize() and compute().
01474     Condest_ = -Teuchos::ScalarTraits<magnitude_type>::one ();
01475     OverlappingMatrix_ = Teuchos::null;
01476     ReorderedLocalizedMatrix_ = Teuchos::null;
01477     innerMatrix_ = Teuchos::null;
01478     SingletonMatrix_ = Teuchos::null;
01479     localMap_ = Teuchos::null;
01480     DistributedImporter_ = Teuchos::null;
01481 
01482     Matrix_ = A;
01483   }
01484 }
01485 
01486 } // namespace Ifpack2
01487 
01488 #define IFPACK2_ADDITIVESCHWARZ_INSTANT(S,LO,GO,N) \
01489   template class Ifpack2::AdditiveSchwarz< Tpetra::CrsMatrix<S, LO, GO, N> >; \
01490   template class Ifpack2::AdditiveSchwarz< Tpetra::CrsMatrix<S, LO, GO, N>, \
01491                                   Ifpack2::ILUT<Tpetra::CrsMatrix< S, LO, GO, N > > >;
01492 
01493 #endif // IFPACK2_ADDITIVESCHWARZ_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends