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