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   out << ", Combine mode: \"";
01051   if (CombineMode_ == Tpetra::INSERT) {
01052     out << "INSERT";
01053   } else if (CombineMode_ == Tpetra::ADD) {
01054     out << "ADD";
01055   } else if (CombineMode_ == Tpetra::REPLACE) {
01056     out << "REPLACE";
01057   } else if (CombineMode_ == Tpetra::ABSMAX) {
01058     out << "ABSMAX";
01059   } else if (CombineMode_ == Tpetra::ZERO) {
01060     out << "ZERO";
01061   }
01062   out << "\"";
01063   if (Matrix_.is_null ()) {
01064     out << ", Matrix: null";
01065   }
01066   else {
01067     out << ", Global matrix dimensions: ["
01068         << Matrix_->getGlobalNumRows () << ", "
01069         << Matrix_->getGlobalNumCols () << "]";
01070   }
01071   out << ", Inner solver: ";
01072   if (!Inverse_.is_null ())
01073     out << "{" << Inverse_->description() << "}";
01074   else
01075     out << "null";
01076 
01077   out << "}";
01078   return out.str ();
01079 }
01080 
01081 
01082 template<class MatrixType,class LocalInverseType>
01083 void
01084 AdditiveSchwarz<MatrixType,LocalInverseType>::
01085 describe (Teuchos::FancyOStream& out,
01086           const Teuchos::EVerbosityLevel verbLevel) const
01087 {
01088   using Teuchos::OSTab;
01089   using Teuchos::TypeNameTraits;
01090   using std::endl;
01091 
01092   const int myRank = Matrix_->getComm ()->getRank ();
01093   const int numProcs = Matrix_->getComm ()->getSize ();
01094   const Teuchos::EVerbosityLevel vl =
01095     (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
01096 
01097   if (vl > Teuchos::VERB_NONE) {
01098     // describe() starts with a tab, by convention.
01099     OSTab tab0 (out);
01100     if (myRank == 0) {
01101       out << "\"Ifpack2::AdditiveSchwarz\":";
01102     }
01103     OSTab tab1 (out);
01104     if (myRank == 0) {
01105       out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
01106       out << "LocalInverseType: " << TypeNameTraits<LocalInverseType>::name () << endl;
01107       if (this->getObjectLabel () != "") {
01108         out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
01109       }
01110 
01111       out << "Overlap level: " << OverlapLevel_ << endl
01112           << "Combine mode: \"";
01113       if (CombineMode_ == Tpetra::INSERT) {
01114         out << "INSERT";
01115       } else if (CombineMode_ == Tpetra::ADD) {
01116         out << "ADD";
01117       } else if (CombineMode_ == Tpetra::REPLACE) {
01118         out << "REPLACE";
01119       } else if (CombineMode_ == Tpetra::ABSMAX) {
01120         out << "ABSMAX";
01121       } else if (CombineMode_ == Tpetra::ZERO) {
01122         out << "ZERO";
01123       }
01124       out << "\"" << endl
01125           << "Subdomain reordering: \"" << ReorderingAlgorithm_ << "\"" << endl;
01126     }
01127 
01128     if (Matrix_.is_null ()) {
01129       if (myRank == 0) {
01130         out << "Matrix: null" << endl;
01131       }
01132     }
01133     else {
01134       if (myRank == 0) {
01135         out << "Matrix:" << endl;
01136         std::flush (out);
01137       }
01138       Matrix_->getComm ()->barrier (); // wait for output to finish
01139       Matrix_->describe (out, Teuchos::VERB_LOW);
01140     }
01141 
01142     if (myRank == 0) {
01143       out << "Number of initialize calls: " << getNumInitialize () << endl
01144           << "Number of compute calls: " << getNumCompute () << endl
01145           << "Number of apply calls: " << getNumApply () << endl
01146           << "Total time in seconds for initialize: " << getInitializeTime () << endl
01147           << "Total time in seconds for compute: " << getComputeTime () << endl
01148           << "Total time in seconds for apply: " << getApplyTime () << endl;
01149     }
01150 
01151     if (Inverse_.is_null ()) {
01152       if (myRank == 0) {
01153         out << "Subdomain solver: null" << endl;
01154       }
01155     }
01156     else {
01157       if (vl < Teuchos::VERB_EXTREME) {
01158         if (myRank == 0) {
01159           out << "Subdomain solver: not null" << endl;
01160         }
01161       }
01162       else { // vl >= Teuchos::VERB_EXTREME
01163         for (int p = 0; p < numProcs; ++p) {
01164           if (p == myRank) {
01165             out << "Subdomain solver on Process " << myRank << ":";
01166             if (Inverse_.is_null ()) {
01167               out << "null" << endl;
01168             } else {
01169               out << endl;
01170               Inverse_->describe (out, vl);
01171             }
01172           }
01173           Matrix_->getComm ()->barrier ();
01174           Matrix_->getComm ()->barrier ();
01175           Matrix_->getComm ()->barrier (); // wait for output to finish
01176         }
01177       }
01178     }
01179 
01180     Matrix_->getComm ()->barrier (); // wait for output to finish
01181   }
01182 }
01183 
01184 
01185 template<class MatrixType,class LocalInverseType>
01186 std::ostream& AdditiveSchwarz<MatrixType,LocalInverseType>::print(std::ostream& os) const
01187 {
01188   Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
01189   fos.setOutputToRootOnly(0);
01190   describe(fos);
01191   return(os);
01192 }
01193 
01194 
01195 template<class MatrixType,class LocalInverseType>
01196 int AdditiveSchwarz<MatrixType,LocalInverseType>::getOverlapLevel() const
01197 {
01198   return OverlapLevel_;
01199 }
01200 
01201 
01202 template<class MatrixType,class LocalInverseType>
01203 void AdditiveSchwarz<MatrixType,LocalInverseType>::setup ()
01204 {
01205 #ifdef HAVE_MPI
01206   using Teuchos::MpiComm;
01207 #endif // HAVE_MPI
01208   using Teuchos::ArrayRCP;
01209   using Teuchos::RCP;
01210   using Teuchos::rcp;
01211   using Teuchos::rcp_dynamic_cast;
01212   using Teuchos::rcpFromRef;
01213 
01214 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
01215   typedef Xpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> XpetraMatrixType;
01216   typedef Xpetra::TpetraRowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> XpetraTpetraMatrixType;
01217 #endif
01218 
01219   TEUCHOS_TEST_FOR_EXCEPTION(
01220     Matrix_.is_null (), std::runtime_error, "Ifpack2::AdditiveSchwarz::"
01221     "initialize: The matrix to precondition is null.  You must either pass "
01222     "a nonnull matrix to the constructor, or call setMatrix() with a nonnull "
01223     "input, before you may call this method.");
01224 
01225   // Localized version of Matrix_ or OverlappingMatrix_.
01226   RCP<row_matrix_type> LocalizedMatrix;
01227 
01228   // The "most current local matrix."  At the end of this method, this
01229   // will be handed off to the inner solver.
01230   RCP<row_matrix_type> ActiveMatrix;
01231 
01232   // Create localized matrix.
01233   if (! OverlappingMatrix_.is_null ()) {
01234     if (UseSubdomain_) {
01235       //      int sid = List_.get("subdomain id",-1);
01236       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
01237         "Ifpack2::AdditiveSchwarz::setup: subdomain code not yet supported.");
01238       //
01239       // FIXME (mfh 18 Nov 2013) btw what's the difference between
01240       // Ifpack_NodeFilter and Ifpack_LocalFilter?  The former's
01241       // documentation sounds a lot like what Ifpack2::LocalFilter
01242       // does.
01243       //
01244       //Ifpack2_NodeFilter *tt = new Ifpack2_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
01245       //LocalizedMatrix = Teuchos::rcp(tt);
01246     }
01247     else
01248       LocalizedMatrix = rcp (new LocalFilter<row_matrix_type> (OverlappingMatrix_));
01249   }
01250   else {
01251     if (UseSubdomain_) {
01252       //      int sid = List_.get("subdomain id",-1);
01253       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
01254         "Ifpack2::AdditiveSchwarz::setup: subdomain code not yet supported.");
01255     }
01256     else {
01257       LocalizedMatrix = rcp (new LocalFilter<row_matrix_type> (Matrix_));
01258     }
01259   }
01260 
01261   // Sanity check; I don't trust the logic above to have created LocalizedMatrix.
01262   TEUCHOS_TEST_FOR_EXCEPTION(
01263     LocalizedMatrix.is_null (), std::logic_error,
01264     "Ifpack2::AdditiveSchwarz::setup: LocalizedMatrix is null, after the code "
01265     "that claimed to have created it.  This should never be the case.  Please "
01266     "report this bug to the Ifpack2 developers.");
01267 
01268   // Mark localized matrix as active
01269   ActiveMatrix = LocalizedMatrix;
01270 
01271   // Singleton Filtering
01272   if (FilterSingletons_) {
01273     SingletonMatrix_ = rcp (new SingletonFilter<row_matrix_type> (LocalizedMatrix));
01274     ActiveMatrix = SingletonMatrix_;
01275   }
01276 
01277   // Do reordering
01278   if (UseReordering_) {
01279 #if defined(HAVE_IFPACK2_XPETRA) && defined(HAVE_IFPACK2_ZOLTAN2)
01280     // Unlike Ifpack, Zoltan2 does all the dirty work here.
01281     Teuchos::ParameterList zlist = List_.sublist ("schwarz: reordering list");
01282 
01283     // FIXME (mfh 18 Nov 2013) Shouldn't this come from the Zoltan2 sublist?
01284     ReorderingAlgorithm_ = List_.get<std::string> ("order_method", "rcm");
01285     XpetraTpetraMatrixType XpetraMatrix (ActiveMatrix);
01286     typedef Zoltan2::XpetraRowMatrixAdapter<XpetraMatrixType> z2_adapter_type;
01287     z2_adapter_type Zoltan2Matrix (rcpFromRef (XpetraMatrix));
01288     typedef Zoltan2::OrderingProblem<z2_adapter_type> ordering_problem_type;
01289 #ifdef HAVE_MPI
01290     // Grab the MPI Communicator and build the ordering problem with that
01291     MPI_Comm myRawComm;
01292 
01293     RCP<const MpiComm<int> > mpicomm =
01294       rcp_dynamic_cast<const MpiComm<int> > (ActiveMatrix->getComm ());
01295     if (mpicomm == Teuchos::null) {
01296       myRawComm = MPI_COMM_SELF;
01297     } else {
01298       myRawComm = * (mpicomm->getRawMpiComm ());
01299     }
01300     ordering_problem_type MyOrderingProblem (&Zoltan2Matrix, &zlist, myRawComm);
01301 #else
01302     ordering_problem_type MyOrderingProblem (&Zoltan2Matrix, &zlist);
01303 #endif
01304     MyOrderingProblem.solve ();
01305 
01306     // Now create the reordered matrix & mark it as active
01307     {
01308       typedef ReorderFilter<row_matrix_type> reorder_filter_type;
01309       typedef Zoltan2::OrderingSolution<global_ordinal_type,
01310         local_ordinal_type> ordering_solution_type;
01311 
01312       ordering_solution_type sol (*MyOrderingProblem.getSolution ());
01313 
01314       // perm[i] gives the where OLD index i shows up in the NEW
01315       // ordering.  revperm[i] gives the where NEW index i shows
01316       // up in the OLD ordering.  Note that perm is actually the
01317       // "inverse permutation," in Zoltan2 terms.
01318       ArrayRCP<local_ordinal_type> perm = sol.getPermutationRCPConst (true);
01319       ArrayRCP<local_ordinal_type> revperm = sol.getPermutationRCPConst ();
01320 
01321       ReorderedLocalizedMatrix_ =
01322         rcp (new reorder_filter_type (ActiveMatrix, perm, revperm));
01323 
01324       ActiveMatrix = ReorderedLocalizedMatrix_;
01325     }
01326 #else
01327     // This is a logic_error, not a runtime_error, because
01328     // setParameters() should have excluded this case already.
01329     TEUCHOS_TEST_FOR_EXCEPTION(
01330       true, std::logic_error, "Ifpack2::AdditiveSchwarz::setup: "
01331       "The Zoltan2 and Xpetra packages must be enabled in order "
01332       "to support reordering.");
01333 #endif
01334   }
01335 
01336   innerMatrix_ = ActiveMatrix;
01337 
01338   TEUCHOS_TEST_FOR_EXCEPTION(
01339     innerMatrix_.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
01340     "setup: Inner matrix is null right before constructing inner solver.  "
01341     "Please report this bug to the Ifpack2 developers.");
01342 
01343   // Construct the inner solver if necessary.
01344   if (Inverse_.is_null ()) {
01345     const std::string innerName = innerPrecName ();
01346     TEUCHOS_TEST_FOR_EXCEPTION(
01347       innerName == "INVALID", std::logic_error,
01348       "Ifpack2::AdditiveSchwarz::initialize: AdditiveSchwarz doesn't "
01349       "know how to create an instance of your LocalInverseType \""
01350       << Teuchos::TypeNameTraits<LocalInverseType>::name () << "\".  If "
01351       "LocalInverseType is a single-level preconditioner (does not take an "
01352       "inner solver), then you can fix this in one of two ways.  Either (a) "
01353       "create the LocalInverseType instance yourself and give it to "
01354       "AdditiveSchwarz by calling setInnerPreconditioner(), before calling "
01355       "initialize(), or (b) teach Details::OneLevelFactory how to create an "
01356       "inner preconditioner of that type.  "
01357       "Please talk to the Ifpack2 developers for details.");
01358 
01359     TEUCHOS_TEST_FOR_EXCEPTION(
01360       innerName == "CUSTOM", std::runtime_error, "Ifpack2::AdditiveSchwarz::"
01361       "initialize: If the \"inner preconditioner name\" parameter (or any "
01362       "alias thereof) has the value \"CUSTOM\", then you must first call "
01363       "setInnerPreconditioner with a nonnull inner preconditioner input before "
01364       "you may call initialize().");
01365 
01366     Details::OneLevelFactory<MatrixType> factory;
01367     RCP<prec_type> innerPrec = factory.create (innerName, innerMatrix_);
01368     TEUCHOS_TEST_FOR_EXCEPTION(
01369       innerPrec.is_null (), std::logic_error,
01370       "Ifpack2::AdditiveSchwarz::setup: Failed to create inner preconditioner "
01371       "with name \"" << innerName << "\".");
01372 
01373     // Extract and apply the sublist of parameters to give to the
01374     // inner solver, if there is such a sublist of parameters.
01375     std::pair<Teuchos::ParameterList, bool> result = innerPrecParams ();
01376     if (result.second) {
01377       innerPrec->setParameters (result.first);
01378     }
01379     Inverse_ = innerPrec; // "Commit" the inner solver.
01380   }
01381   else if (Inverse_->getMatrix ().getRawPtr () != innerMatrix_.getRawPtr ()) {
01382     // The new inner matrix is different from the inner
01383     // preconditioner's current matrix, so give the inner
01384     // preconditioner the new inner matrix.  First make sure that the
01385     // inner solver knows how to have its matrix changed.
01386     typedef Details::CanChangeMatrix<row_matrix_type> can_change_type;
01387     can_change_type* innerSolver =
01388       dynamic_cast<can_change_type*> (Inverse_.getRawPtr ());
01389     TEUCHOS_TEST_FOR_EXCEPTION(
01390       innerSolver == NULL, std::invalid_argument, "Ifpack2::AdditiveSchwarz::"
01391       "setup: The current inner preconditioner does not implement the "
01392       "setMatrix() feature.  Only preconditioners that inherit from "
01393       "Ifpack2::Details::CanChangeMatrix implement this feature.");
01394 
01395     // Give the new inner matrix to the inner preconditioner.
01396     innerSolver->setMatrix (innerMatrix_);
01397   }
01398   TEUCHOS_TEST_FOR_EXCEPTION(
01399     Inverse_.is_null (), std::logic_error, "Ifpack2::AdditiveSchwarz::"
01400     "setup: Inverse_ is null right after we were supposed to have created it."
01401     "  Please report this bug to the Ifpack2 developers.");
01402 
01403   // We don't have to call setInnerPreconditioner() here, because we
01404   // had the inner matrix (innerMatrix_) before creation of the inner
01405   // preconditioner.  Calling setInnerPreconditioner here would be
01406   // legal, but it would require an unnecessary reset of the inner
01407   // preconditioner (i.e., calling initialize() and compute() again).
01408 }
01409 
01410 
01411 template<class MatrixType, class LocalInverseType>
01412 void AdditiveSchwarz<MatrixType, LocalInverseType>::
01413 setInnerPreconditioner (const Teuchos::RCP<Preconditioner<scalar_type,
01414                                                           local_ordinal_type,
01415                                                           global_ordinal_type,
01416                                                           node_type> >& innerPrec)
01417 {
01418   if (! innerPrec.is_null ()) {
01419     // Make sure that the new inner solver knows how to have its matrix changed.
01420     typedef Details::CanChangeMatrix<row_matrix_type> can_change_type;
01421     can_change_type* innerSolver = dynamic_cast<can_change_type*> (&*innerPrec);
01422     TEUCHOS_TEST_FOR_EXCEPTION(
01423       innerSolver == NULL, std::invalid_argument, "Ifpack2::AdditiveSchwarz::"
01424       "setInnerPreconditioner: The input preconditioner does not implement the "
01425       "setMatrix() feature.  Only input preconditioners that inherit from "
01426       "Ifpack2::Details::CanChangeMatrix implement this feature.");
01427 
01428     // If users provide an inner solver, we assume that
01429     // AdditiveSchwarz's current inner solver parameters no longer
01430     // apply.  (In fact, we will remove those parameters from
01431     // AdditiveSchwarz's current list below.)  Thus, we do /not/ apply
01432     // the current sublist of inner solver parameters to the input
01433     // inner solver.
01434 
01435     // mfh 03 Jan 2014: Thanks to Paul Tsuji for pointing out that it's
01436     // perfectly legal for innerMatrix_ to be null here.  This can
01437     // happen if initialize() has not been called yet.  For example,
01438     // when Factory creates an AdditiveSchwarz instance, it calls
01439     // setInnerPreconditioner() without first calling initialize().
01440 
01441     // Give the local matrix to the new inner solver.
01442     innerSolver->setMatrix (innerMatrix_);
01443 
01444     // If the user previously specified a parameter for the inner
01445     // preconditioner's type, then clear out that parameter and its
01446     // associated sublist.  Replace the inner preconditioner's type with
01447     // "CUSTOM", to make it obvious that AdditiveSchwarz's ParameterList
01448     // does not necessarily describe the current inner preconditioner.
01449     // We have to remove all allowed aliases of "inner preconditioner
01450     // name" before we may set it to "CUSTOM".  Users may also set this
01451     // parameter to "CUSTOM" themselves, but this is not required.
01452     removeInnerPrecName ();
01453     removeInnerPrecParams ();
01454     List_.set ("inner preconditioner name", "CUSTOM");
01455 
01456     // Bring the new inner solver's current status (initialized or
01457     // computed) in line with AdditiveSchwarz's current status.
01458     if (isInitialized ()) {
01459       innerPrec->initialize ();
01460     }
01461     if (isComputed ()) {
01462       innerPrec->compute ();
01463     }
01464   }
01465 
01466   // If the new inner solver is null, we don't change the initialized
01467   // or computed status of AdditiveSchwarz.  That way, AdditiveSchwarz
01468   // won't have to recompute innerMatrix_ if the inner solver changes.
01469   // This does introduce a new error condition in compute() and
01470   // apply(), but that's OK.
01471 
01472   // Set the new inner solver.
01473   Inverse_ = innerPrec;
01474 }
01475 
01476 template<class MatrixType, class LocalInverseType>
01477 void AdditiveSchwarz<MatrixType, LocalInverseType>::
01478 setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
01479 {
01480   // Don't set the matrix unless it is different from the current one.
01481   if (A.getRawPtr () != Matrix_.getRawPtr ()) {
01482     IsInitialized_ = false;
01483     IsComputed_ = false;
01484 
01485     // Reset all the state computed in initialize() and compute().
01486     Condest_ = -Teuchos::ScalarTraits<magnitude_type>::one ();
01487     OverlappingMatrix_ = Teuchos::null;
01488     ReorderedLocalizedMatrix_ = Teuchos::null;
01489     innerMatrix_ = Teuchos::null;
01490     SingletonMatrix_ = Teuchos::null;
01491     localMap_ = Teuchos::null;
01492     DistributedImporter_ = Teuchos::null;
01493 
01494     Matrix_ = A;
01495   }
01496 }
01497 
01498 } // namespace Ifpack2
01499 
01500 #define IFPACK2_ADDITIVESCHWARZ_INSTANT(S,LO,GO,N) \
01501   template class Ifpack2::AdditiveSchwarz< Tpetra::CrsMatrix<S, LO, GO, N> >; \
01502   template class Ifpack2::AdditiveSchwarz< Tpetra::CrsMatrix<S, LO, GO, N>, \
01503                                   Ifpack2::ILUT<Tpetra::CrsMatrix< S, LO, GO, N > > >;
01504 
01505 #endif // IFPACK2_ADDITIVESCHWARZ_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends