Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Details_Chebyshev_decl.hpp
Go to the documentation of this file.
00001 // ***********************************************************************
00002 // 
00003 //      Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
00004 //                 Copyright (2004) Sandia Corporation
00005 // 
00006 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00007 // license for use of this work by or on behalf of the U.S. Government.
00008 // 
00009 // This library is free software; you can redistribute it and/or modify
00010 // it under the terms of the GNU Lesser General Public License as
00011 // published by the Free Software Foundation; either version 2.1 of the
00012 // License, or (at your option) any later version.
00013 //  
00014 // This library is distributed in the hope that it will be useful, but
00015 // WITHOUT ANY WARRANTY; without even the implied warranty of
00016 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 // Lesser General Public License for more details.
00018 //  
00019 // You should have received a copy of the GNU Lesser General Public
00020 // License along with this library; if not, write to the Free Software
00021 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00022 // USA
00023 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00024 // 
00025 // ***********************************************************************
00026 
00027 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
00028 #define IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
00029 
00036 
00037 #include <Ifpack2_ConfigDefs.hpp>
00038 #include <Tpetra_CrsMatrix.hpp>
00039 
00040 namespace Ifpack2 {
00047 namespace Details {
00082 template<class ScalarType, class MV, class MAT>
00083 class Chebyshev {
00084 public:
00085   typedef ScalarType ST;
00086   typedef Teuchos::ScalarTraits<ScalarType> STS;
00087   typedef typename STS::magnitudeType MT;
00088   typedef Tpetra::Operator<typename MV::scalar_type,
00089          typename MV::local_ordinal_type,
00090          typename MV::global_ordinal_type,
00091          typename MV::node_type> OP;
00092   typedef Tpetra::Vector<typename MV::scalar_type,
00093        typename MV::local_ordinal_type,
00094        typename MV::global_ordinal_type,
00095        typename MV::node_type> V;
00096   typedef Tpetra::Map<typename MV::local_ordinal_type,
00097           typename MV::global_ordinal_type,
00098           typename MV::node_type> map_type;
00099 
00104   Chebyshev (Teuchos::RCP<const MAT> A);
00105 
00112   Chebyshev (Teuchos::RCP<const MAT> A, Teuchos::ParameterList& params);
00113 
00206   void setParameters (Teuchos::ParameterList& plist);
00207 
00223   void compute ();
00224 
00241   MT apply (const MV& B, MV& X);
00242 
00244   Teuchos::RCP<const MAT> getMatrix () const;
00245 
00247   bool hasTransposeApply () const;
00248 
00250   void print (std::ostream& out);
00251 
00252 private:
00254 
00255 
00256   Teuchos::RCP<const MAT> A_; 
00257 
00268   Teuchos::RCP<const V> D_;
00269 
00271 
00272 
00273 
00276   Teuchos::RCP<MV> V_;
00279   Teuchos::RCP<MV> W_;
00280 
00284   ST computedLambdaMax_;
00288   ST computedLambdaMin_;
00289 
00291 
00292 
00293 
00295   Teuchos::RCP<const V> userInvDiag_;
00296 
00299   ST lambdaMax_; 
00302   ST lambdaMin_; 
00305   ST eigRatio_;  
00309   ST minDiagVal_;
00311   int numIters_;
00313   int eigMaxIters_;
00315   bool zeroStartingSolution_;
00317   bool textbookAlgorithm_;
00319   bool computeMaxResNorm_;
00320 
00322 
00323 
00324 
00326   void checkConstructorInput () const;
00327 
00340   void
00341   makeTempMultiVectors (Teuchos::RCP<MV>& V1,
00342       Teuchos::RCP<MV>& W,
00343       const MV& X);
00344 
00346   static void 
00347   computeResidual (MV& R, const MV& B, const MAT& A, const MV& X, 
00348        const Teuchos::ETransp mode = Teuchos::NO_TRANS);
00349 
00355   static void solve (MV& Z, const V& D_inv, const MV& R);
00356 
00362   static void solve (MV& Z, const ST alpha, const V& D_inv, const MV& R);
00363 
00365   Teuchos::RCP<V> makeInverseDiagonal (const MAT& A) const;
00366 
00371   Teuchos::RCP<V> makeRangeMapVector (const V& D) const;
00372 
00391   void
00392   textbookApplyImpl (const MAT& A,
00393          const MV& B,
00394          MV& X,
00395          const int numIters,
00396          const ST lambdaMax,
00397          const ST lambdaMin,
00398          const ST eigRatio,
00399          const V& D_inv) const;
00400 
00423   void
00424   ifpackApplyImpl (const MAT& A,
00425        const MV& B,
00426        MV& X,
00427        const int numIters,
00428        const ST lambdaMax,
00429        const ST lambdaMin,
00430        const ST eigRatio,
00431        const V& D_inv);
00432 
00435   static ST
00436   powerMethod (const MAT& A, const V& D_inv, const int numIters);
00437 
00439   static MT maxNormInf (const MV& X);
00440 
00441   // mfh 22 Jan 2013: This code builds correctly, but does not
00442   // converge.  I'm leaving it in place in case someone else wants to
00443   // work on it.
00444 #if 0
00445 
00446 
00447 
00448 
00449 
00450 
00451 
00452 
00453 
00454 
00455 
00456 
00457 
00458 
00459 
00460 
00461 
00462 
00463 
00464 
00465 
00466 
00467   void
00468   mlApplyImpl (const MAT& A,
00469          const MV& B,
00470          MV& X,
00471          const int numIters,
00472          const ST lambdaMax,
00473          const ST lambdaMin,
00474          const ST eigRatio,
00475          const V& D_inv) 
00476   {
00477     const ST zero = Teuchos::as<ST> (0);
00478     const ST one = Teuchos::as<ST> (1);
00479     const ST two = Teuchos::as<ST> (2);
00480 
00481     MV pAux (B.getMap (), B.getNumVectors ()); // Result of A*X
00482     MV dk (B.getMap (), B.getNumVectors ()); // Solution update
00483     MV R (B.getMap (), B.getNumVectors ()); // Not in original ML; need for B - pAux
00484 
00485     ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
00486     ST alpha = lambdaMax / eigRatio;
00487 
00488     ST delta = (beta - alpha) / two;
00489     ST theta = (beta + alpha) / two;
00490     ST s1 = theta / delta;
00491     ST rhok = one / s1;
00492 
00493     // Diagonal: ML replaces entries containing 0 with 1.  We
00494     // shouldn't have any entries like that in typical test problems,
00495     // so it's OK not to do that here.
00496 
00497     // The (scaled) matrix is the identity: set X = D_inv * B.  (If A
00498     // is the identity, then certainly D_inv is too.  D_inv comes from
00499     // A, so if D_inv * A is the identity, then we still need to apply
00500     // the "preconditioner" D_inv to B as well, to get X.)
00501     if (lambdaMin == one && lambdaMin == lambdaMax) {
00502       solve (X, D_inv, B);
00503       return;
00504     }
00505 
00506     // The next bit of code is a direct translation of code from ML's
00507     // ML_Cheby function, in the "normal point scaling" section, which
00508     // is in lines 7365-7392 of ml_smoother.c.
00509 
00510     if (! zeroStartingSolution_) {
00511       // dk = (1/theta) * D_inv * (B - (A*X))
00512       A.apply (X, pAux); // pAux = A * X
00513       R = B;
00514       R.update (-one, pAux, one); // R = B - pAux
00515       dk.elementWiseMultiply (one/theta, D_inv, R, zero); // dk = (1/theta)*D_inv*R
00516       X.update (one, dk, one); // X = X + dk
00517     } else {
00518       dk.elementWiseMultiply (one/theta, D_inv, B, zero); // dk = (1/theta)*D_inv*B
00519       X = dk;
00520     }
00521 
00522     ST rhokp1, dtemp1, dtemp2;
00523     for (int k = 0; k < numIters-1; ++k) {
00524       A.apply (X, pAux);
00525       rhokp1 = one / (two*s1 - rhok);
00526       dtemp1 = rhokp1*rhok;
00527       dtemp2 = two*rhokp1/delta;
00528       rhok = rhokp1;
00529 
00530       R = B;
00531       R.update (-one, pAux, one); // R = B - pAux
00532       // dk = dtemp1 * dk + dtemp2 * D_inv * (B - pAux)
00533       dk.elementWiseMultiply (dtemp2, D_inv, B, dtemp1);
00534       X.update (one, dk, one); // X = X + dk
00535     }
00536   }
00537 #endif // 0
00538 
00539 };
00540 
00541 } // namespace Details
00542 } // namespace Ifpack2
00543 
00544 #endif // IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends