Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_BlockRelaxation_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_BLOCKRELAXATION_DEF_HPP
00044 #define IFPACK2_BLOCKRELAXATION_DEF_HPP
00045 
00046 #include "Ifpack2_BlockRelaxation_decl.hpp"
00047 #include "Ifpack2_LinearPartitioner.hpp"
00048 #include "Ifpack2_Details_UserPartitioner_decl.hpp"
00049 #include "Ifpack2_Details_UserPartitioner_def.hpp"
00050 #include <Ifpack2_Condest.hpp>
00051 #include <Ifpack2_Parameters.hpp>
00052 
00053 namespace Ifpack2 {
00054 
00055 template<class MatrixType,class ContainerType>
00056 void BlockRelaxation<MatrixType,ContainerType>::
00057 setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
00058 {
00059   if (A.getRawPtr () != A_.getRawPtr ()) { // it's a different matrix
00060     IsInitialized_ = false;
00061     IsComputed_ = false;
00062 
00063     Condest_ = 0;
00064 
00065     Partitioner_ = Teuchos::null;
00066     Importer_ = Teuchos::null;
00067     W_ = Teuchos::null;
00068 
00069     if (! A.is_null ()) {
00070       IsParallel_ = (A->getRowMap ()->getComm ()->getSize () > 1);
00071     }
00072 
00073     std::vector<Teuchos::RCP<ContainerType> > emptyVec;
00074     std::swap (Containers_, emptyVec);
00075     NumLocalBlocks_ = 0;
00076 
00077     A_ = A;
00078   }
00079 }
00080 
00081 template<class MatrixType,class ContainerType>
00082 BlockRelaxation<MatrixType,ContainerType>::
00083 BlockRelaxation (const Teuchos::RCP<const row_matrix_type>& A)
00084 : A_ (A),
00085   Time_ (Teuchos::rcp (new Teuchos::Time ("Ifpack2::BlockRelaxation"))),
00086   OverlapLevel_ (0),
00087   PartitionerType_ ("linear"),
00088   NumSweeps_ (1),
00089   NumLocalBlocks_(0),
00090   PrecType_ (Ifpack2::Details::JACOBI),
00091   MinDiagonalValue_ (STS::zero ()),
00092   DampingFactor_ (STS::one ()),
00093   IsParallel_ (false),
00094   ZeroStartingSolution_ (true),
00095   DoBackwardGS_ (false),
00096   Condest_ (STS::real(-STS::one ())),
00097   IsInitialized_ (false),
00098   IsComputed_ (false),
00099   NumInitialize_ (0),
00100   NumCompute_ (0),
00101   NumApply_ (0),
00102   InitializeTime_ (0.0),
00103   ComputeTime_ (0.0),
00104   ApplyTime_ (0.0),
00105   ComputeFlops_ (0.0),
00106   ApplyFlops_ (0.0),
00107   NumMyRows_ (0),
00108   NumGlobalRows_ (0),
00109   NumGlobalNonzeros_ (0)
00110 {
00111   TEUCHOS_TEST_FOR_EXCEPTION(
00112     A_.is_null (), std::invalid_argument,
00113     Teuchos::typeName(*this) << "::BlockRelaxation(): input matrix is null.");
00114 }
00115 
00116 //==========================================================================
00117 template<class MatrixType,class ContainerType>
00118 BlockRelaxation<MatrixType,ContainerType>::~BlockRelaxation() {}
00119 
00120 //==========================================================================
00121 template<class MatrixType,class ContainerType>
00122 void
00123 BlockRelaxation<MatrixType,ContainerType>::
00124 setParameters (const Teuchos::ParameterList& List)
00125 {
00126   Teuchos::ParameterList validparams;
00127   Ifpack2::getValidParameters (validparams);
00128   List.validateParameters (validparams);
00129 
00130   std::string PT;
00131   if (PrecType_ == Ifpack2::Details::JACOBI) {
00132     PT = "Jacobi";
00133   } else if (PrecType_ == Ifpack2::Details::GS) {
00134     PT = "Gauss-Seidel";
00135   } else if (PrecType_ == Ifpack2::Details::SGS) {
00136     PT = "Symmetric Gauss-Seidel";
00137   }
00138 
00139   Ifpack2::getParameter (List, "relaxation: type", PT);
00140 
00141   if (PT == "Jacobi") {
00142     PrecType_ = Ifpack2::Details::JACOBI;
00143   }
00144   else if (PT == "Gauss-Seidel") {
00145     PrecType_ = Ifpack2::Details::GS;
00146   }
00147   else if (PT == "Symmetric Gauss-Seidel") {
00148     PrecType_ = Ifpack2::Details::SGS;
00149   }
00150   else {
00151     TEUCHOS_TEST_FOR_EXCEPTION(
00152       true, std::invalid_argument, "Ifpack2::BlockRelaxation::setParameters: "
00153       "Invalid parameter value \"" << PT << "\" for parameter \"relaxation: type\".");
00154   }
00155 
00156   Ifpack2::getParameter (List, "relaxation: sweeps",NumSweeps_);
00157   Ifpack2::getParameter (List, "relaxation: damping factor", DampingFactor_);
00158   Ifpack2::getParameter (List, "relaxation: min diagonal value", MinDiagonalValue_);
00159   Ifpack2::getParameter (List, "relaxation: zero starting solution", ZeroStartingSolution_);
00160   Ifpack2::getParameter (List, "relaxation: backward mode",DoBackwardGS_);
00161   Ifpack2::getParameter (List, "partitioner: type",PartitionerType_);
00162   Ifpack2::getParameter (List, "partitioner: local parts",NumLocalBlocks_);
00163   Ifpack2::getParameter (List, "partitioner: overlap",OverlapLevel_);
00164 
00165   // check parameters
00166   if (PrecType_ != Ifpack2::Details::JACOBI) {
00167     OverlapLevel_ = 0;
00168   }
00169   if (NumLocalBlocks_ < 0) {
00170     NumLocalBlocks_ = A_->getNodeNumRows() / (-NumLocalBlocks_);
00171   }
00172   // other checks are performed in Partitioner_
00173 
00174   // NTS: Sanity check to be removed at a later date when Backward mode is enabled
00175   TEUCHOS_TEST_FOR_EXCEPTION(
00176     DoBackwardGS_, std::runtime_error,
00177     "Ifpack2::BlockRelaxation:setParameters: \"relaxation: backward mode\" == "
00178     "true is not supported yet.");
00179 
00180   // copy the list as each subblock's constructor will
00181   // require it later
00182   List_ = List;
00183 }
00184 
00185 //==========================================================================
00186 template<class MatrixType,class ContainerType>
00187 Teuchos::RCP<const Teuchos::Comm<int> >
00188 BlockRelaxation<MatrixType,ContainerType>::getComm() const{
00189   return A_->getComm();
00190 }
00191 
00192 //==========================================================================
00193 template<class MatrixType,class ContainerType>
00194 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
00195                                      typename MatrixType::local_ordinal_type,
00196                                      typename MatrixType::global_ordinal_type,
00197                                      typename MatrixType::node_type> >
00198 BlockRelaxation<MatrixType,ContainerType>::getMatrix() const {
00199   return(A_);
00200 }
00201 
00202 //==========================================================================
00203 template<class MatrixType,class ContainerType>
00204 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00205                                typename MatrixType::global_ordinal_type,
00206                                typename MatrixType::node_type> >
00207 BlockRelaxation<MatrixType,ContainerType>::getDomainMap() const {
00208   return A_->getDomainMap();
00209 }
00210 
00211 //==========================================================================
00212 template<class MatrixType,class ContainerType>
00213 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00214                                typename MatrixType::global_ordinal_type,
00215                                typename MatrixType::node_type> >
00216 BlockRelaxation<MatrixType,ContainerType>::getRangeMap() const {
00217   return A_->getRangeMap();
00218 }
00219 
00220 //==========================================================================
00221 template<class MatrixType,class ContainerType>
00222 bool BlockRelaxation<MatrixType,ContainerType>::hasTransposeApply() const {
00223   return true;
00224 }
00225 
00226 //==========================================================================
00227 template<class MatrixType,class ContainerType>
00228 int BlockRelaxation<MatrixType,ContainerType>::getNumInitialize() const {
00229   return(NumInitialize_);
00230 }
00231 
00232 //==========================================================================
00233 template<class MatrixType,class ContainerType>
00234 int BlockRelaxation<MatrixType,ContainerType>::getNumCompute() const {
00235   return(NumCompute_);
00236 }
00237 
00238 //==========================================================================
00239 template<class MatrixType,class ContainerType>
00240 int BlockRelaxation<MatrixType,ContainerType>::getNumApply() const {
00241   return(NumApply_);
00242 }
00243 
00244 //==========================================================================
00245 template<class MatrixType,class ContainerType>
00246 double BlockRelaxation<MatrixType,ContainerType>::getInitializeTime() const {
00247   return(InitializeTime_);
00248 }
00249 
00250 //==========================================================================
00251 template<class MatrixType,class ContainerType>
00252 double BlockRelaxation<MatrixType,ContainerType>::getComputeTime() const {
00253   return(ComputeTime_);
00254 }
00255 
00256 //==========================================================================
00257 template<class MatrixType,class ContainerType>
00258 double BlockRelaxation<MatrixType,ContainerType>::getApplyTime() const {
00259   return(ApplyTime_);
00260 }
00261 
00262 //==========================================================================
00263 template<class MatrixType,class ContainerType>
00264 double BlockRelaxation<MatrixType,ContainerType>::getComputeFlops() const {
00265   return(ComputeFlops_);
00266 }
00267 
00268 //==========================================================================
00269 template<class MatrixType,class ContainerType>
00270 double BlockRelaxation<MatrixType,ContainerType>::getApplyFlops() const {
00271   return ApplyFlops_;
00272 }
00273 
00274 //==========================================================================
00275 template<class MatrixType,class ContainerType>
00276 typename BlockRelaxation<MatrixType, ContainerType>::magnitude_type
00277 BlockRelaxation<MatrixType,ContainerType>::getCondEst () const
00278 {
00279   return Condest_;
00280 }
00281 
00282 //==========================================================================
00283 template<class MatrixType,class ContainerType>
00284 typename BlockRelaxation<MatrixType, ContainerType>::magnitude_type
00285 BlockRelaxation<MatrixType,ContainerType>::
00286 computeCondEst (CondestType CT,
00287                 typename MatrixType::local_ordinal_type MaxIters,
00288                 magnitude_type Tol,
00289                 const Teuchos::Ptr<const row_matrix_type>& matrix)
00290 {
00291   if (! isComputed ()) {// cannot compute right now
00292     return -STM::one ();
00293   }
00294 
00295   // always compute it. Call Condest() with no parameters to get
00296   // the previous estimate.
00297   Condest_ = Ifpack2::Condest (*this, CT, MaxIters, Tol, matrix);
00298 
00299   return Condest_;
00300 }
00301 
00302 //==========================================================================
00303 template<class MatrixType,class ContainerType>
00304 void
00305 BlockRelaxation<MatrixType,ContainerType>::
00306 apply (const Tpetra::MultiVector<typename MatrixType::scalar_type,
00307                                  typename MatrixType::local_ordinal_type,
00308                                  typename MatrixType::global_ordinal_type,
00309                                  typename MatrixType::node_type>& X,
00310        Tpetra::MultiVector<typename MatrixType::scalar_type,
00311                            typename MatrixType::local_ordinal_type,
00312                            typename MatrixType::global_ordinal_type,
00313                            typename MatrixType::node_type>& Y,
00314        Teuchos::ETransp mode,
00315        scalar_type alpha,
00316        scalar_type beta) const
00317 {
00318   TEUCHOS_TEST_FOR_EXCEPTION(
00319     ! isComputed (), std::runtime_error, "Ifpack2::BlockRelaxation::apply: "
00320     "isComputed() must be true prior to calling apply.");
00321   TEUCHOS_TEST_FOR_EXCEPTION(
00322     X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
00323     "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = "
00324     << X.getNumVectors () << " != Y.getNumVectors() = "
00325     << Y.getNumVectors () << ".");
00326   TEUCHOS_TEST_FOR_EXCEPTION(
00327     mode != Teuchos::NO_TRANS, std::logic_error, "Ifpack2::BlockRelaxation::"
00328     "apply: This method currently only implements the case mode == Teuchos::"
00329     "NO_TRANS.  Transposed modes are not currently supported.");
00330   TEUCHOS_TEST_FOR_EXCEPTION(
00331     alpha != Teuchos::ScalarTraits<scalar_type>::one (), std::logic_error,
00332     "Ifpack2::BlockRelaxation::apply: This method currently only implements "
00333     "the case alpha == 1.  You specified alpha = " << alpha << ".");
00334   TEUCHOS_TEST_FOR_EXCEPTION(
00335     beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
00336     "Ifpack2::BlockRelaxation::apply: This method currently only implements "
00337     "the case beta == 0.  You specified beta = " << beta << ".");
00338 
00339   Time_->start(true);
00340 
00341   // If X and Y are pointing to the same memory location,
00342   // we need to create an auxiliary vector, Xcopy
00343   Teuchos::RCP<const MV> X_copy;
00344   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues()) {
00345     X_copy = Teuchos::rcp (new MV (createCopy(X)));
00346   } else {
00347     X_copy = Teuchos::rcpFromRef (X);
00348   }
00349 
00350   if (ZeroStartingSolution_) {
00351     Y.putScalar (STS::zero ());
00352   }
00353 
00354   // Flops are updated in each of the following.
00355   switch (PrecType_) {
00356   case Ifpack2::Details::JACOBI:
00357     ApplyInverseJacobi(*X_copy,Y);
00358     break;
00359   case Ifpack2::Details::GS:
00360     ApplyInverseGS(*X_copy,Y);
00361     break;
00362   case Ifpack2::Details::SGS:
00363     ApplyInverseSGS(*X_copy,Y);
00364     break;
00365   default:
00366     throw std::runtime_error("Ifpack2::BlockRelaxation::apply internal logic error.");
00367   }
00368 
00369   ++NumApply_;
00370   Time_->stop();
00371   ApplyTime_ += Time_->totalElapsedTime();
00372 }
00373 
00374 //==========================================================================
00375 template<class MatrixType,class ContainerType>
00376 void BlockRelaxation<MatrixType,ContainerType>::applyMat(
00377           const Tpetra::MultiVector<typename MatrixType::scalar_type,
00378                                     typename MatrixType::local_ordinal_type,
00379                                     typename MatrixType::global_ordinal_type,
00380                                     typename MatrixType::node_type>& X,
00381                 Tpetra::MultiVector<typename MatrixType::scalar_type,
00382                                     typename MatrixType::local_ordinal_type,
00383                                     typename MatrixType::global_ordinal_type,
00384                                     typename MatrixType::node_type>& Y,
00385              Teuchos::ETransp mode) const
00386 {
00387   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == false, std::runtime_error,
00388      "Ifpack2::BlockRelaxation::applyMat() ERROR: isComputed() must be true prior to calling applyMat().");
00389   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00390      "Ifpack2::BlockRelaxation::applyMat() ERROR: X.getNumVectors() != Y.getNumVectors().");
00391   A_->apply (X, Y, mode);
00392 }
00393 
00394 //==========================================================================
00395 template<class MatrixType,class ContainerType>
00396 void BlockRelaxation<MatrixType,ContainerType>::initialize() {
00397   using Teuchos::rcp;
00398   typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
00399     row_graph_type;
00400 
00401   IsInitialized_ = false;
00402   Time_->start (true);
00403 
00404   NumMyRows_         = A_->getNodeNumRows ();
00405   NumGlobalRows_     = A_->getGlobalNumRows ();
00406   NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
00407 
00408   // NTS: Will need to add support for Zoltan2 partitions later Also,
00409   // will need a better way of handling the Graph typing issue.
00410   // Especially with ordinal types w.r.t the container.
00411 
00412   if (PartitionerType_ == "linear") {
00413     Partitioner_ =
00414       rcp (new Ifpack2::LinearPartitioner<row_graph_type> (A_->getGraph ()));
00415   } else if (PartitionerType_ == "user") {
00416     Partitioner_ =
00417       rcp (new Ifpack2::Details::UserPartitioner<row_graph_type> (A_->getGraph () ) );
00418   } else {
00419     // We should have checked for this in setParameters(), so it's a
00420     // logic_error, not an invalid_argument or runtime_error.
00421     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00422       "Ifpack2::BlockRelaxation::initialize, invalid partitioner type.");
00423   }
00424 
00425   // need to partition the graph of A
00426   Partitioner_->setParameters (List_);
00427   Partitioner_->compute ();
00428 
00429   // get actual number of partitions
00430   NumLocalBlocks_ = Partitioner_->numLocalParts ();
00431 
00432   // Note: Unlike Ifpack, we'll punt on computing W_ until compute(), which is where
00433   // we assume that the type of relaxation has been chosen.
00434 
00435   if (A_->getComm()->getSize() != 1) {
00436     IsParallel_ = true;
00437   } else {
00438     IsParallel_ = false;
00439   }
00440 
00441   ++NumInitialize_;
00442   Time_->stop ();
00443   InitializeTime_ += Time_->totalElapsedTime ();
00444   IsInitialized_ = true;
00445 }
00446 
00447 //==========================================================================
00448 template<class MatrixType,class ContainerType>
00449 void BlockRelaxation<MatrixType,ContainerType>::compute()
00450 {
00451   using Teuchos::rcp;
00452   typedef Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> vector_type;
00453   typedef Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
00454 
00455   // We should have checked for this in setParameters(), so it's a
00456   // logic_error, not an invalid_argument or runtime_error.
00457   TEUCHOS_TEST_FOR_EXCEPTION(NumSweeps_ < 0, std::logic_error,
00458     "Ifpack2::BlockRelaxation::compute, NumSweeps_ must be >= 0");
00459 
00460   if (! isInitialized ()) {
00461     initialize ();
00462   }
00463   Time_->start (true);
00464 
00465   // reset values
00466   IsComputed_ = false;
00467   Condest_ = STS::real(-STS::one ());
00468 
00469   // Extract the submatrices
00470   ExtractSubmatrices ();
00471 
00472   // Compute the weight vector if we're doing overlapped Jacobi (and
00473   // only if we're doing overlapped Jacobi).
00474   if (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0) {
00475     // weight of each vertex
00476     W_ = rcp (new vector_type (A_->getRowMap ()));
00477     W_->putScalar (STS::zero ());
00478     Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
00479 
00480     for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; ++i) {
00481       for (size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
00482         int LID = (*Partitioner_)(i,j);
00483         w_ptr[LID]+= STS::one();
00484       }
00485     }
00486     W_->reciprocal (*W_);
00487   }
00488 
00489   // We need to import data from external processors. Here I create a
00490   // Tpetra::Import object if needed (stealing from A_ if possible)
00491   // Marzio's comment:
00492   // Note that I am doing some strange stuff to set the components of Y
00493   // from Y2 (to save some time).
00494   //
00495   if (IsParallel_ && (PrecType_ == Ifpack2::Details::GS ||
00496                       PrecType_ == Ifpack2::Details::SGS)) {
00497     Importer_ = A_->getGraph ()->getImporter ();
00498     if (Importer_.is_null ()) {
00499       Importer_ = rcp (new import_type (A_->getDomainMap (), A_->getColMap ()));
00500     }
00501   }
00502 
00503   ++NumCompute_;
00504   Time_->stop ();
00505   ComputeTime_ += Time_->totalElapsedTime();
00506   IsComputed_ = true;
00507 }
00508 
00509 //==============================================================================
00510 template<class MatrixType,class ContainerType>
00511 void BlockRelaxation<MatrixType,ContainerType>::ExtractSubmatrices()
00512 {
00513   TEUCHOS_TEST_FOR_EXCEPTION(
00514     Partitioner_.is_null (), std::runtime_error,
00515     "Ifpack2::BlockRelaxation::ExtractSubmatrices: Partitioner object is null.");
00516 
00517   NumLocalBlocks_ = Partitioner_->numLocalParts ();
00518   Containers_.resize (NumLocalBlocks_);
00519 
00520   for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
00521     const size_t numRows = Partitioner_->numRowsInPart (i);
00522 
00523     // Extract a list of the indices of each partitioner row.
00524     Teuchos::Array<local_ordinal_type> localRows (numRows);
00525     for (size_t j = 0; j < numRows; ++j) {
00526       localRows[j] = (*Partitioner_) (i,j);
00527     }
00528 
00529     Containers_[i] = Teuchos::rcp (new ContainerType (A_, localRows ()));
00530     Containers_[i]->setParameters (List_);
00531     Containers_[i]->initialize ();
00532     Containers_[i]->compute ();
00533   }
00534 }
00535 
00536 
00537 
00538 //==========================================================================
00539 template<class MatrixType,class ContainerType>
00540 void BlockRelaxation<MatrixType,ContainerType>::ApplyInverseJacobi (const MV& X, MV& Y) const
00541 {
00542   const size_t NumVectors = X.getNumVectors ();
00543   MV AY (Y.getMap (), NumVectors);
00544 
00545   // Initial matvec not needed
00546   int starting_iteration = 0;
00547   if (ZeroStartingSolution_) {
00548     DoJacobi(X,Y);
00549     starting_iteration = 1;
00550   }
00551 
00552   const scalar_type ONE = STS::one ();
00553   for (int j = starting_iteration; j < NumSweeps_; ++j) {
00554     applyMat (Y, AY);
00555     AY.update (ONE, X, -ONE);
00556     DoJacobi (AY, Y);
00557 
00558     // Flops for matrix apply & update
00559     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
00560   }
00561 
00562 }
00563 
00564 
00565 //==========================================================================
00566 template<class MatrixType,class ContainerType>
00567 void BlockRelaxation<MatrixType,ContainerType>::DoJacobi(const MV& X, MV& Y) const
00568 {
00569   const size_t NumVectors = X.getNumVectors();
00570   const scalar_type one = STS::one();
00571   // Note: Flop counts copied naively from Ifpack.
00572 
00573   if (OverlapLevel_ == 0) {
00574     // Non-overlapping Jacobi
00575     for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
00576       // may happen that a partition is empty
00577       if (Containers_[i]->getNumRows () == 0) {
00578         continue;
00579       }
00580       Containers_[i]->apply (X, Y, Teuchos::NO_TRANS, DampingFactor_, one);
00581       ApplyFlops_ += NumVectors * 2 * NumGlobalRows_;
00582     }
00583   }
00584   else {
00585     // Overlapping Jacobi
00586     for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; i++) {
00587       // may happen that a partition is empty
00588       if (Containers_[i]->getNumRows() == 0) continue;
00589 
00590       try {
00591         Containers_[i]->weightedApply(X,Y,*W_,Teuchos::NO_TRANS,DampingFactor_,one);
00592       } catch (std::exception& e) {
00593         std::cerr << "BlockRelaxation::DoJacobi: Containers_[" << i
00594                   << "]->weightedApply() threw an exception: " << e.what ()
00595                   << std::endl;
00596         throw;
00597       }
00598       // NOTE: do not count (for simplicity) the flops due to overlapping rows
00599       ApplyFlops_ += NumVectors * 4 * NumGlobalRows_;
00600     }
00601   }
00602 }
00603 
00604 //==========================================================================
00605 template<class MatrixType,class ContainerType>
00606 void BlockRelaxation<MatrixType,ContainerType>::
00607 ApplyInverseGS (const MV& X, MV& Y) const
00608 {
00609   MV Xcopy (X);
00610   for (int j = 0; j < NumSweeps_ ; j++) {
00611     DoGaussSeidel (Xcopy, Y);
00612     if (j != NumSweeps_ - 1) {
00613       Xcopy = X;
00614     }
00615   }
00616 }
00617 
00618 
00619 //==============================================================================
00620 template<class MatrixType,class ContainerType>
00621 void BlockRelaxation<MatrixType,ContainerType>::
00622 DoGaussSeidel (MV& X, MV& Y) const
00623 {
00624   using Teuchos::Array;
00625   using Teuchos::ArrayRCP;
00626   using Teuchos::ArrayView;
00627   using Teuchos::RCP;
00628   using Teuchos::rcp;
00629   using Teuchos::rcpFromRef;
00630 
00631   // Note: Flop counts copied naively from Ifpack.
00632 
00633   const scalar_type    one =  STS::one ();
00634   int Length = A_->getNodeMaxNumRowEntries();
00635   const size_t NumVectors = X.getNumVectors();
00636   Array<scalar_type>         Values;
00637   Array<local_ordinal_type>   Indices;
00638   Values.resize(Length);
00639   Indices.resize(Length);
00640 
00641   // an additonal vector is needed by parallel computations
00642   // (note that applications through Ifpack2_AdditiveSchwarz
00643   // are always seen are serial)
00644   RCP<MV> Y2;
00645   if (IsParallel_) {
00646     Y2 = rcp (new MV (Importer_->getTargetMap(), NumVectors));
00647   } else {
00648     Y2 = rcpFromRef (Y);
00649   }
00650 
00651   // I think I decided I need two extra vectors:
00652   // One to store the sum of the corrections (initialized to zero)
00653   // One to store the temporary residual (doesn't matter if it is zeroed or not)
00654   // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
00655   MV Residual(X.getMap(),NumVectors,false);
00656 
00657   ArrayRCP<ArrayRCP<scalar_type> >           x_ptr = X.get2dViewNonConst();
00658   ArrayRCP<ArrayRCP<scalar_type> >           y_ptr = Y.get2dViewNonConst();
00659   ArrayRCP<ArrayRCP<scalar_type> >          y2_ptr = Y2->get2dViewNonConst();
00660   ArrayRCP<ArrayRCP<scalar_type> >    residual_ptr = Residual.get2dViewNonConst();
00661 
00662   // data exchange is here, once per sweep
00663   if (IsParallel_)  Y2->doImport(Y,*Importer_,Tpetra::INSERT);
00664 
00665   for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
00666     if (Containers_[i]->getNumRows () == 0) {
00667       continue; // Skip empty partitions
00668     }
00669 
00670     // update from previous block
00671     ArrayView<const local_ordinal_type> localRows =
00672       Containers_[i]->getLocalRows ();
00673     const size_t localNumRows = Containers_[i]->getNumRows ();
00674     for (size_t j = 0; j < localNumRows; ++j) {
00675       const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
00676       size_t NumEntries;
00677       A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
00678 
00679       for (size_t m = 0; m < NumVectors; ++m) {
00680         ArrayView<const scalar_type> x_local = (x_ptr())[m]();
00681         ArrayView<scalar_type>      y2_local = (y2_ptr())[m]();
00682         ArrayView<scalar_type>       r_local = (residual_ptr())[m]();
00683 
00684         r_local[LID] = x_local[LID];
00685         for (size_t k = 0; k < NumEntries; ++k) {
00686           const local_ordinal_type col = Indices[k];
00687           r_local[LID] -= Values[k] * y2_local[col];
00688         }
00689       }
00690     }
00691     // solve with this block
00692     //
00693     // Note: I'm abusing the ordering information, knowing that X/Y
00694     // and Y2 have the same ordering for on-proc unknowns.
00695     //
00696     // Note: Add flop counts for inverse apply
00697     Containers_[i]->apply (Residual, *Y2, Teuchos::NO_TRANS,
00698                            DampingFactor_,one);
00699 
00700     // operations for all getrow's
00701     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
00702   } // end for NumLocalBlocks_
00703 
00704   // Attention: this is delicate... Not all combinations
00705   // of Y2 and Y will always work (tough for ML it should be ok)
00706   if (IsParallel_) {
00707     for (size_t m = 0; m < NumVectors; ++m) {
00708       ArrayView<scalar_type>      y2_local = (y2_ptr())[m]();
00709       ArrayView<scalar_type>      y_local  = (y_ptr())[m]();
00710       for (size_t i = 0; i < NumMyRows_; ++i) {
00711         y_local[i] = y2_local[i];
00712       }
00713     }
00714   }
00715 }
00716 
00717 //==========================================================================
00718 template<class MatrixType,class ContainerType>
00719 void
00720 BlockRelaxation<MatrixType,ContainerType>::
00721 ApplyInverseSGS (const MV& X, MV& Y) const
00722 {
00723   MV Xcopy (X);
00724   for (int j = 0; j < NumSweeps_; ++j) {
00725     DoSGS (Xcopy, Y);
00726     if (j != NumSweeps_ - 1) {
00727       Xcopy = X;
00728     }
00729   }
00730 }
00731 
00732 //==========================================================================
00733 template<class MatrixType,class ContainerType>
00734 void
00735 BlockRelaxation<MatrixType,ContainerType>::DoSGS (MV& X, MV& Y) const
00736 {
00737   using Teuchos::Array;
00738   using Teuchos::ArrayRCP;
00739   using Teuchos::ArrayView;
00740   using Teuchos::RCP;
00741   using Teuchos::rcp;
00742   using Teuchos::rcpFromRef;
00743 
00744   const scalar_type    one =  STS::one ();
00745   int Length = A_->getNodeMaxNumRowEntries();
00746   const size_t NumVectors = X.getNumVectors();
00747   Array<scalar_type>         Values;
00748   Array<local_ordinal_type>   Indices;
00749   Values.resize(Length);
00750   Indices.resize(Length);
00751 
00752   // an additonal vector is needed by parallel computations
00753   // (note that applications through Ifpack2_AdditiveSchwarz
00754   // are always seen are serial)
00755   RCP<MV> Y2;
00756   if (IsParallel_) {
00757     Y2 = rcp (new MV (Importer_->getTargetMap (), NumVectors));
00758   } else {
00759     Y2 = rcpFromRef (Y);
00760   }
00761 
00762   // I think I decided I need two extra vectors:
00763   // One to store the sum of the corrections (initialized to zero)
00764   // One to store the temporary residual (doesn't matter if it is zeroed or not)
00765   // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
00766   MV Residual(X.getMap(),NumVectors,false);
00767 
00768   ArrayRCP<ArrayRCP<scalar_type> >     x_ptr       = X.get2dViewNonConst();
00769   ArrayRCP<ArrayRCP<scalar_type> >     y_ptr       = Y.get2dViewNonConst();
00770   ArrayRCP<ArrayRCP<scalar_type> >     y2_ptr      = Y2->get2dViewNonConst();
00771   ArrayRCP<ArrayRCP<scalar_type> >    residual_ptr = Residual.get2dViewNonConst();
00772 
00773   // data exchange is here, once per sweep
00774   if (IsParallel_) {
00775     Y2->doImport (Y, *Importer_, Tpetra::INSERT);
00776   }
00777 
00778   // Forward Sweep
00779   for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
00780     if (Containers_[i]->getNumRows () == 0) {
00781       continue; // Skip empty partitions
00782     }
00783     // update from previous block
00784     ArrayView<const local_ordinal_type> localRows =
00785       Containers_[i]->getLocalRows ();
00786     for (size_t j = 0; j < Containers_[i]->getNumRows (); ++j) {
00787       const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
00788       size_t NumEntries;
00789       A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
00790 
00791       //set tmpresid = initresid - A*correction
00792       for (size_t m = 0; m < NumVectors; ++m) {
00793         ArrayView<const scalar_type> x_local = (x_ptr())[m]();
00794         ArrayView<scalar_type>      y2_local = (y2_ptr())[m]();
00795         ArrayView<scalar_type>       r_local = (residual_ptr())[m]();
00796 
00797         r_local[LID] = x_local[LID];
00798         for (size_t k = 0 ; k < NumEntries ; k++) {
00799           local_ordinal_type col = Indices[k];
00800           r_local[LID] -= Values[k] * y2_local[col];
00801         }
00802       }
00803     }
00804     // solve with this block
00805     //
00806     // Note: I'm abusing the ordering information, knowing that X/Y
00807     // and Y2 have the same ordering for on-proc unknowns.
00808     //
00809     // Note: Add flop counts for inverse apply
00810     Containers_[i]->apply (Residual, *Y2, Teuchos::NO_TRANS,
00811                            DampingFactor_, one);
00812 
00813     // operations for all getrow's
00814     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
00815   }// end forward sweep
00816 
00817   // Reverse Sweep
00818   //
00819   // mfh 12 July 2013: The unusual iteration bounds, and the use of
00820   // i-1 rather than i in the loop body, ensure correctness even if
00821   // local_ordinal_type is unsigned.  "i = NumLocalBlocks_-1; i >= 0;
00822   // i--" will loop forever if local_ordinal_type is unsigned, because
00823   // unsigned integers are (trivially) always nonnegative.
00824   for (local_ordinal_type i = NumLocalBlocks_; i > 0; --i) {
00825     if (Containers_[i-1]->getNumRows () == 0) {
00826       continue; // Skip empty partitions
00827     }
00828     // update from previous block
00829     ArrayView<const local_ordinal_type> localRows =
00830       Containers_[i-1]->getLocalRows ();
00831     for (size_t j = 0; j < Containers_[i-1]->getNumRows (); ++j) {
00832       const local_ordinal_type LID = localRows[j]; // Containers_[i-1]->ID (j);
00833       size_t NumEntries;
00834       A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
00835 
00836       //set tmpresid = initresid - A*correction
00837       for (size_t m = 0; m < NumVectors; ++m) {
00838         ArrayView<const scalar_type> x_local = (x_ptr())[m]();
00839         ArrayView<scalar_type>      y2_local = (y2_ptr())[m]();
00840         ArrayView<scalar_type>       r_local = (residual_ptr())[m]();
00841 
00842         r_local [LID] = x_local[LID];
00843         for (size_t k = 0; k < NumEntries; ++k)  {
00844           local_ordinal_type col = Indices[k];
00845           r_local[LID] -= Values[k] * y2_local[col];
00846         }
00847       }
00848     }
00849 
00850     // solve with this block
00851     //
00852     // Note: I'm abusing the ordering information, knowing that X/Y
00853     // and Y2 have the same ordering for on-proc unknowns.
00854     //
00855     // Note: Add flop counts for inverse apply
00856     Containers_[i-1]->apply (Residual, *Y2, Teuchos::NO_TRANS,
00857                            DampingFactor_, one);
00858 
00859     // operations for all getrow's
00860     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
00861   } //end reverse sweep
00862 
00863   // Attention: this is delicate... Not all combinations
00864   // of Y2 and Y will always work (though for ML it should be ok)
00865   if (IsParallel_) {
00866     for (size_t m = 0; m < NumVectors; ++m) {
00867       ArrayView<scalar_type>       y_local = (y_ptr())[m]();
00868       ArrayView<scalar_type>      y2_local = (y2_ptr())[m]();
00869       for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00870         y_local[i] = y2_local[i];
00871       }
00872     }
00873   }
00874 }
00875 
00876 template<class MatrixType, class ContainerType>
00877 std::string BlockRelaxation<MatrixType,ContainerType>::description () const
00878 {
00879   std::ostringstream out;
00880 
00881   // Output is a valid YAML dictionary in flow style.  If you don't
00882   // like everything on a single line, you should call describe()
00883   // instead.
00884   out << "\"Ifpack2::BlockRelaxation\": {";
00885   if (this->getObjectLabel () != "") {
00886     out << "Label: \"" << this->getObjectLabel () << "\", ";
00887   }
00888   out << "Initialized: " << (isInitialized () ? "true" : "false") << ", ";
00889   out << "Computed: " << (isComputed () ? "true" : "false") << ", ";
00890 
00891   if (A_.is_null ()) {
00892     out << "Matrix: null";
00893   }
00894   else {
00895     out << "Matrix: not null"
00896         << ", Global matrix dimensions: ["
00897         << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]";
00898   }
00899 
00900   // It's useful to print this instance's relaxation method.  If you
00901   // want more info than that, call describe() instead.
00902   out << "\"relaxation: type\": ";
00903   if (PrecType_ == Ifpack2::Details::JACOBI) {
00904     out << "Block Jacobi";
00905   } else if (PrecType_ == Ifpack2::Details::GS) {
00906     out << "Block Gauss-Seidel";
00907   } else if (PrecType_ == Ifpack2::Details::SGS) {
00908     out << "Block Symmetric Gauss-Seidel";
00909   } else {
00910     out << "INVALID";
00911   }
00912 
00913   out << "}";
00914   return out.str();
00915 }
00916 
00917 
00918 template<class MatrixType,class ContainerType>
00919 void BlockRelaxation<MatrixType,ContainerType>::
00920 describe (Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
00921 {
00922   using std::endl;
00923   using std::setw;
00924   using Teuchos::TypeNameTraits;
00925   using Teuchos::VERB_DEFAULT;
00926   using Teuchos::VERB_NONE;
00927   using Teuchos::VERB_LOW;
00928   using Teuchos::VERB_MEDIUM;
00929   using Teuchos::VERB_HIGH;
00930   using Teuchos::VERB_EXTREME;
00931 
00932   Teuchos::EVerbosityLevel vl = verbLevel;
00933   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00934   const int myImageID = A_->getComm()->getRank();
00935 
00936   // Convention requires that describe() begin with a tab.
00937   Teuchos::OSTab tab (out);
00938 
00939   //    none: print nothing
00940   //     low: print O(1) info from node 0
00941   //  medium:
00942   //    high:
00943   // extreme:
00944   if (vl != VERB_NONE && myImageID == 0) {
00945     out << "Ifpack2::BlockRelaxation<"
00946         << TypeNameTraits<MatrixType>::name () << ", "
00947         << TypeNameTraits<ContainerType>::name () << " >:";
00948     Teuchos::OSTab tab1 (out);
00949 
00950     if (this->getObjectLabel () != "") {
00951       out << "label: \"" << this->getObjectLabel () << "\"" << endl;
00952     }
00953     out << "initialized: " << (isInitialized () ? "true" : "false") << endl
00954         << "computed: " << (isComputed () ? "true" : "false") << endl;
00955 
00956     std::string precType;
00957     if (PrecType_ == Ifpack2::Details::JACOBI) {
00958       precType = "Block Jacobi";
00959     } else if (PrecType_ == Ifpack2::Details::GS) {
00960       precType = "Block Gauss-Seidel";
00961     } else if (PrecType_ == Ifpack2::Details::SGS) {
00962       precType = "Block Symmetric Gauss-Seidel";
00963     }
00964     out << "type: " << precType << endl
00965         << "global number of rows: " << A_->getGlobalNumRows () << endl
00966         << "global number of columns: " << A_->getGlobalNumCols () << endl
00967         << "number of sweeps: " << NumSweeps_ << endl
00968         << "damping factor: " << DampingFactor_ << endl
00969         << "backwards mode: "
00970         << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ? "true" : "false")
00971         << endl
00972         << "zero starting solution: "
00973         << (ZeroStartingSolution_ ? "true" : "false") << endl
00974         << "condition number estimate: " << Condest_ << endl;
00975 
00976     out << "===============================================================================" << endl;
00977     out << "Phase           # calls    Total Time (s)     Total MFlops      MFlops/s       " << endl;
00978     out << "------------    -------    ---------------    ---------------   ---------------" << endl;
00979     out << setw(12) << "initialize()" << setw(5) << getNumInitialize() << "    " << setw(15) << getInitializeTime() << endl;
00980     out << setw(12) << "compute()" << setw(5) << getNumCompute()    << "    " << setw(15) << getComputeTime() << "    "
00981         << setw(15) << getComputeFlops() << "    "
00982         << setw(15) << (getComputeTime() != 0.0 ? getComputeFlops() / getComputeTime() * 1.0e-6 : 0.0) << endl;
00983     out << setw(12) << "apply()" << setw(5) << getNumApply()    << "    " << setw(15) << getApplyTime() << "    "
00984         << setw(15) << getApplyFlops() << "    "
00985         << setw(15) << (getApplyTime() != 0.0 ? getApplyFlops() / getApplyTime() * 1.0e-6 : 0.0) << endl;
00986     out << "===============================================================================" << endl;
00987     out << endl;
00988   }
00989 }
00990 
00991 }//namespace Ifpack2
00992 
00993 
00994 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION
00995 
00996 // For ETI
00997 #include "Ifpack2_DenseContainer_decl.hpp"
00998 #include "Ifpack2_SparseContainer_decl.hpp"
00999 #include "Ifpack2_ILUT_decl.hpp"
01000 
01001 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \
01002   template \
01003   class Ifpack2::BlockRelaxation<      \
01004     Tpetra::CrsMatrix<S, LO, GO, N>, \
01005     Ifpack2::SparseContainer<       \
01006       Tpetra::CrsMatrix<S, LO, GO, N>, \
01007       Ifpack2::ILUT< ::Tpetra::CrsMatrix<S,LO,GO,N> > > >; \
01008   template \
01009   class Ifpack2::BlockRelaxation<      \
01010     Tpetra::CrsMatrix<S, LO, GO, N>, \
01011     Ifpack2::DenseContainer<        \
01012       Tpetra::CrsMatrix<S, LO, GO, N>, \
01013       S > >;
01014 
01015 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION
01016 
01017 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends