Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Details_Chebyshev_decl.hpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ***********************************************************************
4//
5// Ifpack2: Templated 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#ifndef IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
45#define IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
46
53
54#include "Ifpack2_ConfigDefs.hpp"
55#include "Teuchos_VerbosityLevel.hpp"
56#include "Teuchos_Describable.hpp"
57#include "Tpetra_CrsMatrix.hpp"
58
59namespace Ifpack2 {
60namespace Details {
61
62#ifndef DOXYGEN_SHOULD_SKIP_THIS
63template<class TpetraOperatorType>
64class ChebyshevKernel; // forward declaration
65#endif // DOXYGEN_SHOULD_SKIP_THIS
66
105template<class ScalarType, class MV>
106class Chebyshev : public Teuchos::Describable {
107public:
109
110
112 typedef ScalarType ST;
114 typedef Teuchos::ScalarTraits<ScalarType> STS;
116 typedef typename STS::magnitudeType MT;
118 typedef Tpetra::Operator<typename MV::scalar_type,
119 typename MV::local_ordinal_type,
120 typename MV::global_ordinal_type,
121 typename MV::node_type> op_type;
123 typedef Tpetra::RowMatrix<typename MV::scalar_type,
124 typename MV::local_ordinal_type,
125 typename MV::global_ordinal_type,
126 typename MV::node_type> row_matrix_type;
128 typedef Tpetra::Vector<typename MV::scalar_type,
129 typename MV::local_ordinal_type,
130 typename MV::global_ordinal_type,
131 typename MV::node_type> V;
133 typedef Tpetra::Map<typename MV::local_ordinal_type,
134 typename MV::global_ordinal_type,
135 typename MV::node_type> map_type;
137
145 Chebyshev (Teuchos::RCP<const row_matrix_type> A);
146
156 Chebyshev (Teuchos::RCP<const row_matrix_type> A, Teuchos::ParameterList& params);
157
266 void setParameters (Teuchos::ParameterList& plist);
267
268 void setZeroStartingSolution (bool zeroStartingSolution) { zeroStartingSolution_ = zeroStartingSolution; }
269
303 void compute ();
304
321 MT apply (const MV& B, MV& X);
322
323 ST getLambdaMaxForApply() const;
324
326 Teuchos::RCP<const row_matrix_type> getMatrix () const;
327
333 void setMatrix (const Teuchos::RCP<const row_matrix_type>& A);
334
336 bool hasTransposeApply () const;
337
339 void print (std::ostream& out);
340
342
344
346 std::string description() const;
347
349 void
350 describe (Teuchos::FancyOStream& out,
351 const Teuchos::EVerbosityLevel verbLevel =
352 Teuchos::Describable::verbLevel_default) const;
354private:
356
357
364 Teuchos::RCP<const row_matrix_type> A_;
365
367 Teuchos::RCP<ChebyshevKernel<op_type> > ck_;
368
379 Teuchos::RCP<const V> D_;
380
382 typedef Kokkos::View<size_t*, typename MV::node_type::device_type> offsets_type;
383
389 offsets_type diagOffsets_;
390
398 bool savedDiagOffsets_;
399
401
403
407 Teuchos::RCP<MV> W_;
408
414 ST computedLambdaMax_;
415
421 ST computedLambdaMin_;
422
424
426
429 ST lambdaMaxForApply_;
430
437 MT boostFactor_;
440 ST lambdaMinForApply_;
443 ST eigRatioForApply_;
444
446
448
453 Teuchos::RCP<const V> userInvDiag_;
454
458 ST userLambdaMax_;
459
463 ST userLambdaMin_;
464
468 ST userEigRatio_;
469
474 ST minDiagVal_;
475
477 int numIters_;
478
480 int eigMaxIters_;
481
483 MT eigRelTolerance_;
484
486 bool eigKeepVectors_;
487
489 std::string eigenAnalysisType_;
490
492 Teuchos::RCP<V> eigVector_;
493 Teuchos::RCP<V> eigVector2_;
494
496 int eigNormalizationFreq_;
497
499 bool zeroStartingSolution_;
500
507 bool assumeMatrixUnchanged_;
508
510 bool textbookAlgorithm_;
511
513 bool computeMaxResNorm_;
514
516 bool computeSpectralRadius_;
517
520 bool ckUseNativeSpMV_;
521
525 Teuchos::RCP<Teuchos::FancyOStream> out_;
526
528 bool debug_;
529
531
533
535 void checkConstructorInput () const;
536
538 void checkInputMatrix () const;
539
547 void reset ();
548
554 Teuchos::RCP<MV> makeTempMultiVector (const MV& B);
555
557 void
558 firstIterationWithZeroStartingSolution
559 (MV& W,
560 const ScalarType& alpha,
561 const V& D_inv,
562 const MV& B,
563 MV& X);
564
566 static void
567 computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
568 const Teuchos::ETransp mode = Teuchos::NO_TRANS);
569
575 static void solve (MV& Z, const V& D_inv, const MV& R);
576
582 static void solve (MV& Z, const ST alpha, const V& D_inv, const MV& R);
583
591 Teuchos::RCP<const V>
592 makeInverseDiagonal (const row_matrix_type& A, const bool useDiagOffsets=false) const;
593
603 Teuchos::RCP<V> makeRangeMapVector (const Teuchos::RCP<V>& D) const;
604
606 Teuchos::RCP<const V>
607 makeRangeMapVectorConst (const Teuchos::RCP<const V>& D) const;
608
627 void
628 textbookApplyImpl (const op_type& A,
629 const MV& B,
630 MV& X,
631 const int numIters,
632 const ST lambdaMax,
633 const ST lambdaMin,
634 const ST eigRatio,
635 const V& D_inv) const;
636
659 void
660 ifpackApplyImpl (const op_type& A,
661 const MV& B,
662 MV& X,
663 const int numIters,
664 const ST lambdaMax,
665 const ST lambdaMin,
666 const ST eigRatio,
667 const V& D_inv);
668
681 ST
682 cgMethodWithInitGuess (const op_type& A, const V& D_inv, const int numIters, V& x);
683
693 ST
694 cgMethod (const op_type& A, const V& D_inv, const int numIters);
695
697 static MT maxNormInf (const MV& X);
698
699 // mfh 22 Jan 2013: This code builds correctly, but does not
700 // converge. I'm leaving it in place in case someone else wants to
701 // work on it.
702#if 0
725 void
726 mlApplyImpl (const op_type& A,
727 const MV& B,
728 MV& X,
729 const int numIters,
730 const ST lambdaMax,
731 const ST lambdaMin,
732 const ST eigRatio,
733 const V& D_inv)
734 {
735 const ST zero = Teuchos::as<ST> (0);
736 const ST one = Teuchos::as<ST> (1);
737 const ST two = Teuchos::as<ST> (2);
738
739 MV pAux (B.getMap (), B.getNumVectors ()); // Result of A*X
740 MV dk (B.getMap (), B.getNumVectors ()); // Solution update
741 MV R (B.getMap (), B.getNumVectors ()); // Not in original ML; need for B - pAux
742
743 ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
744 ST alpha = lambdaMax / eigRatio;
745
746 ST delta = (beta - alpha) / two;
747 ST theta = (beta + alpha) / two;
748 ST s1 = theta / delta;
749 ST rhok = one / s1;
750
751 // Diagonal: ML replaces entries containing 0 with 1. We
752 // shouldn't have any entries like that in typical test problems,
753 // so it's OK not to do that here.
754
755 // The (scaled) matrix is the identity: set X = D_inv * B. (If A
756 // is the identity, then certainly D_inv is too. D_inv comes from
757 // A, so if D_inv * A is the identity, then we still need to apply
758 // the "preconditioner" D_inv to B as well, to get X.)
759 if (lambdaMin == one && lambdaMin == lambdaMax) {
760 solve (X, D_inv, B);
761 return;
762 }
763
764 // The next bit of code is a direct translation of code from ML's
765 // ML_Cheby function, in the "normal point scaling" section, which
766 // is in lines 7365-7392 of ml_smoother.c.
767
768 if (! zeroStartingSolution_) {
769 // dk = (1/theta) * D_inv * (B - (A*X))
770 A.apply (X, pAux); // pAux = A * X
771 R = B;
772 R.update (-one, pAux, one); // R = B - pAux
773 dk.elementWiseMultiply (one/theta, D_inv, R, zero); // dk = (1/theta)*D_inv*R
774 X.update (one, dk, one); // X = X + dk
775 } else {
776 dk.elementWiseMultiply (one/theta, D_inv, B, zero); // dk = (1/theta)*D_inv*B
777 X = dk;
778 }
779
780 ST rhokp1, dtemp1, dtemp2;
781 for (int k = 0; k < numIters-1; ++k) {
782 A.apply (X, pAux);
783 rhokp1 = one / (two*s1 - rhok);
784 dtemp1 = rhokp1*rhok;
785 dtemp2 = two*rhokp1/delta;
786 rhok = rhokp1;
787
788 R = B;
789 R.update (-one, pAux, one); // R = B - pAux
790 // dk = dtemp1 * dk + dtemp2 * D_inv * (B - pAux)
791 dk.elementWiseMultiply (dtemp2, D_inv, B, dtemp1);
792 X.update (one, dk, one); // X = X + dk
793 }
794 }
795#endif // 0
797};
798
799} // namespace Details
800} // namespace Ifpack2
801
802#endif // IFPACK2_DETAILS_CHEBYSHEV_DECL_HPP
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:106
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:131
Teuchos::ScalarTraits< ScalarType > STS
Traits class for ST.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:114
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:796
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:353
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:755
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1560
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:135
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:126
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:1580
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1042
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:121
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:112
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:116
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:989
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1536
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1529
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74