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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #ifndef IFPACK2_BLOCKRELAXATION_DEF_HPP
00031 #define IFPACK2_BLOCKRELAXATION_DEF_HPP
00032 
00033 #include "Ifpack2_BlockRelaxation_decl.hpp"
00034 #include "Ifpack2_LinearPartitioner_decl.hpp"
00035 
00036 namespace Ifpack2 {
00037 
00038 //==========================================================================
00039 template<class MatrixType,class ContainerType>
00040 BlockRelaxation<MatrixType,ContainerType>::BlockRelaxation(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A)
00041 : A_(A),
00042   Time_( Teuchos::rcp( new Teuchos::Time("Ifpack2::BlockRelaxation") ) ),
00043   OverlapLevel_(0),
00044   PartitionerType_("linear"),
00045   NumSweeps_(1),
00046   PrecType_(Ifpack2::JACOBI),
00047   MinDiagonalValue_(0.0),
00048   DampingFactor_(1.0),
00049   IsParallel_(false),
00050   ZeroStartingSolution_(true),
00051   DoBackwardGS_(false),
00052   Condest_(-1.0),
00053   IsInitialized_(false),
00054   IsComputed_(false),
00055   NumInitialize_(0),
00056   NumCompute_(0),
00057   NumApply_(0),
00058   InitializeTime_(0.0),
00059   ComputeTime_(0.0),
00060   ApplyTime_(0.0),
00061   ComputeFlops_(0.0),
00062   ApplyFlops_(0.0),
00063   NumMyRows_(0),
00064   NumGlobalRows_(0),
00065   NumGlobalNonzeros_(0)
00066 {
00067   TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error, 
00068       Teuchos::typeName(*this) << "::BlockRelaxation(): input matrix reference was null.");
00069 }
00070 
00071 //==========================================================================
00072 template<class MatrixType,class ContainerType>
00073 BlockRelaxation<MatrixType,ContainerType>::~BlockRelaxation() {
00074 }
00075 
00076 //==========================================================================
00077 template<class MatrixType,class ContainerType>
00078 void BlockRelaxation<MatrixType,ContainerType>::setParameters(const Teuchos::ParameterList& List)
00079 {
00080   Teuchos::ParameterList validparams;
00081   Ifpack2::getValidParameters(validparams);
00082   List.validateParameters(validparams);
00083 
00084   std::string PT;
00085   if (PrecType_ == Ifpack2::JACOBI)
00086     PT = "Jacobi";
00087   else if (PrecType_ == Ifpack2::GS)
00088     PT = "Gauss-Seidel";
00089   else if (PrecType_ == Ifpack2::SGS)
00090     PT = "Symmetric Gauss-Seidel";
00091 
00092   Ifpack2::getParameter(List, "relaxation: type", PT);
00093 
00094   if (PT == "Jacobi")
00095     PrecType_ = Ifpack2::JACOBI;
00096   else if (PT == "Gauss-Seidel")
00097     PrecType_ = Ifpack2::GS;
00098   else if (PT == "Symmetric Gauss-Seidel")
00099     PrecType_ = Ifpack2::SGS;
00100   else {
00101     std::ostringstream osstr;
00102     osstr << "Ifpack2::BlockRelaxation::setParameters: unsupported parameter-value for 'relaxation: type' (" << PT << ")";
00103     std::string str = osstr.str();
00104     throw std::runtime_error(str);
00105   }
00106 
00107   Ifpack2::getParameter(List, "relaxation: sweeps",NumSweeps_);
00108   Ifpack2::getParameter(List, "relaxation: damping factor", DampingFactor_);
00109   Ifpack2::getParameter(List, "relaxation: min diagonal value", MinDiagonalValue_);
00110   Ifpack2::getParameter(List, "relaxation: zero starting solution", ZeroStartingSolution_);
00111   Ifpack2::getParameter(List, "relaxation: backward mode",DoBackwardGS_);
00112   Ifpack2::getParameter(List, "partitioner: type",PartitionerType_);
00113   Ifpack2::getParameter(List, "partitioner: local parts",NumLocalBlocks_);
00114   Ifpack2::getParameter(List, "partitioner: overlap",OverlapLevel_);
00115 
00116   // check parameters
00117   if (PrecType_ != Ifpack2::JACOBI)
00118     OverlapLevel_ = 0;
00119   if (NumLocalBlocks_ < 0)
00120     NumLocalBlocks_ = A_->getNodeNumRows() / (-NumLocalBlocks_);
00121   // other checks are performed in Partitioner_
00122   
00123 
00124   // NTS: Sanity check to be removed at a later date when Backward mode is enabled
00125   TEUCHOS_TEST_FOR_EXCEPTION(DoBackwardGS_ == true, std::runtime_error,
00126            "Ifpack2::BlockRelaxation:setParameters ERROR: 'relaxation: backward mode' == true is not supported yet.");
00127 
00128 
00129   // copy the list as each subblock's constructor will
00130   // require it later
00131   List_ = List;
00132 }
00133 
00134 //==========================================================================
00135 template<class MatrixType,class ContainerType>
00136 const Teuchos::RCP<const Teuchos::Comm<int> > & 
00137 BlockRelaxation<MatrixType,ContainerType>::getComm() const{
00138   return A_->getComm();
00139 }
00140 
00141 //==========================================================================
00142 template<class MatrixType,class ContainerType>
00143 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
00144                                      typename MatrixType::local_ordinal_type,
00145                                      typename MatrixType::global_ordinal_type,
00146                                      typename MatrixType::node_type> >
00147 BlockRelaxation<MatrixType,ContainerType>::getMatrix() const {
00148   return(A_);
00149 }
00150 
00151 //==========================================================================
00152 template<class MatrixType,class ContainerType>
00153 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00154                                      typename MatrixType::global_ordinal_type,
00155                                      typename MatrixType::node_type> >&
00156 BlockRelaxation<MatrixType,ContainerType>::getDomainMap() const {
00157   return A_->getDomainMap();
00158 }
00159 
00160 //==========================================================================
00161 template<class MatrixType,class ContainerType>
00162 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00163                                      typename MatrixType::global_ordinal_type,
00164                                      typename MatrixType::node_type> >&
00165 BlockRelaxation<MatrixType,ContainerType>::getRangeMap() const {
00166   return A_->getRangeMap();
00167 }
00168 
00169 //==========================================================================
00170 template<class MatrixType,class ContainerType>
00171 bool BlockRelaxation<MatrixType,ContainerType>::hasTransposeApply() const {
00172   return true;
00173 }
00174 
00175 //==========================================================================
00176 template<class MatrixType,class ContainerType>
00177 int BlockRelaxation<MatrixType,ContainerType>::getNumInitialize() const {
00178   return(NumInitialize_);
00179 }
00180 
00181 //==========================================================================
00182 template<class MatrixType,class ContainerType>
00183 int BlockRelaxation<MatrixType,ContainerType>::getNumCompute() const {
00184   return(NumCompute_);
00185 }
00186 
00187 //==========================================================================
00188 template<class MatrixType,class ContainerType>
00189 int BlockRelaxation<MatrixType,ContainerType>::getNumApply() const {
00190   return(NumApply_);
00191 }
00192 
00193 //==========================================================================
00194 template<class MatrixType,class ContainerType>
00195 double BlockRelaxation<MatrixType,ContainerType>::getInitializeTime() const {
00196   return(InitializeTime_);
00197 }
00198 
00199 //==========================================================================
00200 template<class MatrixType,class ContainerType>
00201 double BlockRelaxation<MatrixType,ContainerType>::getComputeTime() const {
00202   return(ComputeTime_);
00203 }
00204 
00205 //==========================================================================
00206 template<class MatrixType,class ContainerType>
00207 double BlockRelaxation<MatrixType,ContainerType>::getApplyTime() const {
00208   return(ApplyTime_);
00209 }
00210 
00211 //==========================================================================
00212 template<class MatrixType,class ContainerType>
00213 double BlockRelaxation<MatrixType,ContainerType>::getComputeFlops() const {
00214   return(ComputeFlops_);
00215 }
00216 
00217 //==========================================================================
00218 template<class MatrixType,class ContainerType>
00219 double BlockRelaxation<MatrixType,ContainerType>::getApplyFlops() const {
00220   return(ApplyFlops_);
00221 }
00222 
00223 //==========================================================================
00224 template<class MatrixType,class ContainerType>
00225 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00226 BlockRelaxation<MatrixType,ContainerType>::getCondEst() const
00227 {
00228   return(Condest_);
00229 }
00230 
00231 //==========================================================================
00232 template<class MatrixType,class ContainerType>
00233 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00234 BlockRelaxation<MatrixType,ContainerType>::computeCondEst(
00235                      CondestType CT,
00236                      typename MatrixType::local_ordinal_type MaxIters, 
00237                      magnitudeType Tol,
00238      const Teuchos::Ptr<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
00239                                                 typename MatrixType::local_ordinal_type,
00240                                                 typename MatrixType::global_ordinal_type,
00241                                                 typename MatrixType::node_type> > &matrix)
00242 {
00243   if (!isComputed()) // cannot compute right now
00244     return(-1.0);
00245 
00246   // always compute it. Call Condest() with no parameters to get
00247   // the previous estimate.
00248   Condest_ = Ifpack2::Condest(*this, CT, MaxIters, Tol, matrix);
00249 
00250   return(Condest_);
00251 }
00252 
00253 //==========================================================================
00254 template<class MatrixType,class ContainerType>
00255 void BlockRelaxation<MatrixType,ContainerType>::apply(
00256           const Tpetra::MultiVector<typename MatrixType::scalar_type,
00257                                     typename MatrixType::local_ordinal_type,
00258                                     typename MatrixType::global_ordinal_type,
00259                                     typename MatrixType::node_type>& X,
00260                 Tpetra::MultiVector<typename MatrixType::scalar_type,
00261                                     typename MatrixType::local_ordinal_type,
00262                                     typename MatrixType::global_ordinal_type,
00263                                     typename MatrixType::node_type>& Y,
00264                 Teuchos::ETransp mode,
00265                  Scalar alpha,
00266                  Scalar beta) const
00267 {
00268   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == false, std::runtime_error,
00269      "Ifpack2::BlockRelaxation::apply ERROR: isComputed() must be true prior to calling apply.");
00270 
00271   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00272      "Ifpack2::BlockRelaxation::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
00273 
00274   TEUCHOS_TEST_FOR_EXCEPTION(mode != Teuchos::NO_TRANS, std::runtime_error,
00275            "Ifpack2::BlockRelaxation::apply ERORR: transpose modes not supported.");
00276 
00277   Time_->start(true);
00278 
00279   // If X and Y are pointing to the same memory location,
00280   // we need to create an auxiliary vector, Xcopy
00281   Teuchos::RCP< const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xcopy;
00282   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues())
00283     Xcopy = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X) );
00284   else
00285     Xcopy = Teuchos::rcp( &X, false );
00286 
00287   if (ZeroStartingSolution_)
00288     Y.putScalar(0.0);
00289 
00290   // Flops are updated in each of the following.
00291   switch (PrecType_) {
00292   case Ifpack2::JACOBI:
00293     ApplyInverseJacobi(*Xcopy,Y);
00294     break;
00295   case Ifpack2::GS:
00296     ApplyInverseGS(*Xcopy,Y);
00297     break;
00298   case Ifpack2::SGS:
00299     ApplyInverseSGS(*Xcopy,Y);
00300     break;
00301   default:
00302     throw std::runtime_error("Ifpack2::BlockRelaxation::apply internal logic error.");
00303   }
00304 
00305   ++NumApply_;
00306   Time_->stop();
00307   ApplyTime_ += Time_->totalElapsedTime();
00308 }
00309 
00310 //==========================================================================
00311 template<class MatrixType,class ContainerType>
00312 void BlockRelaxation<MatrixType,ContainerType>::applyMat(
00313           const Tpetra::MultiVector<typename MatrixType::scalar_type,
00314                                     typename MatrixType::local_ordinal_type,
00315                                     typename MatrixType::global_ordinal_type,
00316                                     typename MatrixType::node_type>& X,
00317                 Tpetra::MultiVector<typename MatrixType::scalar_type,
00318                                     typename MatrixType::local_ordinal_type,
00319                                     typename MatrixType::global_ordinal_type,
00320                                     typename MatrixType::node_type>& Y,
00321              Teuchos::ETransp mode) const
00322 {
00323   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == false, std::runtime_error,
00324      "Ifpack2::BlockRelaxation::applyMat() ERROR: isComputed() must be true prior to calling applyMat().");
00325   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00326      "Ifpack2::BlockRelaxation::applyMat() ERROR: X.getNumVectors() != Y.getNumVectors().");
00327   A_->apply(X, Y, mode);
00328 }
00329 
00330 //==========================================================================
00331 template<class MatrixType,class ContainerType>
00332 void BlockRelaxation<MatrixType,ContainerType>::initialize() {
00333   IsInitialized_ = false;
00334 
00335   TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error,
00336     "Ifpack2::BlockRelaxation::Initialize ERROR, Matrix is NULL");
00337 
00338   Time_->start(true);
00339 
00340   NumMyRows_         = A_->getNodeNumRows();
00341   NumGlobalRows_     = A_->getGlobalNumRows();
00342   NumGlobalNonzeros_ = A_->getGlobalNumEntries();
00343 
00344   // NTS: Will need to add support for Zoltan2 partitions later
00345   // Also, will need a better way of handling the Graph typing issue.  Especially with ordinal types w.r.t the container
00346 
00347   if (PartitionerType_ == "linear")
00348     Partitioner_ = Teuchos::rcp( new Ifpack2::LinearPartitioner<Tpetra::RowGraph<LocalOrdinal,GlobalOrdinal,Node> >(A_->getGraph()) );
00349   else
00350     TEUCHOS_TEST_FOR_EXCEPTION(0, std::runtime_error,"Ifpack2::BlockRelaxation::initialize, invalid partitioner type.");
00351 
00352   // need to partition the graph of A
00353   Partitioner_->setParameters(List_);
00354   Partitioner_->compute();
00355 
00356   // get actual number of partitions
00357   NumLocalBlocks_ = Partitioner_->numLocalParts();
00358   
00359   // Note: Unlike Ifpack, we'll punt on computing W_ until compute(), which is where
00360   // we assume that the type of relaxation has been chosen.
00361   
00362   if (A_->getComm()->getSize() != 1)
00363     IsParallel_ = true;
00364   else
00365     IsParallel_ = false;
00366 
00367   ++NumInitialize_;
00368   Time_->stop();
00369   InitializeTime_ += Time_->totalElapsedTime();
00370   IsInitialized_ = true;
00371 }
00372 
00373 //==========================================================================
00374 template<class MatrixType,class ContainerType>
00375 void BlockRelaxation<MatrixType,ContainerType>::compute()
00376 {
00377   if (!isInitialized()) {
00378     initialize();
00379   }
00380 
00381   Time_->start(true);
00382 
00383   // reset values
00384   IsComputed_ = false;
00385   Condest_ = -1.0;
00386 
00387   TEUCHOS_TEST_FOR_EXCEPTION(NumSweeps_ < 0, std::runtime_error,
00388     "Ifpack2::BlockRelaxation::compute, NumSweeps_ must be >= 0");
00389 
00390   // Extract the submatrices
00391   ExtractSubmatrices();
00392 
00393   // Compute the weight vector if we're doing overlapped Jacobi (and only if we're doing overlapped Jacobi).
00394   if (PrecType_ == Ifpack2::JACOBI && OverlapLevel_ > 0) {
00395     // weight of each vertex
00396     W_ = Teuchos::rcp( new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A_->getRowMap()) );
00397     W_->putScalar(Teuchos::ScalarTraits<Scalar>::zero());
00398     Teuchos::ArrayRCP<Scalar > w_ptr = W_->getDataNonConst(0);
00399 
00400     for (LocalOrdinal i = 0 ; i < NumLocalBlocks_ ; ++i) {    
00401       for (size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
00402   int LID = (*Partitioner_)(i,j);
00403   w_ptr[LID]+= Teuchos::ScalarTraits<Scalar>::one();
00404       }
00405     }
00406     W_->reciprocal(*W_);
00407   }
00408 
00409   // We need to import data from external processors. Here I create a
00410   // Tpetra::Import object if needed (stealing from A_ if possible) 
00411   // Marzio's comment:
00412   // Note that I am doing some strange stuff to set the components of Y
00413   // from Y2 (to save some time).
00414   //
00415   if (IsParallel_ && ((PrecType_ == Ifpack2::GS) || (PrecType_ == Ifpack2::SGS))) {
00416     Importer_=A_->getGraph()->getImporter();
00417     if(Importer_==Teuchos::null)
00418       Importer_ = Teuchos::rcp( new Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>(A_->getDomainMap(),
00419                         A_->getColMap()) );
00420 
00421     TEUCHOS_TEST_FOR_EXCEPTION(Importer_ == Teuchos::null, std::runtime_error,
00422       "Ifpack2::BlockRelaxation::compute ERROR failed to create Importer_");
00423   }
00424 
00425   ++NumCompute_;
00426   Time_->stop();
00427   ComputeTime_ += Time_->totalElapsedTime();
00428   IsComputed_ = true;
00429 }
00430 
00431 //==============================================================================
00432 template<class MatrixType,class ContainerType>
00433 void BlockRelaxation<MatrixType,ContainerType>::ExtractSubmatrices()
00434 {
00435   TEUCHOS_TEST_FOR_EXCEPTION(Partitioner_==Teuchos::null, std::runtime_error,
00436            "Ifpack2::BlockRelaxation::ExtractSubmatrices, partitioner is null.");
00437 
00438   NumLocalBlocks_ = Partitioner_->numLocalParts();
00439 
00440   Containers_.resize(NumLocalBlocks_);
00441 
00442   for (LocalOrdinal i = 0 ; i < NumLocalBlocks_ ; ++i) {
00443     size_t rows = Partitioner_->numRowsInPart(i);
00444     Containers_[i] = Teuchos::rcp( new ContainerType(rows) );
00445     TEUCHOS_TEST_FOR_EXCEPTION(Containers_[i]==Teuchos::null, std::runtime_error,
00446            "Ifpack2::BlockRelaxation::ExtractSubmatrices, container consturctor failed.");
00447     
00448     Containers_[i]->setParameters(List_);
00449     Containers_[i]->initialize();
00450     // flops in initialize() will be computed on-the-fly in method initializeFlops().
00451 
00452     // set "global" ID of each partitioner row
00453     for (size_t j = 0 ; j < rows ; ++j) {
00454       Containers_[i]->ID(j) = (*Partitioner_)(i,j);
00455     }
00456 
00457     Containers_[i]->compute(A_);
00458     // flops in compute() will be computed on-the-fly in method computeFlops().
00459   }
00460 }
00461 
00462 
00463 
00464 //==========================================================================
00465 template<class MatrixType,class ContainerType>
00466 void BlockRelaxation<MatrixType,ContainerType>::ApplyInverseJacobi(
00467         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00468               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00469 {
00470   size_t NumVectors = X.getNumVectors();
00471   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> AY( Y.getMap(),NumVectors );
00472   
00473   // Initial matvec not needed
00474   int starting_iteration=0;
00475   if(ZeroStartingSolution_) {
00476     DoJacobi(X,Y);
00477     starting_iteration=1;
00478   }
00479 
00480   for (int j = starting_iteration; j < NumSweeps_ ; j++) {       
00481     applyMat(Y,AY);
00482     AY.update(1.0,X,-1.0);
00483     DoJacobi(AY,Y);
00484 
00485     // Flops for matrix apply & update
00486     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
00487   }
00488 
00489 }
00490 
00491 
00492 //==========================================================================
00493 template<class MatrixType,class ContainerType>
00494 void BlockRelaxation<MatrixType,ContainerType>::DoJacobi(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00495 {
00496   size_t NumVectors = X.getNumVectors();
00497   Scalar one=Teuchos::ScalarTraits<Scalar>::one();
00498   // Note: Flop counts copied naively from Ifpack.
00499 
00500   if (OverlapLevel_ == 0) {
00501     // Non-overlapping Jacobi
00502     for (LocalOrdinal i = 0 ; i < NumLocalBlocks_ ; i++) {     
00503       // may happen that a partition is empty
00504       if (Containers_[i]->getNumRows() == 0) continue;
00505       Containers_[i]->apply(X,Y,Teuchos::NO_TRANS,DampingFactor_,one);     
00506       ApplyFlops_ += NumVectors * 2 * NumGlobalRows_;
00507     }
00508   }
00509   else {
00510     // Overlapping Jacobi
00511     for (LocalOrdinal i = 0 ; i < NumLocalBlocks_ ; i++) {
00512       // may happen that a partition is empty
00513       if (Containers_[i]->getNumRows() == 0) continue;
00514       Containers_[i]->weightedApply(X,Y,*W_,Teuchos::NO_TRANS,DampingFactor_,one);
00515       // NOTE: do not count (for simplicity) the flops due to overlapping rows
00516       ApplyFlops_ += NumVectors * 4 * NumGlobalRows_;
00517     }
00518   }
00519 }
00520 
00521 //==========================================================================
00522 template<class MatrixType,class ContainerType>
00523 void BlockRelaxation<MatrixType,ContainerType>::ApplyInverseGS(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00524                      Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00525 {
00526 
00527   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>  Xcopy(X);
00528   for (int j = 0; j < NumSweeps_ ; j++) {
00529     DoGaussSeidel(Xcopy,Y);
00530     if(j != NumSweeps_-1)
00531       Xcopy=X;
00532   }
00533 }
00534 
00535 
00536 //==============================================================================
00537 template<class MatrixType,class ContainerType>
00538 void BlockRelaxation<MatrixType,ContainerType>::DoGaussSeidel(Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00539                     Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00540 {
00541   // Note: Flop counts copied naively from Ifpack.
00542   
00543   // allocations
00544   Scalar one=Teuchos::ScalarTraits<Scalar>::one();
00545   int Length = A_->getNodeMaxNumRowEntries(); 
00546   int NumVectors = X.getNumVectors();
00547   Teuchos::Array<Scalar>         Values;
00548   Teuchos::Array<LocalOrdinal>   Indices;
00549   Values.resize(Length);
00550   Indices.resize(Length);
00551 
00552   // an additonal vector is needed by parallel computations
00553   // (note that applications through Ifpack2_AdditiveSchwarz
00554   // are always seen are serial)
00555   Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node > > Y2;
00556   if (IsParallel_)
00557     Y2 = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Importer_->getTargetMap(), NumVectors) );
00558   else
00559     Y2 = Teuchos::rcp( &Y, false );
00560 
00561 
00562   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       x_ptr   = X.get2dViewNonConst();
00563   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       y_ptr   = Y.get2dViewNonConst();
00564   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       y2_ptr  = Y2->get2dViewNonConst();
00565 
00566   // data exchange is here, once per sweep
00567   if (IsParallel_)  Y2->doImport(Y,*Importer_,Tpetra::INSERT);
00568 
00569   for (LocalOrdinal i = 0 ; i < NumLocalBlocks_ ; i++) {
00570     // may happen that a partition is empty
00571     if (Containers_[i]->getNumRows() == 0) continue;
00572     LocalOrdinal LID;
00573 
00574     // update from previous block
00575     for (size_t j = 0 ; j < Containers_[i]->getNumRows(); j++) {
00576       LID = Containers_[i]->ID(j);
00577       size_t NumEntries;
00578       A_->getLocalRowCopy(LID,Indices(),Values(),NumEntries);
00579 
00580       for (size_t k = 0 ; k < NumEntries ; k++) {
00581         LocalOrdinal col = Indices[k];
00582   for (int kk = 0 ; kk < NumVectors ; kk++)
00583     x_ptr[kk][LID] -= Values[k] * y2_ptr[kk][col];  
00584       }
00585     }
00586     // solve with this block
00587     // Note: I'm abusing the ordering information, knowing that X/Y and Y2 have the same ordering for on-proc unknowns.
00588     // Note: Add flop counts for inverse apply
00589     Containers_[i]->apply(X,*Y2,Teuchos::NO_TRANS,DampingFactor_,one);
00590     
00591     // operations for all getrow's
00592     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);    
00593   } // end for NumLocalBlocks_
00594 
00595   // Attention: this is delicate... Not all combinations
00596   // of Y2 and Y will always work (tough for ML it should be ok)
00597   if (IsParallel_)
00598     for (int m = 0 ; m < NumVectors ; ++m) 
00599       for (size_t i = 0 ; i < NumMyRows_ ; ++i)
00600         y_ptr[m][i] = y2_ptr[m][i];
00601 
00602 }
00603 
00604 //==========================================================================
00605 template<class MatrixType,class ContainerType>
00606 void BlockRelaxation<MatrixType,ContainerType>::ApplyInverseSGS(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00607                 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00608 {
00609   
00610   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>  Xcopy(X);
00611   for (int j = 0; j < NumSweeps_ ; j++) {
00612     DoSGS(X,Xcopy,Y);
00613     if(j != NumSweeps_-1)
00614       Xcopy=X;
00615   }
00616 }
00617 
00618 //==========================================================================
00619 template<class MatrixType,class ContainerType>
00620 void BlockRelaxation<MatrixType,ContainerType>::DoSGS(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00621                   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Xcopy, 
00622                   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00623 {
00624   Scalar one=Teuchos::ScalarTraits<Scalar>::one();
00625   int Length = A_->getNodeMaxNumRowEntries(); 
00626   int NumVectors = X.getNumVectors();
00627   Teuchos::Array<Scalar>         Values;
00628   Teuchos::Array<LocalOrdinal>   Indices;
00629   Values.resize(Length);
00630   Indices.resize(Length);
00631 
00632   // an additonal vector is needed by parallel computations
00633   // (note that applications through Ifpack2_AdditiveSchwarz
00634   // are always seen are serial)
00635   Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node > > Y2;
00636   if (IsParallel_)
00637     Y2 = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Importer_->getTargetMap(), NumVectors) );
00638   else
00639     Y2 = Teuchos::rcp( &Y, false );
00640 
00641   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr       = X.get2dView();
00642   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       xcopy_ptr   = Xcopy.get2dViewNonConst();
00643   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       y_ptr       = Y.get2dViewNonConst();
00644   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       y2_ptr      = Y2->get2dViewNonConst();
00645 
00646   // data exchange is here, once per sweep
00647   if (IsParallel_)  Y2->doImport(Y,*Importer_,Tpetra::INSERT);
00648 
00649   // Forward Sweep
00650   for(LocalOrdinal i = 0 ; i < NumLocalBlocks_ ; i++) {
00651     // may happen that a partition is empty
00652     if (Containers_[i]->getNumRows() == 0) continue;
00653     LocalOrdinal LID;
00654 
00655     // update from previous block
00656     for(size_t j = 0 ; j < Containers_[i]->getNumRows(); j++) {
00657       LID = Containers_[i]->ID(j);
00658       size_t NumEntries;
00659       A_->getLocalRowCopy(LID,Indices(),Values(),NumEntries);
00660 
00661       for (size_t k = 0 ; k < NumEntries ; k++) {
00662         LocalOrdinal col = Indices[k];
00663         for (int kk = 0 ; kk < NumVectors ; kk++) 
00664           xcopy_ptr[kk][LID] -= Values[k] * y2_ptr[kk][col];
00665       }
00666     }
00667     // solve with this block
00668     // Note: I'm abusing the ordering information, knowing that X/Y and Y2 have the same ordering for on-proc unknowns.
00669     // Note: Add flop counts for inverse apply
00670     Containers_[i]->apply(Xcopy,*Y2,Teuchos::NO_TRANS,DampingFactor_,one);
00671 
00672     // operations for all getrow's
00673     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
00674   }// end forward sweep
00675 
00676   // Reverse Sweep
00677   Xcopy = X;
00678   for(LocalOrdinal i = NumLocalBlocks_-1; i >=0 ; i--) {
00679     // may happen that a partition is empty
00680     if (Containers_[i]->getNumRows() == 0) continue;
00681     LocalOrdinal LID;
00682 
00683     // update from previous block
00684     for(size_t j = 0 ; j < Containers_[i]->getNumRows(); j++) {
00685       LID = Containers_[i]->ID(j);
00686       size_t NumEntries;
00687       A_->getLocalRowCopy(LID,Indices(),Values(),NumEntries);
00688 
00689       for (size_t k = 0 ; k < NumEntries ; k++) {
00690         LocalOrdinal col = Indices[k];
00691         for (int kk = 0 ; kk < NumVectors ; kk++) 
00692           xcopy_ptr[kk][LID] -= Values[k] * y2_ptr[kk][col];
00693       }
00694     }
00695     
00696     // solve with this block
00697     // Note: I'm abusing the ordering information, knowing that X/Y and Y2 have the same ordering for on-proc unknowns.
00698     // Note: Add flop counts for inverse apply
00699     Containers_[i]->apply(Xcopy,*Y2,Teuchos::NO_TRANS,DampingFactor_,one);
00700 
00701     // operations for all getrow's
00702     ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
00703   } //end reverse sweep
00704 
00705 
00706   // Attention: this is delicate... Not all combinations
00707   // of Y2 and Y will always work (tough for ML it should be ok)
00708   if (IsParallel_)
00709     for (int m = 0 ; m < NumVectors ; ++m) 
00710       for (size_t i = 0 ; i < NumMyRows_ ; ++i)
00711         y_ptr[m][i] = y2_ptr[m][i];
00712 }
00713 
00714 //==========================================================================
00715 template<class MatrixType,class ContainerType>
00716 std::string BlockRelaxation<MatrixType,ContainerType>::description() const {
00717   std::ostringstream oss;
00718   oss << Teuchos::Describable::description();
00719   if (isInitialized()) {
00720     if (isComputed()) {
00721       oss << "{status = initialized, computed";
00722     }
00723     else {
00724       oss << "{status = initialized, not computed";
00725     }
00726   }
00727   else {
00728     oss << "{status = not initialized, not computed";
00729   }
00730   //
00731   if (PrecType_ == Ifpack2::JACOBI)   oss << "Type = Block Jacobi, " << std::endl;
00732   else if (PrecType_ == Ifpack2::GS)  oss << "Type = Block Gauss-Seidel, " << std::endl;
00733   else if (PrecType_ == Ifpack2::SGS) oss << "Type = Block Sym. Gauss-Seidel, " << std::endl;
00734   //
00735   oss << ", global rows = " << A_->getGlobalNumRows()
00736       << ", global cols = " << A_->getGlobalNumCols();
00737 
00738   oss << "}";
00739   return oss.str();
00740 }
00741 
00742 //==========================================================================
00743 template<class MatrixType,class ContainerType>
00744 void BlockRelaxation<MatrixType,ContainerType>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00745   using std::endl;
00746   using std::setw;
00747   using Teuchos::VERB_DEFAULT;
00748   using Teuchos::VERB_NONE;
00749   using Teuchos::VERB_LOW;
00750   using Teuchos::VERB_MEDIUM;
00751   using Teuchos::VERB_HIGH;
00752   using Teuchos::VERB_EXTREME;
00753   Teuchos::EVerbosityLevel vl = verbLevel;
00754   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00755   const int myImageID = A_->getComm()->getRank();
00756   Teuchos::OSTab tab(out);
00757 
00758   //    none: print nothing
00759   //     low: print O(1) info from node 0
00760   //  medium: 
00761   //    high: 
00762   // extreme: 
00763   if (vl != VERB_NONE && myImageID == 0) {
00764     out << this->description() << endl;
00765     out << endl;
00766     out << "===============================================================================" << endl;
00767     out << "Sweeps         = " << NumSweeps_ << endl;
00768     out << "damping factor = " << DampingFactor_ << endl;
00769     if (PrecType_ == Ifpack2::GS && DoBackwardGS_) {
00770       out << "Using backward mode (BGS only)" << endl;
00771     }
00772     if   (ZeroStartingSolution_) { out << "Using zero starting solution" << endl; }
00773     else                         { out << "Using input starting solution" << endl; }
00774     if   (Condest_ == -1.0) { out << "Condition number estimate       = N/A" << endl; }
00775     else                    { out << "Condition number estimate       = " << Condest_ << endl; }
00776     out << endl;
00777     out << "Phase           # calls    Total Time (s)     Total MFlops      MFlops/s       " << endl;
00778     out << "------------    -------    ---------------    ---------------   ---------------" << endl;
00779     out << setw(12) << "initialize()" << setw(5) << getNumInitialize() << "    " << setw(15) << getInitializeTime() << endl;
00780     out << setw(12) << "compute()" << setw(5) << getNumCompute()    << "    " << setw(15) << getComputeTime() << "    " 
00781         << setw(15) << getComputeFlops() << "    " 
00782         << setw(15) << (getComputeTime() != 0.0 ? getComputeFlops() / getComputeTime() * 1.0e-6 : 0.0) << endl;
00783     out << setw(12) << "apply()" << setw(5) << getNumApply()    << "    " << setw(15) << getApplyTime() << "    " 
00784         << setw(15) << getApplyFlops() << "    " 
00785         << setw(15) << (getApplyTime() != 0.0 ? getApplyFlops() / getApplyTime() * 1.0e-6 : 0.0) << endl;
00786     out << "===============================================================================" << endl;
00787     out << endl;
00788   }
00789 }
00790 
00791 }//namespace Ifpack2
00792 
00793 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP
00794 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends