Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_Details_Chebyshev_decl.hpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
6 // Copyright (2009) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 // ***********************************************************************
45 //
46 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
47 // Copyright (2004) Sandia Corporation
48 //
49 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
50 // license for use of this work by or on behalf of the U.S. Government.
51 //
52 // This library is free software; you can redistribute it and/or modify
53 // it under the terms of the GNU Lesser General Public License as
54 // published by the Free Software Foundation; either version 2.1 of the
55 // License, or (at your option) any later version.
56 //
57 // This library is distributed in the hope that it will be useful, but
58 // WITHOUT ANY WARRANTY; without even the implied warranty of
59 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
60 // Lesser General Public License for more details.
61 //
62 // You should have received a copy of the GNU Lesser General Public
63 // License along with this library; if not, write to the Free Software
64 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
65 // USA
66 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
67 //
68 // ***********************************************************************
69 
70 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
71 #define IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
72 
79 
80 #include <Ifpack2_ConfigDefs.hpp>
81 #include <Teuchos_VerbosityLevel.hpp>
82 #include <Teuchos_Describable.hpp>
83 #include <Tpetra_CrsMatrix.hpp>
84 
85 namespace Ifpack2 {
86 
87 namespace Details {
126 template<class ScalarType, class MV>
127 class Chebyshev : public Teuchos::Describable {
128 public:
130 
131 
133  typedef ScalarType ST;
135  typedef Teuchos::ScalarTraits<ScalarType> STS;
137  typedef typename STS::magnitudeType MT;
139  typedef Tpetra::Operator<typename MV::scalar_type,
140  typename MV::local_ordinal_type,
141  typename MV::global_ordinal_type,
142  typename MV::node_type> op_type;
144  typedef Tpetra::RowMatrix<typename MV::scalar_type,
145  typename MV::local_ordinal_type,
146  typename MV::global_ordinal_type,
147  typename MV::node_type> row_matrix_type;
149  typedef Tpetra::Vector<typename MV::scalar_type,
150  typename MV::local_ordinal_type,
151  typename MV::global_ordinal_type,
152  typename MV::node_type> V;
154  typedef Tpetra::Map<typename MV::local_ordinal_type,
155  typename MV::global_ordinal_type,
156  typename MV::node_type> map_type;
158 
166  Chebyshev (Teuchos::RCP<const row_matrix_type> A);
167 
177  Chebyshev (Teuchos::RCP<const row_matrix_type> A, Teuchos::ParameterList& params);
178 
284  void setParameters (Teuchos::ParameterList& plist);
285 
319  void compute ();
320 
337  MT apply (const MV& B, MV& X);
338 
339  ST getLambdaMaxForApply() const;
340 
342  Teuchos::RCP<const row_matrix_type> getMatrix () const;
343 
349  void setMatrix (const Teuchos::RCP<const row_matrix_type>& A);
350 
352  bool hasTransposeApply () const;
353 
355  void print (std::ostream& out);
356 
358 
360 
362  std::string description() const;
363 
365  void
366  describe (Teuchos::FancyOStream& out,
367  const Teuchos::EVerbosityLevel verbLevel =
368  Teuchos::Describable::verbLevel_default) const;
370 private:
372 
373 
380  Teuchos::RCP<const row_matrix_type> A_;
381 
392  Teuchos::RCP<const V> D_;
393 
399  Teuchos::ArrayRCP<size_t> diagOffsets_;
400 
408  bool savedDiagOffsets_;
409 
411 
413 
417  Teuchos::RCP<MV> V_;
418 
422  Teuchos::RCP<MV> W_;
423 
429  ST computedLambdaMax_;
430 
436  ST computedLambdaMin_;
437 
439 
441 
444  ST lambdaMaxForApply_;
447  ST lambdaMinForApply_;
450  ST eigRatioForApply_;
451 
453 
455 
460  Teuchos::RCP<const V> userInvDiag_;
461 
465  ST userLambdaMax_;
466 
470  ST userLambdaMin_;
471 
475  ST userEigRatio_;
476 
481  ST minDiagVal_;
482 
484  int numIters_;
485 
487  int eigMaxIters_;
488 
490  bool zeroStartingSolution_;
491 
498  bool assumeMatrixUnchanged_;
499 
501  bool textbookAlgorithm_;
502 
504  bool computeMaxResNorm_;
505 
507 
509 
511  void checkConstructorInput () const;
512 
514  void checkInputMatrix () const;
515 
523  void reset ();
524 
537  void
538  makeTempMultiVectors (Teuchos::RCP<MV>& V1,
539  Teuchos::RCP<MV>& W,
540  const MV& X);
541 
543  static void
544  computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
545  const Teuchos::ETransp mode = Teuchos::NO_TRANS);
546 
552  static void solve (MV& Z, const V& D_inv, const MV& R);
553 
559  static void solve (MV& Z, const ST alpha, const V& D_inv, const MV& R);
560 
568  Teuchos::RCP<const V>
569  makeInverseDiagonal (const row_matrix_type& A, const bool useDiagOffsets=false) const;
570 
580  Teuchos::RCP<V> makeRangeMapVector (const Teuchos::RCP<V>& D) const;
581 
583  Teuchos::RCP<const V>
584  makeRangeMapVectorConst (const Teuchos::RCP<const V>& D) const;
585 
604  void
605  textbookApplyImpl (const op_type& A,
606  const MV& B,
607  MV& X,
608  const int numIters,
609  const ST lambdaMax,
610  const ST lambdaMin,
611  const ST eigRatio,
612  const V& D_inv) const;
613 
636  void
637  ifpackApplyImpl (const op_type& A,
638  const MV& B,
639  MV& X,
640  const int numIters,
641  const ST lambdaMax,
642  const ST lambdaMin,
643  const ST eigRatio,
644  const V& D_inv);
645 
658  static ST
659  powerMethodWithInitGuess (const op_type& A, const V& D_inv, const int numIters, V& x);
660 
670  static ST
671  powerMethod (const op_type& A, const V& D_inv, const int numIters);
672 
674  static MT maxNormInf (const MV& X);
675 
676  // mfh 22 Jan 2013: This code builds correctly, but does not
677  // converge. I'm leaving it in place in case someone else wants to
678  // work on it.
679 #if 0
680  void
703  mlApplyImpl (const op_type& A,
704  const MV& B,
705  MV& X,
706  const int numIters,
707  const ST lambdaMax,
708  const ST lambdaMin,
709  const ST eigRatio,
710  const V& D_inv)
711  {
712  const ST zero = Teuchos::as<ST> (0);
713  const ST one = Teuchos::as<ST> (1);
714  const ST two = Teuchos::as<ST> (2);
715 
716  MV pAux (B.getMap (), B.getNumVectors ()); // Result of A*X
717  MV dk (B.getMap (), B.getNumVectors ()); // Solution update
718  MV R (B.getMap (), B.getNumVectors ()); // Not in original ML; need for B - pAux
719 
720  ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
721  ST alpha = lambdaMax / eigRatio;
722 
723  ST delta = (beta - alpha) / two;
724  ST theta = (beta + alpha) / two;
725  ST s1 = theta / delta;
726  ST rhok = one / s1;
727 
728  // Diagonal: ML replaces entries containing 0 with 1. We
729  // shouldn't have any entries like that in typical test problems,
730  // so it's OK not to do that here.
731 
732  // The (scaled) matrix is the identity: set X = D_inv * B. (If A
733  // is the identity, then certainly D_inv is too. D_inv comes from
734  // A, so if D_inv * A is the identity, then we still need to apply
735  // the "preconditioner" D_inv to B as well, to get X.)
736  if (lambdaMin == one && lambdaMin == lambdaMax) {
737  solve (X, D_inv, B);
738  return;
739  }
740 
741  // The next bit of code is a direct translation of code from ML's
742  // ML_Cheby function, in the "normal point scaling" section, which
743  // is in lines 7365-7392 of ml_smoother.c.
744 
745  if (! zeroStartingSolution_) {
746  // dk = (1/theta) * D_inv * (B - (A*X))
747  A.apply (X, pAux); // pAux = A * X
748  R = B;
749  R.update (-one, pAux, one); // R = B - pAux
750  dk.elementWiseMultiply (one/theta, D_inv, R, zero); // dk = (1/theta)*D_inv*R
751  X.update (one, dk, one); // X = X + dk
752  } else {
753  dk.elementWiseMultiply (one/theta, D_inv, B, zero); // dk = (1/theta)*D_inv*B
754  X = dk;
755  }
756 
757  ST rhokp1, dtemp1, dtemp2;
758  for (int k = 0; k < numIters-1; ++k) {
759  A.apply (X, pAux);
760  rhokp1 = one / (two*s1 - rhok);
761  dtemp1 = rhokp1*rhok;
762  dtemp2 = two*rhokp1/delta;
763  rhok = rhokp1;
764 
765  R = B;
766  R.update (-one, pAux, one); // R = B - pAux
767  // dk = dtemp1 * dk + dtemp2 * D_inv * (B - pAux)
768  dk.elementWiseMultiply (dtemp2, D_inv, B, dtemp1);
769  X.update (one, dk, one); // X = X + dk
770  }
771  }
772 #endif // 0
773 
774 };
775 
776 } // namespace Details
777 } // namespace Ifpack2
778 
779 #endif // IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:316
MT apply(const MV &B, MV &X)
Solve Ax=b for x with Chebyshev iteration with left diagonal scaling.
Definition: Ifpack2_Details_Chebyshev_def.hpp:825
Tpetra::Operator< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > op_type
Specialization of Tpetra::Operator.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:142
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print a description of the Chebyshev solver to out.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1348
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:652
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1298
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:266
Ifpack2 implementation details.
Tpetra::Map< typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > map_type
Specialization of Tpetra::Map.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:156
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:872
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:137
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:133
Tpetra::RowMatrix< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > row_matrix_type
Specialization of Tpetra::RowMatrix.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:147
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:127
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1331
Teuchos::ScalarTraits< ScalarType > STS
Traits class for ST.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:135
bool hasTransposeApply() const
Whether it&#39;s possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1305
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
void compute()
(Re)compute the left scaling D_inv, and estimate min and max eigenvalues of D_inv * A...
Definition: Ifpack2_Details_Chebyshev_def.hpp:665
Tpetra::Vector< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > V
Specialization of Tpetra::Vector.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:152