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