Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Details_Chebyshev_decl.hpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00006 //                 Copyright (2009) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 // ***********************************************************************
00045 //
00046 //      Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
00047 //                 Copyright (2004) Sandia Corporation
00048 //
00049 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00050 // license for use of this work by or on behalf of the U.S. Government.
00051 //
00052 // This library is free software; you can redistribute it and/or modify
00053 // it under the terms of the GNU Lesser General Public License as
00054 // published by the Free Software Foundation; either version 2.1 of the
00055 // License, or (at your option) any later version.
00056 //
00057 // This library is distributed in the hope that it will be useful, but
00058 // WITHOUT ANY WARRANTY; without even the implied warranty of
00059 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00060 // Lesser General Public License for more details.
00061 //
00062 // You should have received a copy of the GNU Lesser General Public
00063 // License along with this library; if not, write to the Free Software
00064 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00065 // USA
00066 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00067 //
00068 // ***********************************************************************
00069 
00070 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
00071 #define IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
00072 
00079 
00080 #include <Ifpack2_ConfigDefs.hpp>
00081 #include <Teuchos_VerbosityLevel.hpp>
00082 #include <Teuchos_Describable.hpp>
00083 #include <Tpetra_CrsMatrix.hpp>
00084 
00085 namespace Ifpack2 {
00086 
00087 namespace Details {
00126 template<class ScalarType, class MV>
00127 class Chebyshev : public Teuchos::Describable {
00128 public:
00130 
00131 
00133   typedef ScalarType ST;
00135   typedef Teuchos::ScalarTraits<ScalarType> STS;
00137   typedef typename STS::magnitudeType MT;
00139   typedef Tpetra::Operator<typename MV::scalar_type,
00140                            typename MV::local_ordinal_type,
00141                            typename MV::global_ordinal_type,
00142                            typename MV::node_type> op_type;
00144   typedef Tpetra::RowMatrix<typename MV::scalar_type,
00145                            typename MV::local_ordinal_type,
00146                            typename MV::global_ordinal_type,
00147                            typename MV::node_type> row_matrix_type;
00149   typedef Tpetra::Vector<typename MV::scalar_type,
00150                          typename MV::local_ordinal_type,
00151                          typename MV::global_ordinal_type,
00152                          typename MV::node_type> V;
00154   typedef Tpetra::Map<typename MV::local_ordinal_type,
00155                       typename MV::global_ordinal_type,
00156                       typename MV::node_type> map_type;
00158 
00166   Chebyshev (Teuchos::RCP<const row_matrix_type> A);
00167 
00177   Chebyshev (Teuchos::RCP<const row_matrix_type> A, Teuchos::ParameterList& params);
00178 
00284   void setParameters (Teuchos::ParameterList& plist);
00285 
00319   void compute ();
00320 
00337   MT apply (const MV& B, MV& X);
00338 
00339   ST getLambdaMaxForApply() const;
00340 
00342   Teuchos::RCP<const row_matrix_type> getMatrix () const;
00343 
00349   void setMatrix (const Teuchos::RCP<const row_matrix_type>& A);
00350 
00352   bool hasTransposeApply () const;
00353 
00355   void print (std::ostream& out);
00356 
00358 
00359 
00360 
00362   std::string description() const;
00363 
00365   void
00366   describe (Teuchos::FancyOStream& out,
00367             const Teuchos::EVerbosityLevel verbLevel =
00368             Teuchos::Describable::verbLevel_default) const;
00370 private:
00372 
00373 
00380   Teuchos::RCP<const row_matrix_type> A_;
00381 
00392   Teuchos::RCP<const V> D_;
00393 
00399   Teuchos::ArrayRCP<size_t> diagOffsets_;
00400 
00408   bool savedDiagOffsets_;
00409 
00411 
00412 
00413 
00417   Teuchos::RCP<MV> V_;
00418 
00422   Teuchos::RCP<MV> W_;
00423 
00429   ST computedLambdaMax_;
00430 
00436   ST computedLambdaMin_;
00437 
00439 
00440 
00441 
00444   ST lambdaMaxForApply_;
00447   ST lambdaMinForApply_;
00450   ST eigRatioForApply_;
00451 
00453 
00454 
00455 
00460   Teuchos::RCP<const V> userInvDiag_;
00461 
00465   ST userLambdaMax_;
00466 
00470   ST userLambdaMin_;
00471 
00475   ST userEigRatio_;
00476 
00481   ST minDiagVal_;
00482 
00484   int numIters_;
00485 
00487   int eigMaxIters_;
00488 
00490   bool zeroStartingSolution_;
00491 
00498   bool assumeMatrixUnchanged_;
00499 
00501   bool textbookAlgorithm_;
00502 
00504   bool computeMaxResNorm_;
00505 
00507 
00508 
00509 
00511   void checkConstructorInput () const;
00512 
00514   void checkInputMatrix () const;
00515 
00523   void reset ();
00524 
00537   void
00538   makeTempMultiVectors (Teuchos::RCP<MV>& V1,
00539                         Teuchos::RCP<MV>& W,
00540                         const MV& X);
00541 
00543   static void
00544   computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
00545                    const Teuchos::ETransp mode = Teuchos::NO_TRANS);
00546 
00552   static void solve (MV& Z, const V& D_inv, const MV& R);
00553 
00559   static void solve (MV& Z, const ST alpha, const V& D_inv, const MV& R);
00560 
00568   Teuchos::RCP<const V>
00569   makeInverseDiagonal (const row_matrix_type& A, const bool useDiagOffsets=false) const;
00570 
00580   Teuchos::RCP<V> makeRangeMapVector (const Teuchos::RCP<V>& D) const;
00581 
00583   Teuchos::RCP<const V>
00584   makeRangeMapVectorConst (const Teuchos::RCP<const V>& D) const;
00585 
00604   void
00605   textbookApplyImpl (const op_type& A,
00606                      const MV& B,
00607                      MV& X,
00608                      const int numIters,
00609                      const ST lambdaMax,
00610                      const ST lambdaMin,
00611                      const ST eigRatio,
00612                      const V& D_inv) const;
00613 
00636   void
00637   ifpackApplyImpl (const op_type& A,
00638                    const MV& B,
00639                    MV& X,
00640                    const int numIters,
00641                    const ST lambdaMax,
00642                    const ST lambdaMin,
00643                    const ST eigRatio,
00644                    const V& D_inv);
00645 
00648   static ST
00649   powerMethod (const op_type& A, const V& D_inv, const int numIters);
00650 
00652   static MT maxNormInf (const MV& X);
00653 
00654   // mfh 22 Jan 2013: This code builds correctly, but does not
00655   // converge.  I'm leaving it in place in case someone else wants to
00656   // work on it.
00657 #if 0
00658 
00659 
00660 
00661 
00662 
00663 
00664 
00665 
00666 
00667 
00668 
00669 
00670 
00671 
00672 
00673 
00674 
00675 
00676 
00677 
00678 
00679 
00680   void
00681   mlApplyImpl (const op_type& A,
00682                const MV& B,
00683                MV& X,
00684                const int numIters,
00685                const ST lambdaMax,
00686                const ST lambdaMin,
00687                const ST eigRatio,
00688                const V& D_inv)
00689   {
00690     const ST zero = Teuchos::as<ST> (0);
00691     const ST one = Teuchos::as<ST> (1);
00692     const ST two = Teuchos::as<ST> (2);
00693 
00694     MV pAux (B.getMap (), B.getNumVectors ()); // Result of A*X
00695     MV dk (B.getMap (), B.getNumVectors ()); // Solution update
00696     MV R (B.getMap (), B.getNumVectors ()); // Not in original ML; need for B - pAux
00697 
00698     ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
00699     ST alpha = lambdaMax / eigRatio;
00700 
00701     ST delta = (beta - alpha) / two;
00702     ST theta = (beta + alpha) / two;
00703     ST s1 = theta / delta;
00704     ST rhok = one / s1;
00705 
00706     // Diagonal: ML replaces entries containing 0 with 1.  We
00707     // shouldn't have any entries like that in typical test problems,
00708     // so it's OK not to do that here.
00709 
00710     // The (scaled) matrix is the identity: set X = D_inv * B.  (If A
00711     // is the identity, then certainly D_inv is too.  D_inv comes from
00712     // A, so if D_inv * A is the identity, then we still need to apply
00713     // the "preconditioner" D_inv to B as well, to get X.)
00714     if (lambdaMin == one && lambdaMin == lambdaMax) {
00715       solve (X, D_inv, B);
00716       return;
00717     }
00718 
00719     // The next bit of code is a direct translation of code from ML's
00720     // ML_Cheby function, in the "normal point scaling" section, which
00721     // is in lines 7365-7392 of ml_smoother.c.
00722 
00723     if (! zeroStartingSolution_) {
00724       // dk = (1/theta) * D_inv * (B - (A*X))
00725       A.apply (X, pAux); // pAux = A * X
00726       R = B;
00727       R.update (-one, pAux, one); // R = B - pAux
00728       dk.elementWiseMultiply (one/theta, D_inv, R, zero); // dk = (1/theta)*D_inv*R
00729       X.update (one, dk, one); // X = X + dk
00730     } else {
00731       dk.elementWiseMultiply (one/theta, D_inv, B, zero); // dk = (1/theta)*D_inv*B
00732       X = dk;
00733     }
00734 
00735     ST rhokp1, dtemp1, dtemp2;
00736     for (int k = 0; k < numIters-1; ++k) {
00737       A.apply (X, pAux);
00738       rhokp1 = one / (two*s1 - rhok);
00739       dtemp1 = rhokp1*rhok;
00740       dtemp2 = two*rhokp1/delta;
00741       rhok = rhokp1;
00742 
00743       R = B;
00744       R.update (-one, pAux, one); // R = B - pAux
00745       // dk = dtemp1 * dk + dtemp2 * D_inv * (B - pAux)
00746       dk.elementWiseMultiply (dtemp2, D_inv, B, dtemp1);
00747       X.update (one, dk, one); // X = X + dk
00748     }
00749   }
00750 #endif // 0
00751 
00752 };
00753 
00754 } // namespace Details
00755 } // namespace Ifpack2
00756 
00757 #endif // IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends