Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Relaxation_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_RELAXATION_DEF_HPP
00031 #define IFPACK2_RELAXATION_DEF_HPP
00032 
00033 #include "Ifpack2_Relaxation_decl.hpp"
00034 
00035 namespace Ifpack2 {
00036 
00037 //==========================================================================
00038 template<class MatrixType>
00039 Relaxation<MatrixType>::Relaxation(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A)
00040 : A_(A),
00041   Comm_(A->getRowMap()->getComm()),
00042   Time_( Teuchos::rcp( new Teuchos::Time("Ifpack2::Relaxation") ) ),
00043   NumSweeps_(1),
00044   PrecType_(Ifpack2::JACOBI),
00045   MinDiagonalValue_(0.0),
00046   DampingFactor_(1.0),
00047   IsParallel_(false),
00048   ZeroStartingSolution_(true),
00049   DoBackwardGS_(false),
00050   DoL1Method_(false),
00051   L1Eta_(1.5),
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) << "::Relaxation(): input matrix reference was null.");
00069 }
00070 
00071 //==========================================================================
00072 template<class MatrixType>
00073 Relaxation<MatrixType>::~Relaxation() {
00074 }
00075 
00076 //==========================================================================
00077 template<class MatrixType>
00078 void Relaxation<MatrixType>::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::Relaxation::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, "relaxation: use l1",DoL1Method_);
00113   Ifpack2::getParameter(List, "relaxation: l1 eta",L1Eta_);
00114 }
00115 
00116 //==========================================================================
00117 template<class MatrixType>
00118 const Teuchos::RCP<const Teuchos::Comm<int> > & 
00119 Relaxation<MatrixType>::getComm() const{
00120   return(Comm_);
00121 }
00122 
00123 //==========================================================================
00124 template<class MatrixType>
00125 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
00126                                      typename MatrixType::local_ordinal_type,
00127                                      typename MatrixType::global_ordinal_type,
00128                                      typename MatrixType::node_type> >
00129 Relaxation<MatrixType>::getMatrix() const {
00130   return(A_);
00131 }
00132 
00133 //==========================================================================
00134 template<class MatrixType>
00135 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00136                                      typename MatrixType::global_ordinal_type,
00137                                      typename MatrixType::node_type> >&
00138 Relaxation<MatrixType>::getDomainMap() const {
00139   return A_->getDomainMap();
00140 }
00141 
00142 //==========================================================================
00143 template<class MatrixType>
00144 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
00145                                      typename MatrixType::global_ordinal_type,
00146                                      typename MatrixType::node_type> >&
00147 Relaxation<MatrixType>::getRangeMap() const {
00148   return A_->getRangeMap();
00149 }
00150 
00151 //==========================================================================
00152 template<class MatrixType>
00153 bool Relaxation<MatrixType>::hasTransposeApply() const {
00154   return true;
00155 }
00156 
00157 //==========================================================================
00158 template<class MatrixType>
00159 int Relaxation<MatrixType>::getNumInitialize() const {
00160   return(NumInitialize_);
00161 }
00162 
00163 //==========================================================================
00164 template<class MatrixType>
00165 int Relaxation<MatrixType>::getNumCompute() const {
00166   return(NumCompute_);
00167 }
00168 
00169 //==========================================================================
00170 template<class MatrixType>
00171 int Relaxation<MatrixType>::getNumApply() const {
00172   return(NumApply_);
00173 }
00174 
00175 //==========================================================================
00176 template<class MatrixType>
00177 double Relaxation<MatrixType>::getInitializeTime() const {
00178   return(InitializeTime_);
00179 }
00180 
00181 //==========================================================================
00182 template<class MatrixType>
00183 double Relaxation<MatrixType>::getComputeTime() const {
00184   return(ComputeTime_);
00185 }
00186 
00187 //==========================================================================
00188 template<class MatrixType>
00189 double Relaxation<MatrixType>::getApplyTime() const {
00190   return(ApplyTime_);
00191 }
00192 
00193 //==========================================================================
00194 template<class MatrixType>
00195 double Relaxation<MatrixType>::getComputeFlops() const {
00196   return(ComputeFlops_);
00197 }
00198 
00199 //==========================================================================
00200 template<class MatrixType>
00201 double Relaxation<MatrixType>::getApplyFlops() const {
00202   return(ApplyFlops_);
00203 }
00204 
00205 //==========================================================================
00206 template<class MatrixType>
00207 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00208 Relaxation<MatrixType>::getCondEst() const
00209 {
00210   return(Condest_);
00211 }
00212 
00213 //==========================================================================
00214 template<class MatrixType>
00215 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00216 Relaxation<MatrixType>::computeCondEst(
00217                      CondestType CT,
00218                      typename MatrixType::local_ordinal_type MaxIters, 
00219                      magnitudeType Tol,
00220      const Teuchos::Ptr<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
00221                                                 typename MatrixType::local_ordinal_type,
00222                                                 typename MatrixType::global_ordinal_type,
00223                                                 typename MatrixType::node_type> > &matrix)
00224 {
00225   if (!isComputed()) // cannot compute right now
00226     return(-1.0);
00227 
00228   // always compute it. Call Condest() with no parameters to get
00229   // the previous estimate.
00230   Condest_ = Ifpack2::Condest(*this, CT, MaxIters, Tol, matrix);
00231 
00232   return(Condest_);
00233 }
00234 
00235 //==========================================================================
00236 template<class MatrixType>
00237 void Relaxation<MatrixType>::apply(
00238           const Tpetra::MultiVector<typename MatrixType::scalar_type,
00239                                     typename MatrixType::local_ordinal_type,
00240                                     typename MatrixType::global_ordinal_type,
00241                                     typename MatrixType::node_type>& X,
00242                 Tpetra::MultiVector<typename MatrixType::scalar_type,
00243                                     typename MatrixType::local_ordinal_type,
00244                                     typename MatrixType::global_ordinal_type,
00245                                     typename MatrixType::node_type>& Y,
00246                 Teuchos::ETransp mode,
00247                  Scalar alpha,
00248                  Scalar beta) const
00249 {
00250   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == false, std::runtime_error,
00251      "Ifpack2::Relaxation::apply ERROR: isComputed() must be true prior to calling apply.");
00252 
00253   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00254      "Ifpack2::Relaxation::apply ERROR: X.getNumVectors() != Y.getNumVectors().");
00255 
00256   Time_->start(true);
00257 
00258   // If X and Y are pointing to the same memory location,
00259   // we need to create an auxiliary vector, Xcopy
00260   Teuchos::RCP< const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xcopy;
00261   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues())
00262     Xcopy = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X) );
00263   else
00264     Xcopy = Teuchos::rcp( &X, false );
00265 
00266   if (ZeroStartingSolution_)
00267     Y.putScalar(0.0);
00268 
00269   // Flops are updated in each of the following.
00270   switch (PrecType_) {
00271   case Ifpack2::JACOBI:
00272     ApplyInverseJacobi(*Xcopy,Y);
00273     break;
00274   case Ifpack2::GS:
00275     ApplyInverseGS(*Xcopy,Y);
00276     break;
00277   case Ifpack2::SGS:
00278     ApplyInverseSGS(*Xcopy,Y);
00279     break;
00280   default:
00281     throw std::runtime_error("Ifpack2::Relaxation::apply internal logic error.");
00282   }
00283 
00284   ++NumApply_;
00285   Time_->stop();
00286   ApplyTime_ += Time_->totalElapsedTime();
00287 }
00288 
00289 //==========================================================================
00290 template<class MatrixType>
00291 void Relaxation<MatrixType>::applyMat(
00292           const Tpetra::MultiVector<typename MatrixType::scalar_type,
00293                                     typename MatrixType::local_ordinal_type,
00294                                     typename MatrixType::global_ordinal_type,
00295                                     typename MatrixType::node_type>& X,
00296                 Tpetra::MultiVector<typename MatrixType::scalar_type,
00297                                     typename MatrixType::local_ordinal_type,
00298                                     typename MatrixType::global_ordinal_type,
00299                                     typename MatrixType::node_type>& Y,
00300              Teuchos::ETransp mode) const
00301 {
00302   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == false, std::runtime_error,
00303      "Ifpack2::Relaxation::applyMat() ERROR: isComputed() must be true prior to calling applyMat().");
00304   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00305      "Ifpack2::Relaxation::applyMat() ERROR: X.getNumVectors() != Y.getNumVectors().");
00306   A_->apply(X, Y, mode);
00307 }
00308 
00309 //==========================================================================
00310 template<class MatrixType>
00311 void Relaxation<MatrixType>::initialize() {
00312   IsInitialized_ = false;
00313 
00314   TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error,
00315     "Ifpack2::Relaxation::Initialize ERROR, Matrix is NULL");
00316 
00317   Time_->start(true);
00318   
00319   NumMyRows_ = A_->getNodeNumRows();
00320   NumGlobalRows_ = A_->getGlobalNumRows();
00321   NumGlobalNonzeros_ = A_->getGlobalNumEntries();
00322 
00323   if (Comm_->getSize() != 1)
00324     IsParallel_ = true;
00325   else
00326     IsParallel_ = false;
00327 
00328   ++NumInitialize_;
00329   Time_->stop();
00330   InitializeTime_ += Time_->totalElapsedTime();
00331   IsInitialized_ = true;
00332 }
00333 
00334 //==========================================================================
00335 template<class MatrixType>
00336 void Relaxation<MatrixType>::compute()
00337 {
00338   if (!isInitialized()) {
00339     initialize();
00340   }
00341 
00342   Time_->start(true);
00343 
00344   // reset values
00345   IsComputed_ = false;
00346   Condest_ = -1.0;
00347 
00348   TEUCHOS_TEST_FOR_EXCEPTION(NumSweeps_ < 0, std::runtime_error,
00349     "Ifpack2::Relaxation::compute, NumSweeps_ must be >= 0");
00350   
00351   Diagonal_ = Teuchos::rcp( new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A_->getRowMap()) );
00352 
00353   TEUCHOS_TEST_FOR_EXCEPTION(Diagonal_ == Teuchos::null, std::runtime_error,
00354     "Ifpack2::Relaxation::compute, failed to create Diagonal_");
00355 
00356   A_->getLocalDiagCopy(*Diagonal_);
00357   Teuchos::ArrayRCP<Scalar> DiagView = Diagonal_->get1dViewNonConst();
00358   
00359   // Setup for L1 Methods.
00360   // Here we add half the value of the off-processor entries in the row, 
00361   // but only if diagonal isn't sufficiently large.
00362   //
00363   // Note: This is only done in the slower-but-more-general "RowMatrix" mode.  
00364   //
00365   // This follows from Equation (6.5) in:
00366   // Baker, Falgout, Kolev and Yang.  Multigrid Smoothers for Ultraparallel Computing.
00367   // SIAM J. Sci. Comput., Vol. 33, No. 5. (2011), pp. 2864--2887.
00368   if(DoL1Method_ && IsParallel_) {
00369     size_t maxLength = A_->getNodeMaxNumRowEntries();
00370     Teuchos::Array<LocalOrdinal> Indices(maxLength);
00371     Teuchos::Array<Scalar> Values(maxLength);
00372     size_t NumEntries;
00373 
00374     Scalar two=Teuchos::ScalarTraits<Scalar>::one()+Teuchos::ScalarTraits<Scalar>::one();
00375 
00376     for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00377       A_->getLocalRowCopy(i, Indices(), Values(), NumEntries);      
00378       magnitudeType diagonal_boost=Teuchos::ScalarTraits<magnitudeType>::zero();
00379       for (size_t k = 0 ; k < NumEntries ; ++k) 
00380   if((size_t)Indices[k] > i)
00381     diagonal_boost+= Teuchos::ScalarTraits<Scalar>::magnitude(Values[k]/two);  
00382       if (Teuchos::ScalarTraits<Scalar>::magnitude(DiagView[i]) < L1Eta_*diagonal_boost)
00383   DiagView[i]+=diagonal_boost;
00384     }
00385   }
00386 
00387   // check diagonal elements, store the inverses, and verify that
00388   // no zeros are around. If an element is zero, then by default
00389   // its inverse is zero as well (that is, the row is ignored).
00390   for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00391     Scalar& diag = DiagView[i];
00392     if (Teuchos::ScalarTraits<Scalar>::magnitude(diag) < Teuchos::ScalarTraits<Scalar>::magnitude(MinDiagonalValue_))
00393       diag = MinDiagonalValue_;
00394     if (diag != Teuchos::ScalarTraits<Scalar>::zero())
00395       diag = Teuchos::ScalarTraits<Scalar>::one() / diag;
00396   }
00397   ComputeFlops_ += NumMyRows_;
00398 
00399 
00400   // We need to import data from external processors. Here I create a
00401   // Tpetra::Import object if needed (stealing from A_ if possible) 
00402   // Marzio's comment:
00403   // Note that I am doing some strange stuff to set the components of Y
00404   // from Y2 (to save some time).
00405   //
00406   if (IsParallel_ && ((PrecType_ == Ifpack2::GS) || (PrecType_ == Ifpack2::SGS))) {
00407     Importer_=A_->getGraph()->getImporter();
00408     if(Importer_==Teuchos::null)
00409       Importer_ = Teuchos::rcp( new Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>(A_->getDomainMap(),
00410                         A_->getColMap()) );
00411 
00412     TEUCHOS_TEST_FOR_EXCEPTION(Importer_ == Teuchos::null, std::runtime_error,
00413       "Ifpack2::Relaxation::compute ERROR failed to create Importer_");
00414   }
00415 
00416 
00417   ++NumCompute_;
00418   Time_->stop();
00419   ComputeTime_ += Time_->totalElapsedTime();
00420   IsComputed_ = true;
00421 }
00422 
00423 //==========================================================================
00424 template<class MatrixType>
00425 void Relaxation<MatrixType>::ApplyInverseJacobi(
00426         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00427               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00428 {
00429   int NumVectors = X.getNumVectors();
00430   Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> A_times_Y( Y.getMap(),NumVectors );
00431 
00432   for (int j = 0; j < NumSweeps_ ; j++) {
00433     applyMat(Y,A_times_Y);
00434     A_times_Y.update(1.0,X,-1.0);
00435     Y.elementWiseMultiply(DampingFactor_, *Diagonal_, A_times_Y, 1.0);
00436   }
00437 
00438   // Flops:
00439   // - matrix vector              (2 * NumGlobalNonzeros_)
00440   // - update                     (2 * NumGlobalRows_)
00441   // - Multiply:
00442   //   - DampingFactor            (NumGlobalRows_)
00443   //   - Diagonal                 (NumGlobalRows_)
00444   //   - A + B                    (NumGlobalRows_)
00445   //   - 1.0                      (NumGlobalRows_)
00446   ApplyFlops_ += NumVectors * (6 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
00447 }
00448 
00449 //==========================================================================
00450 template<class MatrixType>
00451 void Relaxation<MatrixType>::ApplyInverseGS(
00452         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00453               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00454 {
00455   const MatrixType* CrsMatrix = dynamic_cast<const MatrixType*>(&*A_);
00456   // try to pick the best option; performances may be improved
00457   // if several sweeps are used.
00458   if (CrsMatrix != 0)
00459   {
00460     ApplyInverseGS_CrsMatrix(*CrsMatrix, X, Y);
00461   }
00462   else {
00463     ApplyInverseGS_RowMatrix(X, Y);
00464   }
00465 }
00466 
00467 //==========================================================================
00468 template<class MatrixType>
00469 void Relaxation<MatrixType>::ApplyInverseGS_RowMatrix(
00470         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00471               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00472 {
00473   size_t NumVectors = X.getNumVectors();
00474 
00475   size_t maxLength = A_->getNodeMaxNumRowEntries();
00476   Teuchos::Array<LocalOrdinal> Indices(maxLength);
00477   Teuchos::Array<Scalar> Values(maxLength);
00478 
00479   Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Y2;
00480   if (IsParallel_)
00481     Y2 = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Importer_->getTargetMap(), NumVectors) );
00482   else
00483     Y2 = Teuchos::rcp( &Y, false );
00484 
00485   // extract views (for nicer and faster code)
00486   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y_ptr = Y.get2dViewNonConst();
00487   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y2_ptr = Y2->get2dViewNonConst();
00488   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr =  X.get2dView();
00489   Teuchos::ArrayRCP<const Scalar> d_ptr = Diagonal_->get1dView();
00490  
00491   for (int j = 0; j < NumSweeps_ ; j++) {
00492 
00493     // data exchange is here, once per sweep
00494     if (IsParallel_)
00495       Y2->doImport(Y,*Importer_,Tpetra::INSERT);
00496 
00497     if(!DoBackwardGS_){      
00498       /* Forward Mode */
00499       for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00500         
00501         size_t NumEntries;
00502         A_->getLocalRowCopy(i, Indices(), Values(), NumEntries);
00503         
00504         for (size_t m = 0 ; m < NumVectors ; ++m) {
00505 
00506           Scalar dtemp = 0.0;
00507           for (size_t k = 0 ; k < NumEntries ; ++k) {
00508             LocalOrdinal col = Indices[k];
00509             dtemp += Values[k] * y2_ptr[m][col];
00510           }
00511           
00512           y2_ptr[m][i] += DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
00513         }
00514       }
00515     }
00516     else {
00517       /* Backward Mode */
00518       for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {
00519         size_t NumEntries;
00520         A_->getLocalRowCopy(i, Indices(), Values(), NumEntries);
00521 
00522         for (size_t m = 0 ; m < NumVectors ; ++m) {
00523           
00524           Scalar dtemp = 0.0;
00525           for (size_t k = 0 ; k < NumEntries ; ++k) {
00526             LocalOrdinal col = Indices[k];
00527             dtemp += Values[k] * y2_ptr[m][col];
00528           }
00529 
00530           y2_ptr[m][i] += DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
00531         }
00532       }
00533     }
00534 
00535     // using Export() sounded quite expensive   
00536     if (IsParallel_) {
00537       for (size_t m = 0 ; m < NumVectors ; ++m) 
00538         for (size_t i = 0 ; i < NumMyRows_ ; ++i)
00539           y_ptr[m][i] = y2_ptr[m][i];
00540     }
00541   }
00542 
00543   ApplyFlops_ += NumVectors * (4 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
00544 }
00545 
00546 //==========================================================================
00547 template<class MatrixType>
00548 void Relaxation<MatrixType>::ApplyInverseGS_CrsMatrix(
00549         const MatrixType& A,
00550         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
00551               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00552 {
00553   size_t NumVectors = X.getNumVectors();
00554 
00555   Teuchos::ArrayView<const LocalOrdinal> Indices;
00556   Teuchos::ArrayView<const Scalar> Values;
00557 
00558   Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Y2;
00559   if (IsParallel_) {
00560     Y2 = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Importer_->getTargetMap(), NumVectors) );
00561   }
00562   else
00563     Y2 = Teuchos::rcp( &Y, false );
00564 
00565   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y_ptr = Y.get2dViewNonConst();
00566   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y2_ptr = Y2->get2dViewNonConst();
00567   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr =  X.get2dView();
00568   Teuchos::ArrayRCP<const Scalar> d_ptr = Diagonal_->get1dView();
00569   
00570   for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
00571     
00572     // only one data exchange per sweep
00573     if (IsParallel_)
00574       Y2->doImport(Y,*Importer_,Tpetra::INSERT);
00575 
00576     if(!DoBackwardGS_){  
00577       /* Forward Mode */
00578       for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00579 
00580         LocalOrdinal col;
00581         Scalar diag = d_ptr[i];
00582         
00583         A.getLocalRowView(i, Indices, Values);
00584         size_t NumEntries = Indices.size();
00585         
00586         for (size_t m = 0 ; m < NumVectors ; ++m) {
00587           
00588           Scalar dtemp = 0.0;
00589           
00590           for (size_t k = 0; k < NumEntries; ++k) {
00591             col = Indices[k];
00592             dtemp += Values[k] * y2_ptr[m][col];
00593           }
00594           
00595           y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00596         }
00597       }
00598     }
00599     else {
00600       /* Backward Mode */
00601       for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {
00602 
00603         LocalOrdinal col;
00604         Scalar diag = d_ptr[i];
00605         
00606         A.getLocalRowView(i, Indices, Values);
00607         size_t NumEntries = Indices.size();
00608         
00609         for (size_t m = 0 ; m < NumVectors ; ++m) {
00610           
00611           Scalar dtemp = 0.0;
00612           for (size_t k = 0; k < NumEntries; ++k) {
00613             col = Indices[k];
00614             dtemp += Values[k] * y2_ptr[m][col];
00615           }
00616           
00617           y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00618           
00619         }
00620       }
00621     }
00622     
00623     if (IsParallel_) {
00624       for (size_t m = 0 ; m < NumVectors ; ++m) 
00625         for (size_t i = 0 ; i < NumMyRows_ ; ++i)
00626           y_ptr[m][i] = y2_ptr[m][i];
00627     }
00628   }
00629 
00630   ApplyFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
00631 }
00632 
00633 //==========================================================================
00634 template<class MatrixType>
00635 void Relaxation<MatrixType>::ApplyInverseSGS(
00636         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00637               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00638 {
00639   const MatrixType* CrsMatrix = dynamic_cast<const MatrixType*>(&*A_);
00640   // try to pick the best option; performance may be improved
00641   // if several sweeps are used.
00642   if (CrsMatrix != 0)
00643   {
00644     ApplyInverseSGS_CrsMatrix(*CrsMatrix, X, Y);
00645   }
00646   else
00647     ApplyInverseSGS_RowMatrix(X, Y);
00648 }
00649 
00650 //==========================================================================
00651 template<class MatrixType>
00652 void Relaxation<MatrixType>::ApplyInverseSGS_RowMatrix(
00653         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00654               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00655 {
00656   size_t NumVectors = X.getNumVectors();
00657   size_t maxLength = A_->getNodeMaxNumRowEntries();
00658   Teuchos::Array<LocalOrdinal> Indices(maxLength);
00659   Teuchos::Array<Scalar> Values(maxLength);
00660 
00661   Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Y2;
00662   if (IsParallel_) {
00663     Y2 = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Importer_->getTargetMap(), NumVectors) );
00664   }
00665   else
00666     Y2 = Teuchos::rcp( &Y, false );
00667 
00668   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y_ptr = Y.get2dViewNonConst();
00669   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y2_ptr = Y2->get2dViewNonConst();
00670   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr =  X.get2dView();
00671   Teuchos::ArrayRCP<const Scalar> d_ptr = Diagonal_->get1dView();
00672   
00673   for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
00674     
00675     // only one data exchange per sweep
00676     if (IsParallel_)
00677       Y2->doImport(Y,*Importer_,Tpetra::INSERT);
00678 
00679     for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00680 
00681       size_t NumEntries;
00682       Scalar diag = d_ptr[i];
00683 
00684       A_->getLocalRowCopy(i, Indices(), Values(), NumEntries);
00685 
00686       for (size_t m = 0 ; m < NumVectors ; ++m) {
00687 
00688         Scalar dtemp = 0.0;
00689 
00690         for (size_t k = 0 ; k < NumEntries ; ++k) {
00691           LocalOrdinal col = Indices[k];
00692           dtemp += Values[k] * y2_ptr[m][col];
00693         }
00694 
00695         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00696       }
00697     }
00698 
00699     for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {
00700 
00701       size_t NumEntries;
00702       Scalar diag = d_ptr[i];
00703 
00704       A_->getLocalRowCopy(i, Indices(), Values(), NumEntries);
00705 
00706       for (size_t m = 0 ; m < NumVectors ; ++m) {
00707 
00708         Scalar dtemp = Teuchos::ScalarTraits<Scalar>::zero();
00709         for (size_t k = 0 ; k < NumEntries ; ++k) {
00710           LocalOrdinal col = Indices[k];
00711           dtemp += Values[k] * y2_ptr[m][col];
00712         }
00713 
00714         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00715       }
00716     }
00717 
00718     if (IsParallel_) {
00719       for (size_t m = 0 ; m < NumVectors ; ++m) 
00720         for (size_t i = 0 ; i < NumMyRows_ ; ++i)
00721           y_ptr[m][i] = y2_ptr[m][i];
00722     }
00723   }
00724 
00725   ApplyFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
00726 }
00727 
00728 //==========================================================================
00729 template<class MatrixType>
00730 void Relaxation<MatrixType>::ApplyInverseSGS_CrsMatrix(
00731         const MatrixType& A,
00732         const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00733               Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00734 {
00735   size_t NumVectors = X.getNumVectors();
00736 
00737   Teuchos::ArrayView<const LocalOrdinal> Indices;
00738   Teuchos::ArrayView<const Scalar> Values;
00739 
00740   Teuchos::RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Y2;
00741   if (IsParallel_) {
00742     Y2 = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Importer_->getTargetMap(), NumVectors) );
00743   }
00744   else
00745     Y2 = Teuchos::rcp( &Y, false );
00746 
00747   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y_ptr = Y.get2dViewNonConst();
00748   Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> > y2_ptr = Y2->get2dViewNonConst();
00749   Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr =  X.get2dView();
00750   Teuchos::ArrayRCP<const Scalar> d_ptr = Diagonal_->get1dView();
00751   
00752   for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
00753 
00754     // only one data exchange per sweep
00755     if (IsParallel_)
00756       Y2->doImport(Y,*Importer_,Tpetra::INSERT);
00757 
00758 
00759     for (size_t i = 0 ; i < NumMyRows_ ; ++i) {
00760 
00761       Scalar diag = d_ptr[i];
00762 
00763       A.getLocalRowView(i, Indices, Values);
00764       size_t NumEntries = Indices.size();
00765 
00766       for (size_t m = 0 ; m < NumVectors ; ++m) {
00767 
00768         Scalar dtemp = Teuchos::ScalarTraits<Scalar>::zero();
00769 
00770         for (size_t k = 0; k < NumEntries; ++k) {
00771 
00772           LocalOrdinal col = Indices[k];
00773           dtemp += Values[k] * y2_ptr[m][col];
00774         }
00775 
00776         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00777       }
00778     }
00779 
00780     for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {
00781 
00782       Scalar diag = d_ptr[i];
00783 
00784       A.getLocalRowView(i, Indices, Values);
00785       size_t NumEntries = Indices.size();
00786 
00787       for (size_t m = 0 ; m < NumVectors ; ++m) {
00788 
00789         Scalar dtemp = Teuchos::ScalarTraits<Scalar>::zero();
00790         for (size_t k = 0; k < NumEntries; ++k) {
00791 
00792           LocalOrdinal col = Indices[k];
00793           dtemp += Values[k] * y2_ptr[m][col];
00794         }
00795 
00796         y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
00797 
00798       }
00799     }
00800 
00801     if (IsParallel_) {
00802       for (size_t m = 0 ; m < NumVectors ; ++m) 
00803         for (size_t i = 0 ; i < NumMyRows_ ; ++i)
00804           y_ptr[m][i] = y2_ptr[m][i];
00805     }
00806   }
00807 
00808   ApplyFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
00809 }
00810 
00811 //==========================================================================
00812 template<class MatrixType>
00813 std::string Relaxation<MatrixType>::description() const {
00814   std::ostringstream oss;
00815   oss << Teuchos::Describable::description();
00816   if (isInitialized()) {
00817     if (isComputed()) {
00818       oss << "{status = initialized, computed";
00819     }
00820     else {
00821       oss << "{status = initialized, not computed";
00822     }
00823   }
00824   else {
00825     oss << "{status = not initialized, not computed";
00826   }
00827   //
00828   if (PrecType_ == Ifpack2::JACOBI)   oss << "Type = Jacobi, " << std::endl;
00829   else if (PrecType_ == Ifpack2::GS)  oss << "Type = Gauss-Seidel, " << std::endl;
00830   else if (PrecType_ == Ifpack2::SGS) oss << "Type = Sym. Gauss-Seidel, " << std::endl;
00831   //
00832   oss << ", global rows = " << A_->getGlobalNumRows()
00833       << ", global cols = " << A_->getGlobalNumCols();
00834 
00835   if (DoL1Method_)
00836     oss<<", using L1 correction with eta "<<L1Eta_;
00837 
00838   oss << "}";
00839   return oss.str();
00840 }
00841 
00842 //==========================================================================
00843 template<class MatrixType>
00844 void Relaxation<MatrixType>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00845   using std::endl;
00846   using std::setw;
00847   using Teuchos::VERB_DEFAULT;
00848   using Teuchos::VERB_NONE;
00849   using Teuchos::VERB_LOW;
00850   using Teuchos::VERB_MEDIUM;
00851   using Teuchos::VERB_HIGH;
00852   using Teuchos::VERB_EXTREME;
00853   Teuchos::EVerbosityLevel vl = verbLevel;
00854   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00855   const int myImageID = Comm_->getRank();
00856   Teuchos::OSTab tab(out);
00857 
00858   Scalar MinVal = Teuchos::ScalarTraits<Scalar>::zero();
00859   Scalar MaxVal = Teuchos::ScalarTraits<Scalar>::zero();
00860 
00861   if (IsComputed_) {
00862     Teuchos::ArrayRCP<Scalar> DiagView = Diagonal_->get1dViewNonConst();
00863     Scalar myMinVal = DiagView[0];
00864     Scalar myMaxVal = DiagView[0];
00865     for(typename Teuchos::ArrayRCP<Scalar>::size_type i=0; i<DiagView.size(); ++i) {
00866       if (Teuchos::ScalarTraits<Scalar>::magnitude(myMinVal) > Teuchos::ScalarTraits<Scalar>::magnitude(DiagView[i])) myMinVal = DiagView[i];
00867       if (Teuchos::ScalarTraits<Scalar>::magnitude(myMaxVal) < Teuchos::ScalarTraits<Scalar>::magnitude(DiagView[i])) myMaxVal = DiagView[i];
00868     }
00869 
00870     Teuchos::reduceAll(*Comm_, Teuchos::REDUCE_MIN, 1, &myMinVal, &MinVal);
00871     Teuchos::reduceAll(*Comm_, Teuchos::REDUCE_MAX, 1, &myMaxVal, &MaxVal);
00872   }
00873 
00874   //    none: print nothing
00875   //     low: print O(1) info from node 0
00876   //  medium: 
00877   //    high: 
00878   // extreme: 
00879   if (vl != VERB_NONE && myImageID == 0) {
00880     out << this->description() << endl;
00881     out << endl;
00882     out << "===============================================================================" << endl;
00883     out << "Sweeps         = " << NumSweeps_ << endl;
00884     out << "damping factor = " << DampingFactor_ << endl;
00885     if (PrecType_ == Ifpack2::GS && DoBackwardGS_) {
00886       out << "Using backward mode (GS only)" << endl;
00887     }
00888     if   (ZeroStartingSolution_) { out << "Using zero starting solution" << endl; }
00889     else                         { out << "Using input starting solution" << endl; }
00890     if   (Condest_ == -1.0) { out << "Condition number estimate       = N/A" << endl; }
00891     else                    { out << "Condition number estimate       = " << Condest_ << endl; }
00892     if (IsComputed_) {
00893       out << "Minimum value on stored diagonal = " << MinVal << endl;
00894       out << "Maximum value on stored diagonal = " << MaxVal << endl;
00895     }
00896     out << endl;
00897     out << "Phase           # calls    Total Time (s)     Total MFlops      MFlops/s       " << endl;
00898     out << "------------    -------    ---------------    ---------------   ---------------" << endl;
00899     out << setw(12) << "initialize()" << setw(5) << getNumInitialize() << "    " << setw(15) << getInitializeTime() << endl;
00900     out << setw(12) << "compute()" << setw(5) << getNumCompute()    << "    " << setw(15) << getComputeTime() << "    " 
00901         << setw(15) << getComputeFlops() << "    " 
00902         << setw(15) << (getComputeTime() != 0.0 ? getComputeFlops() / getComputeTime() * 1.0e-6 : 0.0) << endl;
00903     out << setw(12) << "apply()" << setw(5) << getNumApply()    << "    " << setw(15) << getApplyTime() << "    " 
00904         << setw(15) << getApplyFlops() << "    " 
00905         << setw(15) << (getApplyTime() != 0.0 ? getApplyFlops() / getApplyTime() * 1.0e-6 : 0.0) << endl;
00906     out << "===============================================================================" << endl;
00907     out << endl;
00908   }
00909 }
00910 
00911 }//namespace Ifpack2
00912 
00913 #endif // IFPACK2_RELAXATION_DEF_HPP
00914 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends