Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Details_Chebyshev_def.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_DEF_HPP
45#define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
46
53
56// #include "Ifpack2_Details_ScaledDampedResidual.hpp"
57#include "Ifpack2_Details_ChebyshevKernel.hpp"
58#include "Kokkos_ArithTraits.hpp"
59#include "Teuchos_FancyOStream.hpp"
60#include "Teuchos_oblackholestream.hpp"
61#include "Tpetra_Details_residual.hpp"
62#include "Teuchos_LAPACK.hpp"
63#include "Ifpack2_Details_LapackSupportsScalar.hpp"
64#include <cmath>
65#include <iostream>
66
67namespace Ifpack2 {
68namespace Details {
69
70namespace { // (anonymous)
71
72// We use this text a lot in error messages.
73const char computeBeforeApplyReminder[] =
74 "This means one of the following:\n"
75 " - you have not yet called compute() on this instance, or \n"
76 " - you didn't call compute() after calling setParameters().\n\n"
77 "After creating an Ifpack2::Chebyshev instance,\n"
78 "you must _always_ call compute() at least once before calling apply().\n"
79 "After calling compute() once, you do not need to call it again,\n"
80 "unless the matrix has changed or you have changed parameters\n"
81 "(by calling setParameters()).";
82
83} // namespace (anonymous)
84
85// ReciprocalThreshold stuff below needs to be in a namspace visible outside
86// of this file
87template<class XV, class SizeType = typename XV::size_type>
88struct V_ReciprocalThresholdSelfFunctor
89{
90 typedef typename XV::execution_space execution_space;
91 typedef typename XV::non_const_value_type value_type;
92 typedef SizeType size_type;
93 typedef Kokkos::Details::ArithTraits<value_type> KAT;
94 typedef typename KAT::mag_type mag_type;
95
96 XV X_;
97 const value_type minVal_;
98 const mag_type minValMag_;
99
100 V_ReciprocalThresholdSelfFunctor (const XV& X,
101 const value_type& min_val) :
102 X_ (X),
103 minVal_ (min_val),
104 minValMag_ (KAT::abs (min_val))
105 {}
106
107 KOKKOS_INLINE_FUNCTION
108 void operator() (const size_type& i) const
109 {
110 const mag_type X_i_abs = KAT::abs (X_[i]);
111
112 if (X_i_abs < minValMag_) {
113 X_[i] = minVal_;
114 }
115 else {
116 X_[i] = KAT::one () / X_[i];
117 }
118 }
119};
120
121template<class XV, class SizeType = typename XV::size_type>
122struct LocalReciprocalThreshold {
123 static void
124 compute (const XV& X,
125 const typename XV::non_const_value_type& minVal)
126 {
127 typedef typename XV::execution_space execution_space;
128 Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.extent (0));
129 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
130 Kokkos::parallel_for (policy, op);
131 }
132};
133
134template <class TpetraVectorType,
135 const bool classic = TpetraVectorType::node_type::classic>
136struct GlobalReciprocalThreshold {};
137
138template <class TpetraVectorType>
139struct GlobalReciprocalThreshold<TpetraVectorType, true> {
140 static void
141 compute (TpetraVectorType& V,
142 const typename TpetraVectorType::scalar_type& min_val)
143 {
144 typedef typename TpetraVectorType::scalar_type scalar_type;
145 typedef typename TpetraVectorType::mag_type mag_type;
146 typedef Kokkos::Details::ArithTraits<scalar_type> STS;
147
148 const scalar_type ONE = STS::one ();
149 const mag_type min_val_abs = STS::abs (min_val);
150
151 Teuchos::ArrayRCP<scalar_type> V_0 = V.getDataNonConst (0);
152 const size_t lclNumRows = V.getLocalLength ();
153
154 for (size_t i = 0; i < lclNumRows; ++i) {
155 const scalar_type V_0i = V_0[i];
156 if (STS::abs (V_0i) < min_val_abs) {
157 V_0[i] = min_val;
158 } else {
159 V_0[i] = ONE / V_0i;
160 }
161 }
162 }
163};
164
165template <class TpetraVectorType>
166struct GlobalReciprocalThreshold<TpetraVectorType, false> {
167 static void
168 compute (TpetraVectorType& X,
169 const typename TpetraVectorType::scalar_type& minVal)
170 {
171 typedef typename TpetraVectorType::impl_scalar_type value_type;
172
173 const value_type minValS = static_cast<value_type> (minVal);
174 auto X_0 = Kokkos::subview (X.getLocalViewDevice (Tpetra::Access::ReadWrite),
175 Kokkos::ALL (), 0);
176 LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
177 }
178};
179
180// Utility function for inverting diagonal with threshold.
181template <typename S, typename L, typename G, typename N>
182void
183reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V, const S& minVal)
184{
185 GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
186}
187
188
189template<class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
190struct LapackHelper{
191 static
192 ScalarType
193 tri_diag_spectral_radius(Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
194 Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
195 throw std::runtime_error("LAPACK does not support the scalar type.");
196 }
197};
198
199
200template<class ScalarType>
201struct LapackHelper<ScalarType,true> {
202 static
203 ScalarType
204 tri_diag_spectral_radius(Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
205 Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
206 using STS = Teuchos::ScalarTraits<ScalarType>;
207 using MagnitudeType = typename STS::magnitudeType;
208 int info = 0;
209 const int N = diag.size ();
210 ScalarType scalar_dummy;
211 std::vector<MagnitudeType> mag_dummy(4*N);
212 char char_N = 'N';
213
214 // lambdaMin = one;
215 ScalarType lambdaMax = STS::one();
216 if( N > 2 ) {
217 Teuchos::LAPACK<int,ScalarType> lapack;
218 lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
219 &scalar_dummy, 1, &mag_dummy[0], &info);
220 TEUCHOS_TEST_FOR_EXCEPTION
221 (info < 0, std::logic_error, "Ifpack2::Details::LapackHelper::tri_diag_spectral_radius:"
222 "LAPACK's _PTEQR failed with info = "
223 << info << " < 0. This suggests there might be a bug in the way Ifpack2 "
224 "is calling LAPACK. Please report this to the Ifpack2 developers.");
225 // lambdaMin = Teuchos::as<ScalarType> (diag[N-1]);
226 lambdaMax = Teuchos::as<ScalarType> (diag[0]);
227 }
228 return lambdaMax;
229 }
230};
231
232
233template<class ScalarType, class MV>
234void Chebyshev<ScalarType, MV>::checkInputMatrix () const
235{
236 TEUCHOS_TEST_FOR_EXCEPTION(
237 ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
238 std::invalid_argument,
239 "Ifpack2::Chebyshev: The input matrix A must be square. "
240 "A has " << A_->getGlobalNumRows () << " rows and "
241 << A_->getGlobalNumCols () << " columns.");
242
243 // In debug mode, test that the domain and range Maps of the matrix
244 // are the same.
245 if (debug_ && ! A_.is_null ()) {
246 Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
247 Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
248
249 // isSameAs is a collective, but if the two pointers are the same,
250 // isSameAs will assume that they are the same on all processes, and
251 // return true without an all-reduce.
252 TEUCHOS_TEST_FOR_EXCEPTION(
253 ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
254 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
255 "the same (in the sense of isSameAs())." << std::endl << "We only check "
256 "for this in debug mode.");
257 }
258}
259
260
261template<class ScalarType, class MV>
262void
263Chebyshev<ScalarType, MV>::
264checkConstructorInput () const
265{
266 // mfh 12 Aug 2016: The if statement avoids an "unreachable
267 // statement" warning for the checkInputMatrix() call, when
268 // STS::isComplex is false.
269 if (STS::isComplex) {
270 TEUCHOS_TEST_FOR_EXCEPTION
271 (true, std::logic_error, "Ifpack2::Chebyshev: This class' implementation "
272 "of Chebyshev iteration only works for real-valued, symmetric positive "
273 "definite matrices. However, you instantiated this class for ScalarType"
274 " = " << Teuchos::TypeNameTraits<ScalarType>::name () << ", which is a "
275 "complex-valued type. While this may be algorithmically correct if all "
276 "of the complex numbers in the matrix have zero imaginary part, we "
277 "forbid using complex ScalarType altogether in order to remind you of "
278 "the limitations of our implementation (and of the algorithm itself).");
279 }
280 else {
281 checkInputMatrix ();
282 }
283}
284
285template<class ScalarType, class MV>
287Chebyshev (Teuchos::RCP<const row_matrix_type> A) :
288 A_ (A),
289 savedDiagOffsets_ (false),
290 computedLambdaMax_ (STS::nan ()),
291 computedLambdaMin_ (STS::nan ()),
292 lambdaMaxForApply_ (STS::nan ()),
293 lambdaMinForApply_ (STS::nan ()),
294 eigRatioForApply_ (STS::nan ()),
295 userLambdaMax_ (STS::nan ()),
296 userLambdaMin_ (STS::nan ()),
297 userEigRatio_ (Teuchos::as<ST> (30)),
298 minDiagVal_ (STS::eps ()),
299 numIters_ (1),
300 eigMaxIters_ (10),
301 eigRelTolerance_(Teuchos::ScalarTraits<MT>::zero ()),
302 eigKeepVectors_(false),
303 eigenAnalysisType_("power method"),
304 eigNormalizationFreq_(1),
305 zeroStartingSolution_ (true),
306 assumeMatrixUnchanged_ (false),
307 textbookAlgorithm_ (false),
308 computeMaxResNorm_ (false),
309 computeSpectralRadius_(true),
310 ckUseNativeSpMV_(false),
311 debug_ (false)
312{
313 checkConstructorInput ();
314}
315
316template<class ScalarType, class MV>
318Chebyshev (Teuchos::RCP<const row_matrix_type> A,
319 Teuchos::ParameterList& params) :
320 A_ (A),
321 savedDiagOffsets_ (false),
322 computedLambdaMax_ (STS::nan ()),
323 computedLambdaMin_ (STS::nan ()),
324 lambdaMaxForApply_ (STS::nan ()),
325 boostFactor_ (static_cast<MT> (1.1)),
326 lambdaMinForApply_ (STS::nan ()),
327 eigRatioForApply_ (STS::nan ()),
328 userLambdaMax_ (STS::nan ()),
329 userLambdaMin_ (STS::nan ()),
330 userEigRatio_ (Teuchos::as<ST> (30)),
331 minDiagVal_ (STS::eps ()),
332 numIters_ (1),
333 eigMaxIters_ (10),
334 eigRelTolerance_(Teuchos::ScalarTraits<MT>::zero ()),
335 eigKeepVectors_(false),
336 eigenAnalysisType_("power method"),
337 eigNormalizationFreq_(1),
338 zeroStartingSolution_ (true),
339 assumeMatrixUnchanged_ (false),
340 textbookAlgorithm_ (false),
341 computeMaxResNorm_ (false),
342 computeSpectralRadius_(true),
343 ckUseNativeSpMV_(false),
344 debug_ (false)
345{
346 checkConstructorInput ();
347 setParameters (params);
348}
349
350template<class ScalarType, class MV>
351void
353setParameters (Teuchos::ParameterList& plist)
354{
355 using Teuchos::RCP;
356 using Teuchos::rcp;
357 using Teuchos::rcp_const_cast;
358
359 // Note to developers: The logic for this method is complicated,
360 // because we want to accept Ifpack and ML parameters whenever
361 // possible, but we don't want to add their default values to the
362 // user's ParameterList. That's why we do all the isParameter()
363 // checks, instead of using the two-argument version of get()
364 // everywhere. The min and max eigenvalue parameters are also a
365 // special case, because we decide whether or not to do eigenvalue
366 // analysis based on whether the user supplied the max eigenvalue.
367
368 // Default values of all the parameters.
369 const ST defaultLambdaMax = STS::nan ();
370 const ST defaultLambdaMin = STS::nan ();
371 // 30 is Ifpack::Chebyshev's default. ML has a different default
372 // eigRatio for smoothers and the coarse-grid solve (if using
373 // Chebyshev for that). The former uses 20; the latter uses 30.
374 // We're testing smoothers here, so use 20. (However, if you give
375 // ML an Epetra matrix, it will use Ifpack for Chebyshev, in which
376 // case it would defer to Ifpack's default settings.)
377 const ST defaultEigRatio = Teuchos::as<ST> (30);
378 const MT defaultBoostFactor = static_cast<MT> (1.1);
379 const ST defaultMinDiagVal = STS::eps ();
380 const int defaultNumIters = 1;
381 const int defaultEigMaxIters = 10;
382 const MT defaultEigRelTolerance = Teuchos::ScalarTraits<MT>::zero ();
383 const bool defaultEigKeepVectors = false;
384 const int defaultEigNormalizationFreq = 1;
385 const bool defaultZeroStartingSolution = true; // Ifpack::Chebyshev default
386 const bool defaultAssumeMatrixUnchanged = false;
387 const bool defaultTextbookAlgorithm = false;
388 const bool defaultComputeMaxResNorm = false;
389 const bool defaultComputeSpectralRadius = true;
390 const bool defaultCkUseNativeSpMV = false;
391 const bool defaultDebug = false;
392
393 // We'll set the instance data transactionally, after all reads
394 // from the ParameterList. That way, if any of the ParameterList
395 // reads fail (e.g., due to the wrong parameter type), we will not
396 // have left the instance data in a half-changed state.
397 RCP<const V> userInvDiagCopy; // if nonnull: deep copy of user's Vector
398 ST lambdaMax = defaultLambdaMax;
399 ST lambdaMin = defaultLambdaMin;
400 ST eigRatio = defaultEigRatio;
401 MT boostFactor = defaultBoostFactor;
402 ST minDiagVal = defaultMinDiagVal;
403 int numIters = defaultNumIters;
404 int eigMaxIters = defaultEigMaxIters;
405 MT eigRelTolerance = defaultEigRelTolerance;
406 bool eigKeepVectors = defaultEigKeepVectors;
407 int eigNormalizationFreq = defaultEigNormalizationFreq;
408 bool zeroStartingSolution = defaultZeroStartingSolution;
409 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
410 bool textbookAlgorithm = defaultTextbookAlgorithm;
411 bool computeMaxResNorm = defaultComputeMaxResNorm;
412 bool computeSpectralRadius = defaultComputeSpectralRadius;
413 bool ckUseNativeSpMV = defaultCkUseNativeSpMV;
414 bool debug = defaultDebug;
415
416 // Fetch the parameters from the ParameterList. Defer all
417 // externally visible side effects until we have finished all
418 // ParameterList interaction. This makes the method satisfy the
419 // strong exception guarantee.
420
421 if (plist.isType<bool> ("debug")) {
422 debug = plist.get<bool> ("debug");
423 }
424 else if (plist.isType<int> ("debug")) {
425 const int debugInt = plist.get<bool> ("debug");
426 debug = debugInt != 0;
427 }
428
429 // Get the user-supplied inverse diagonal.
430 //
431 // Check for a raw pointer (const V* or V*), for Ifpack
432 // compatibility, as well as for RCP<const V>, RCP<V>, const V, or
433 // V. We'll make a deep copy of the vector at the end of this
434 // method anyway, so its const-ness doesn't matter. We handle the
435 // latter two cases ("const V" or "V") specially (copy them into
436 // userInvDiagCopy first, which is otherwise null at the end of the
437 // long if-then chain) to avoid an extra copy.
438
439 const char opInvDiagLabel[] = "chebyshev: operator inv diagonal";
440 if (plist.isParameter (opInvDiagLabel)) {
441 // Pointer to the user's Vector, if provided.
442 RCP<const V> userInvDiag;
443
444 if (plist.isType<const V*> (opInvDiagLabel)) {
445 const V* rawUserInvDiag =
446 plist.get<const V*> (opInvDiagLabel);
447 // Nonowning reference (we'll make a deep copy below)
448 userInvDiag = rcp (rawUserInvDiag, false);
449 }
450 else if (plist.isType<const V*> (opInvDiagLabel)) {
451 V* rawUserInvDiag = plist.get<V*> (opInvDiagLabel);
452 // Nonowning reference (we'll make a deep copy below)
453 userInvDiag = rcp (const_cast<const V*> (rawUserInvDiag), false);
454 }
455 else if (plist.isType<RCP<const V>> (opInvDiagLabel)) {
456 userInvDiag = plist.get<RCP<const V> > (opInvDiagLabel);
457 }
458 else if (plist.isType<RCP<V>> (opInvDiagLabel)) {
459 RCP<V> userInvDiagNonConst =
460 plist.get<RCP<V> > (opInvDiagLabel);
461 userInvDiag = rcp_const_cast<const V> (userInvDiagNonConst);
462 }
463 else if (plist.isType<const V> (opInvDiagLabel)) {
464 const V& userInvDiagRef = plist.get<const V> (opInvDiagLabel);
465 userInvDiagCopy = rcp (new V (userInvDiagRef, Teuchos::Copy));
466 userInvDiag = userInvDiagCopy;
467 }
468 else if (plist.isType<V> (opInvDiagLabel)) {
469 V& userInvDiagNonConstRef = plist.get<V> (opInvDiagLabel);
470 const V& userInvDiagRef = const_cast<const V&> (userInvDiagNonConstRef);
471 userInvDiagCopy = rcp (new V (userInvDiagRef, Teuchos::Copy));
472 userInvDiag = userInvDiagCopy;
473 }
474
475 // NOTE: If the user's parameter has some strange type that we
476 // didn't test above, userInvDiag might still be null. You may
477 // want to add an error test for this condition. Currently, we
478 // just say in this case that the user didn't give us a Vector.
479
480 // If we have userInvDiag but don't have a deep copy yet, make a
481 // deep copy now.
482 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
483 userInvDiagCopy = rcp (new V (*userInvDiag, Teuchos::Copy));
484 }
485
486 // NOTE: userInvDiag, if provided, is a row Map version of the
487 // Vector. We don't necessarily have a range Map yet. compute()
488 // would be the proper place to compute the range Map version of
489 // userInvDiag.
490 }
491
492 // Load the kernel fuse override from the parameter list
493 if (plist.isParameter ("chebyshev: use native spmv"))
494 ckUseNativeSpMV = plist.get("chebyshev: use native spmv", ckUseNativeSpMV);
495
496 // Don't fill in defaults for the max or min eigenvalue, because
497 // this class uses the existence of those parameters to determine
498 // whether it should do eigenanalysis.
499 if (plist.isParameter ("chebyshev: max eigenvalue")) {
500 if (plist.isType<double>("chebyshev: max eigenvalue"))
501 lambdaMax = plist.get<double> ("chebyshev: max eigenvalue");
502 else
503 lambdaMax = plist.get<ST> ("chebyshev: max eigenvalue");
504 TEUCHOS_TEST_FOR_EXCEPTION(
505 STS::isnaninf (lambdaMax), std::invalid_argument,
506 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
507 "parameter is NaN or Inf. This parameter is optional, but if you "
508 "choose to supply it, it must have a finite value.");
509 }
510 if (plist.isParameter ("chebyshev: min eigenvalue")) {
511 if (plist.isType<double>("chebyshev: min eigenvalue"))
512 lambdaMin = plist.get<double> ("chebyshev: min eigenvalue");
513 else
514 lambdaMin = plist.get<ST> ("chebyshev: min eigenvalue");
515 TEUCHOS_TEST_FOR_EXCEPTION(
516 STS::isnaninf (lambdaMin), std::invalid_argument,
517 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
518 "parameter is NaN or Inf. This parameter is optional, but if you "
519 "choose to supply it, it must have a finite value.");
520 }
521
522 // Only fill in Ifpack2's name for the default parameter, not ML's.
523 if (plist.isParameter ("smoother: Chebyshev alpha")) { // ML compatibility
524 if (plist.isType<double>("smoother: Chebyshev alpha"))
525 eigRatio = plist.get<double> ("smoother: Chebyshev alpha");
526 else
527 eigRatio = plist.get<ST> ("smoother: Chebyshev alpha");
528 }
529 // Ifpack2's name overrides ML's name.
530 eigRatio = plist.get ("chebyshev: ratio eigenvalue", eigRatio);
531 TEUCHOS_TEST_FOR_EXCEPTION(
532 STS::isnaninf (eigRatio), std::invalid_argument,
533 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
534 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
535 "This parameter is optional, but if you choose to supply it, it must have "
536 "a finite value.");
537 // mfh 11 Feb 2013: This class is currently only correct for real
538 // Scalar types, but we still want it to build for complex Scalar
539 // type so that users of Ifpack2::Factory can build their
540 // executables for real or complex Scalar type. Thus, we take the
541 // real parts here, which are always less-than comparable.
542 TEUCHOS_TEST_FOR_EXCEPTION(
543 STS::real (eigRatio) < STS::real (STS::one ()),
544 std::invalid_argument,
545 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
546 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
547 "but you supplied the value " << eigRatio << ".");
548
549 // See Github Issue #234. This parameter may be either MT
550 // (preferred) or double. We check both.
551 {
552 const char paramName[] = "chebyshev: boost factor";
553
554 if (plist.isParameter (paramName)) {
555 if (plist.isType<MT> (paramName)) { // MT preferred
556 boostFactor = plist.get<MT> (paramName);
557 }
558 else if (! std::is_same<double, MT>::value &&
559 plist.isType<double> (paramName)) {
560 const double dblBF = plist.get<double> (paramName);
561 boostFactor = static_cast<MT> (dblBF);
562 }
563 else {
564 TEUCHOS_TEST_FOR_EXCEPTION
565 (true, std::invalid_argument,
566 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
567 "parameter must have type magnitude_type (MT) or double.");
568 }
569 }
570 else { // parameter not in the list
571 // mfh 12 Aug 2016: To preserve current behavior (that fills in
572 // any parameters not in the input ParameterList with their
573 // default values), we call set() here. I don't actually like
574 // this behavior; I prefer the Belos model, where the input
575 // ParameterList is a delta from current behavior. However, I
576 // don't want to break things.
577 plist.set (paramName, defaultBoostFactor);
578 }
579 TEUCHOS_TEST_FOR_EXCEPTION
580 (boostFactor < Teuchos::ScalarTraits<MT>::one (), std::invalid_argument,
581 "Ifpack2::Chebyshev::setParameters: \"" << paramName << "\" parameter "
582 "must be >= 1, but you supplied the value " << boostFactor << ".");
583 }
584
585 // Same name in Ifpack2 and Ifpack.
586 minDiagVal = plist.get ("chebyshev: min diagonal value", minDiagVal);
587 TEUCHOS_TEST_FOR_EXCEPTION(
588 STS::isnaninf (minDiagVal), std::invalid_argument,
589 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
590 "parameter is NaN or Inf. This parameter is optional, but if you choose "
591 "to supply it, it must have a finite value.");
592
593 // Only fill in Ifpack2's name, not ML's or Ifpack's.
594 if (plist.isParameter ("smoother: sweeps")) { // ML compatibility
595 numIters = plist.get<int> ("smoother: sweeps");
596 } // Ifpack's name overrides ML's name.
597 if (plist.isParameter ("relaxation: sweeps")) { // Ifpack compatibility
598 numIters = plist.get<int> ("relaxation: sweeps");
599 } // Ifpack2's name overrides Ifpack's name.
600 numIters = plist.get ("chebyshev: degree", numIters);
601 TEUCHOS_TEST_FOR_EXCEPTION(
602 numIters < 0, std::invalid_argument,
603 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
604 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
605 "nonnegative integer. You gave a value of " << numIters << ".");
606
607 // The last parameter name overrides the first.
608 if (plist.isParameter ("eigen-analysis: iterations")) { // ML compatibility
609 eigMaxIters = plist.get<int> ("eigen-analysis: iterations");
610 } // Ifpack2's name overrides ML's name.
611 eigMaxIters = plist.get ("chebyshev: eigenvalue max iterations", eigMaxIters);
612 TEUCHOS_TEST_FOR_EXCEPTION(
613 eigMaxIters < 0, std::invalid_argument,
614 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
615 "\" parameter (also called \"eigen-analysis: iterations\") must be a "
616 "nonnegative integer. You gave a value of " << eigMaxIters << ".");
617
618 if (plist.isType<double>("chebyshev: eigenvalue relative tolerance"))
619 eigRelTolerance = Teuchos::as<MT>(plist.get<double> ("chebyshev: eigenvalue relative tolerance"));
620 else if (plist.isType<MT>("chebyshev: eigenvalue relative tolerance"))
621 eigRelTolerance = plist.get<MT> ("chebyshev: eigenvalue relative tolerance");
622 else if (plist.isType<ST>("chebyshev: eigenvalue relative tolerance"))
623 eigRelTolerance = Teuchos::ScalarTraits<ST>::magnitude(plist.get<ST> ("chebyshev: eigenvalue relative tolerance"));
624
625 eigKeepVectors = plist.get ("chebyshev: eigenvalue keep vectors", eigKeepVectors);
626
627 eigNormalizationFreq = plist.get ("chebyshev: eigenvalue normalization frequency", eigNormalizationFreq);
628 TEUCHOS_TEST_FOR_EXCEPTION(
629 eigNormalizationFreq < 0, std::invalid_argument,
630 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue normalization frequency"
631 "\" parameter must be a "
632 "nonnegative integer. You gave a value of " << eigNormalizationFreq << ".")
633
634 zeroStartingSolution = plist.get ("chebyshev: zero starting solution",
635 zeroStartingSolution);
636 assumeMatrixUnchanged = plist.get ("chebyshev: assume matrix does not change",
637 assumeMatrixUnchanged);
638
639 // We don't want to fill these parameters in, because they shouldn't
640 // be visible to Ifpack2::Chebyshev users.
641 if (plist.isParameter ("chebyshev: textbook algorithm")) {
642 textbookAlgorithm = plist.get<bool> ("chebyshev: textbook algorithm");
643 }
644 if (plist.isParameter ("chebyshev: compute max residual norm")) {
645 computeMaxResNorm = plist.get<bool> ("chebyshev: compute max residual norm");
646 }
647 if (plist.isParameter ("chebyshev: compute spectral radius")) {
648 computeSpectralRadius = plist.get<bool> ("chebyshev: compute spectral radius");
649 }
650
651 // Test for Ifpack parameters that we won't ever implement here.
652 // Be careful to use the one-argument version of get(), since the
653 // two-argment version adds the parameter if it's not there.
654 TEUCHOS_TEST_FOR_EXCEPTION
655 (plist.isType<bool> ("chebyshev: use block mode") &&
656 ! plist.get<bool> ("chebyshev: use block mode"),
657 std::invalid_argument,
658 "Ifpack2::Chebyshev requires that if you set \"chebyshev: use "
659 "block mode\" at all, you must set it to false. "
660 "Ifpack2::Chebyshev does not implement Ifpack's block mode.");
661 TEUCHOS_TEST_FOR_EXCEPTION
662 (plist.isType<bool> ("chebyshev: solve normal equations") &&
663 ! plist.get<bool> ("chebyshev: solve normal equations"),
664 std::invalid_argument,
665 "Ifpack2::Chebyshev does not and will never implement the Ifpack "
666 "parameter \"chebyshev: solve normal equations\". If you want to "
667 "solve the normal equations, construct a Tpetra::Operator that "
668 "implements A^* A, and use Chebyshev to solve A^* A x = A^* b.");
669
670 // Test for Ifpack parameters that we haven't implemented yet.
671 //
672 // For now, we only check that this ML parameter, if provided, has
673 // the one value that we support. We consider other values "invalid
674 // arguments" rather than "logic errors," because Ifpack does not
675 // implement eigenanalyses other than the power method.
676 std::string eigenAnalysisType ("power-method");
677 if (plist.isParameter ("eigen-analysis: type")) {
678 eigenAnalysisType = plist.get<std::string> ("eigen-analysis: type");
679 TEUCHOS_TEST_FOR_EXCEPTION(
680 eigenAnalysisType != "power-method" &&
681 eigenAnalysisType != "power method" &&
682 eigenAnalysisType != "cg",
683 std::invalid_argument,
684 "Ifpack2::Chebyshev: Ifpack2 only supports \"power method\" and \"cg\" for \"eigen-analysis: type\".");
685 }
686
687 // We've validated all the parameters, so it's safe now to "commit" them.
688 userInvDiag_ = userInvDiagCopy;
689 userLambdaMax_ = lambdaMax;
690 userLambdaMin_ = lambdaMin;
691 userEigRatio_ = eigRatio;
692 boostFactor_ = static_cast<MT> (boostFactor);
693 minDiagVal_ = minDiagVal;
694 numIters_ = numIters;
695 eigMaxIters_ = eigMaxIters;
696 eigRelTolerance_ = eigRelTolerance;
697 eigKeepVectors_ = eigKeepVectors;
698 eigNormalizationFreq_ = eigNormalizationFreq;
699 eigenAnalysisType_ = eigenAnalysisType;
700 zeroStartingSolution_ = zeroStartingSolution;
701 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
702 textbookAlgorithm_ = textbookAlgorithm;
703 computeMaxResNorm_ = computeMaxResNorm;
704 computeSpectralRadius_ = computeSpectralRadius;
705 ckUseNativeSpMV_ = ckUseNativeSpMV;
706 debug_ = debug;
707
708 if (debug_) {
709 // Only print if myRank == 0.
710 int myRank = -1;
711 if (A_.is_null () || A_->getComm ().is_null ()) {
712 // We don't have a communicator (yet), so we assume that
713 // everybody can print. Revise this expectation in setMatrix().
714 myRank = 0;
715 }
716 else {
717 myRank = A_->getComm ()->getRank ();
718 }
719
720 if (myRank == 0) {
721 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
722 }
723 else {
724 using Teuchos::oblackholestream; // prints nothing
725 RCP<oblackholestream> blackHole (new oblackholestream ());
726 out_ = Teuchos::getFancyOStream (blackHole);
727 }
728 }
729 else { // NOT debug
730 // free the "old" output stream, if there was one
731 out_ = Teuchos::null;
732 }
733}
734
735
736template<class ScalarType, class MV>
737void
739{
740 ck_ = Teuchos::null;
741 D_ = Teuchos::null;
742 diagOffsets_ = offsets_type ();
743 savedDiagOffsets_ = false;
744 W_ = Teuchos::null;
745 computedLambdaMax_ = STS::nan ();
746 computedLambdaMin_ = STS::nan ();
747 eigVector_ = Teuchos::null;
748 eigVector2_ = Teuchos::null;
749}
750
751
752template<class ScalarType, class MV>
753void
755setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
756{
757 if (A.getRawPtr () != A_.getRawPtr ()) {
758 if (! assumeMatrixUnchanged_) {
759 reset ();
760 }
761 A_ = A;
762 ck_ = Teuchos::null; // constructed on demand
763
764 // The communicator may have changed, or we may not have had a
765 // communicator before. Thus, we may have to reset the debug
766 // output stream.
767 if (debug_) {
768 // Only print if myRank == 0.
769 int myRank = -1;
770 if (A.is_null () || A->getComm ().is_null ()) {
771 // We don't have a communicator (yet), so we assume that
772 // everybody can print. Revise this expectation in setMatrix().
773 myRank = 0;
774 }
775 else {
776 myRank = A->getComm ()->getRank ();
777 }
778
779 if (myRank == 0) {
780 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
781 }
782 else {
783 Teuchos::RCP<Teuchos::oblackholestream> blackHole (new Teuchos::oblackholestream ());
784 out_ = Teuchos::getFancyOStream (blackHole); // print nothing on other processes
785 }
786 }
787 else { // NOT debug
788 // free the "old" output stream, if there was one
789 out_ = Teuchos::null;
790 }
791 }
792}
793
794template<class ScalarType, class MV>
795void
797{
798 using std::endl;
799 // Some of the optimizations below only work if A_ is a
800 // Tpetra::CrsMatrix. We'll make our best guess about its type
801 // here, since we have no way to get back the original fifth
802 // template parameter.
803 typedef Tpetra::CrsMatrix<typename MV::scalar_type,
804 typename MV::local_ordinal_type,
805 typename MV::global_ordinal_type,
806 typename MV::node_type> crs_matrix_type;
807
808 TEUCHOS_TEST_FOR_EXCEPTION(
809 A_.is_null (), std::runtime_error, "Ifpack2::Chebyshev::compute: The input "
810 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
811 "before calling this method.");
812
813 // If A_ is a CrsMatrix and its graph is constant, we presume that
814 // the user plans to reuse the structure of A_, but possibly change
815 // A_'s values before each compute() call. This is the intended use
816 // case for caching the offsets of the diagonal entries of A_, to
817 // speed up extraction of diagonal entries on subsequent compute()
818 // calls.
819
820 // FIXME (mfh 22 Jan 2013, 10 Feb 2013) In all cases when we use
821 // isnaninf() in this method, we really only want to check if the
822 // number is NaN. Inf means something different. However,
823 // Teuchos::ScalarTraits doesn't distinguish the two cases.
824
825 // makeInverseDiagonal() returns a range Map Vector.
826 if (userInvDiag_.is_null ()) {
827 Teuchos::RCP<const crs_matrix_type> A_crsMat =
828 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
829
830 if (D_.is_null ()) { // We haven't computed D_ before
831 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
832 // It's a CrsMatrix with a const graph; cache diagonal offsets.
833 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
834 if (diagOffsets_.extent (0) < lclNumRows) {
835 diagOffsets_ = offsets_type (); // clear first to save memory
836 diagOffsets_ = offsets_type ("offsets", lclNumRows);
837 }
838 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
839 savedDiagOffsets_ = true;
840 D_ = makeInverseDiagonal (*A_, true);
841 }
842 else { // either A_ is not a CrsMatrix, or its graph is nonconst
843 D_ = makeInverseDiagonal (*A_);
844 }
845 }
846 else if (! assumeMatrixUnchanged_) { // D_ exists but A_ may have changed
847 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
848 // It's a CrsMatrix with a const graph; cache diagonal offsets
849 // if we haven't already.
850 if (! savedDiagOffsets_) {
851 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
852 if (diagOffsets_.extent (0) < lclNumRows) {
853 diagOffsets_ = offsets_type (); // clear first to save memory
854 diagOffsets_ = offsets_type ("offsets", lclNumRows);
855 }
856 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
857 savedDiagOffsets_ = true;
858 }
859 // Now we're guaranteed to have cached diagonal offsets.
860 D_ = makeInverseDiagonal (*A_, true);
861 }
862 else { // either A_ is not a CrsMatrix, or its graph is nonconst
863 D_ = makeInverseDiagonal (*A_);
864 }
865 }
866 }
867 else { // the user provided an inverse diagonal
868 D_ = makeRangeMapVectorConst (userInvDiag_);
869 }
870
871 // Have we estimated eigenvalues before?
872 const bool computedEigenvalueEstimates =
873 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
874
875 // Only recompute the eigenvalue estimates if
876 // - we are supposed to assume that the matrix may have changed, or
877 // - they haven't been computed before, and the user hasn't given
878 // us at least an estimate of the max eigenvalue.
879 //
880 // We at least need an estimate of the max eigenvalue. This is the
881 // most important one if using Chebyshev as a smoother.
882
883 if (! assumeMatrixUnchanged_ ||
884 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
885 ST computedLambdaMax;
886 if ((eigenAnalysisType_ == "power method") || (eigenAnalysisType_ == "power-method")) {
887 Teuchos::RCP<V> x;
888 if (eigVector_.is_null()) {
889 x = Teuchos::rcp(new V(A_->getDomainMap ()));
890 if (eigKeepVectors_)
891 eigVector_ = x;
892 PowerMethod::computeInitialGuessForPowerMethod (*x, false);
893 } else
894 x = eigVector_;
895
896 Teuchos::RCP<V> y;
897 if (eigVector2_.is_null()) {
898 y = rcp(new V(A_->getRangeMap ()));
899 if (eigKeepVectors_)
900 eigVector2_ = y;
901 } else
902 y = eigVector2_;
903
904 Teuchos::RCP<Teuchos::FancyOStream> stream = (debug_ ? out_ : Teuchos::null);
905 computedLambdaMax = PowerMethod::powerMethodWithInitGuess (*A_, *D_, eigMaxIters_, x, y,
906 eigRelTolerance_, eigNormalizationFreq_, stream,
907 computeSpectralRadius_);
908 }
909 else
910 computedLambdaMax = cgMethod (*A_, *D_, eigMaxIters_);
911 TEUCHOS_TEST_FOR_EXCEPTION(
912 STS::isnaninf (computedLambdaMax),
913 std::runtime_error,
914 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
915 "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
916 "the matrix contains Inf or NaN values, or that it is badly scaled.");
917 TEUCHOS_TEST_FOR_EXCEPTION(
918 STS::isnaninf (userEigRatio_),
919 std::logic_error,
920 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
921 << endl << "This should be impossible." << endl <<
922 "Please report this bug to the Ifpack2 developers.");
923
924 // The power method doesn't estimate the min eigenvalue, so we
925 // do our best to provide an estimate. userEigRatio_ has a
926 // reasonable default value, and if the user provided it, we
927 // have already checked that its value is finite and >= 1.
928 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
929
930 // Defer "committing" results until all computations succeeded.
931 computedLambdaMax_ = computedLambdaMax;
932 computedLambdaMin_ = computedLambdaMin;
933 } else {
934 TEUCHOS_TEST_FOR_EXCEPTION(
935 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
936 std::logic_error,
937 "Ifpack2::Chebyshev::compute: " << endl <<
938 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
939 << endl << "This should be impossible." << endl <<
940 "Please report this bug to the Ifpack2 developers.");
941 }
942
944 // Figure out the eigenvalue estimates that apply() will use.
946
947 // Always favor the user's max eigenvalue estimate, if provided.
948 lambdaMaxForApply_ = STS::isnaninf (userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
949
950 // mfh 11 Feb 2013: For now, we imitate Ifpack by ignoring the
951 // user's min eigenvalue estimate, and using the given eigenvalue
952 // ratio to estimate the min eigenvalue. We could instead do this:
953 // favor the user's eigenvalue ratio estimate, but if it's not
954 // provided, use lambdaMax / lambdaMin. However, we want Chebyshev
955 // to have sensible smoother behavior if the user did not provide
956 // eigenvalue estimates. Ifpack's behavior attempts to push down
957 // the error terms associated with the largest eigenvalues, while
958 // expecting that users will only want a small number of iterations,
959 // so that error terms associated with the smallest eigenvalues
960 // won't grow too much. This is sensible behavior for a smoother.
961 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
962 eigRatioForApply_ = userEigRatio_;
963
964 if (! textbookAlgorithm_) {
965 // Ifpack has a special-case modification of the eigenvalue bounds
966 // for the case where the max eigenvalue estimate is close to one.
967 const ST one = Teuchos::as<ST> (1);
968 // FIXME (mfh 20 Nov 2013) Should scale this 1.0e-6 term
969 // appropriately for MT's machine precision.
970 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
971 lambdaMinForApply_ = one;
972 lambdaMaxForApply_ = lambdaMinForApply_;
973 eigRatioForApply_ = one; // Ifpack doesn't include this line.
974 }
975 }
976}
977
978
979template<class ScalarType, class MV>
980ScalarType
982getLambdaMaxForApply() const {
983 return lambdaMaxForApply_;
984}
985
986
987template<class ScalarType, class MV>
990{
991 const char prefix[] = "Ifpack2::Chebyshev::apply: ";
992
993 if (debug_) {
994 *out_ << "apply: " << std::endl;
995 }
996 TEUCHOS_TEST_FOR_EXCEPTION
997 (A_.is_null (), std::runtime_error, prefix << "The input matrix A is null. "
998 " Please call setMatrix() with a nonnull input matrix, and then call "
999 "compute(), before calling this method.");
1000 TEUCHOS_TEST_FOR_EXCEPTION
1001 (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
1002 prefix << "There is no estimate for the max eigenvalue."
1003 << std::endl << computeBeforeApplyReminder);
1004 TEUCHOS_TEST_FOR_EXCEPTION
1005 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
1006 prefix << "There is no estimate for the min eigenvalue."
1007 << std::endl << computeBeforeApplyReminder);
1008 TEUCHOS_TEST_FOR_EXCEPTION
1009 (STS::isnaninf (eigRatioForApply_), std::runtime_error,
1010 prefix <<"There is no estimate for the ratio of the max "
1011 "eigenvalue to the min eigenvalue."
1012 << std::endl << computeBeforeApplyReminder);
1013 TEUCHOS_TEST_FOR_EXCEPTION
1014 (D_.is_null (), std::runtime_error, prefix << "The vector of inverse "
1015 "diagonal entries of the matrix has not yet been computed."
1016 << std::endl << computeBeforeApplyReminder);
1017
1018 if (textbookAlgorithm_) {
1019 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1020 lambdaMinForApply_, eigRatioForApply_, *D_);
1021 }
1022 else {
1023 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1024 lambdaMinForApply_, eigRatioForApply_, *D_);
1025 }
1026
1027 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
1028 MV R (B.getMap (), B.getNumVectors ());
1029 computeResidual (R, B, *A_, X);
1030 Teuchos::Array<MT> norms (B.getNumVectors ());
1031 R.norm2 (norms ());
1032 return *std::max_element (norms.begin (), norms.end ());
1033 }
1034 else {
1035 return Teuchos::ScalarTraits<MT>::zero ();
1036 }
1037}
1038
1039template<class ScalarType, class MV>
1040void
1042print (std::ostream& out)
1043{
1044 using Teuchos::rcpFromRef;
1045 this->describe (* (Teuchos::getFancyOStream (rcpFromRef (out))),
1046 Teuchos::VERB_MEDIUM);
1047}
1048
1049template<class ScalarType, class MV>
1050void
1053 const ScalarType& alpha,
1054 const V& D_inv,
1055 const MV& B,
1056 MV& X)
1057{
1058 solve (W, alpha, D_inv, B); // W = alpha*D_inv*B
1059 Tpetra::deep_copy (X, W); // X = 0 + W
1060}
1061
1062template<class ScalarType, class MV>
1063void
1065computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
1066 const Teuchos::ETransp mode)
1067{
1068 Tpetra::Details::residual(A,X,B,R);
1069}
1070
1071template <class ScalarType, class MV>
1072void
1073Chebyshev<ScalarType, MV>::
1074solve (MV& Z, const V& D_inv, const MV& R) {
1075 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
1076}
1077
1078template<class ScalarType, class MV>
1079void
1080Chebyshev<ScalarType, MV>::
1081solve (MV& Z, const ST alpha, const V& D_inv, const MV& R) {
1082 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1083}
1084
1085template<class ScalarType, class MV>
1086Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1087Chebyshev<ScalarType, MV>::
1088makeInverseDiagonal (const row_matrix_type& A, const bool useDiagOffsets) const
1089{
1090 using Teuchos::RCP;
1091 using Teuchos::rcpFromRef;
1092 using Teuchos::rcp_dynamic_cast;
1093
1094 RCP<V> D_rowMap;
1095 if (!D_.is_null() &&
1096 D_->getMap()->isSameAs(*(A.getRowMap ()))) {
1097 if (debug_)
1098 *out_ << "Reusing pre-existing vector for diagonal extraction" << std::endl;
1099 D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1100 } else {
1101 D_rowMap = Teuchos::rcp(new V (A.getRowMap (), /*zeroOut=*/false));
1102 if (debug_)
1103 *out_ << "Allocated new vector for diagonal extraction" << std::endl;
1104 }
1105 if (useDiagOffsets) {
1106 // The optimizations below only work if A_ is a Tpetra::CrsMatrix.
1107 // We'll make our best guess about its type here, since we have no
1108 // way to get back the original fifth template parameter.
1109 typedef Tpetra::CrsMatrix<typename MV::scalar_type,
1110 typename MV::local_ordinal_type,
1111 typename MV::global_ordinal_type,
1112 typename MV::node_type> crs_matrix_type;
1113 RCP<const crs_matrix_type> A_crsMat =
1114 rcp_dynamic_cast<const crs_matrix_type> (rcpFromRef (A));
1115 if (! A_crsMat.is_null ()) {
1116 TEUCHOS_TEST_FOR_EXCEPTION(
1117 ! savedDiagOffsets_, std::logic_error,
1118 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
1119 "It is not allowed to call this method with useDiagOffsets=true, "
1120 "if you have not previously saved offsets of diagonal entries. "
1121 "This situation should never arise if this class is used properly. "
1122 "Please report this bug to the Ifpack2 developers.");
1123 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1124 }
1125 }
1126 else {
1127 // This always works for a Tpetra::RowMatrix, even if it is not a
1128 // Tpetra::CrsMatrix. We just don't have offsets in this case.
1129 A.getLocalDiagCopy (*D_rowMap);
1130 }
1131 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1132
1133 if (debug_) {
1134 // In debug mode, make sure that all diagonal entries are
1135 // positive, on all processes. Note that *out_ only prints on
1136 // Process 0 of the matrix's communicator.
1137 bool foundNonpositiveValue = false;
1138 {
1139 auto D_lcl = D_rangeMap->getLocalViewHost (Tpetra::Access::ReadOnly);
1140 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1141
1142 typedef typename MV::impl_scalar_type IST;
1143 typedef typename MV::local_ordinal_type LO;
1144 typedef Kokkos::Details::ArithTraits<IST> ATS;
1145 typedef Kokkos::Details::ArithTraits<typename ATS::mag_type> STM;
1146
1147 const LO lclNumRows = static_cast<LO> (D_rangeMap->getLocalLength ());
1148 for (LO i = 0; i < lclNumRows; ++i) {
1149 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1150 foundNonpositiveValue = true;
1151 break;
1152 }
1153 }
1154 }
1155
1156 using Teuchos::outArg;
1157 using Teuchos::REDUCE_MIN;
1158 using Teuchos::reduceAll;
1159
1160 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1161 int gblSuccess = lclSuccess; // to be overwritten
1162 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1163 const Teuchos::Comm<int>& comm = * (D_rangeMap->getMap ()->getComm ());
1164 reduceAll<int, int> (comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1165 }
1166 if (gblSuccess == 1) {
1167 *out_ << "makeInverseDiagonal: The matrix's diagonal entries all have "
1168 "positive real part (this is good for Chebyshev)." << std::endl;
1169 }
1170 else {
1171 *out_ << "makeInverseDiagonal: The matrix's diagonal has at least one "
1172 "entry with nonpositive real part, on at least one process in the "
1173 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1174 }
1175 }
1176
1177 // Invert the diagonal entries, replacing entries less (in
1178 // magnitude) than the user-specified value with that value.
1179 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1180 return Teuchos::rcp_const_cast<const V> (D_rangeMap);
1181}
1182
1183
1184template<class ScalarType, class MV>
1185Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1186Chebyshev<ScalarType, MV>::
1187makeRangeMapVectorConst (const Teuchos::RCP<const V>& D) const
1188{
1189 using Teuchos::RCP;
1190 using Teuchos::rcp;
1191 typedef Tpetra::Export<typename MV::local_ordinal_type,
1192 typename MV::global_ordinal_type,
1193 typename MV::node_type> export_type;
1194 // This throws logic_error instead of runtime_error, because the
1195 // methods that call makeRangeMapVector should all have checked
1196 // whether A_ is null before calling this method.
1197 TEUCHOS_TEST_FOR_EXCEPTION(
1198 A_.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
1199 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1200 "with a nonnull input matrix before calling this method. This is probably "
1201 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1202 TEUCHOS_TEST_FOR_EXCEPTION(
1203 D.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
1204 "makeRangeMapVector: The input Vector D is null. This is probably "
1205 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1206
1207 RCP<const map_type> sourceMap = D->getMap ();
1208 RCP<const map_type> rangeMap = A_->getRangeMap ();
1209 RCP<const map_type> rowMap = A_->getRowMap ();
1210
1211 if (rangeMap->isSameAs (*sourceMap)) {
1212 // The given vector's Map is the same as the matrix's range Map.
1213 // That means we don't need to Export. This is the fast path.
1214 return D;
1215 }
1216 else { // We need to Export.
1217 RCP<const export_type> exporter;
1218 // Making an Export object from scratch is expensive enough that
1219 // it's worth the O(1) global reductions to call isSameAs(), to
1220 // see if we can avoid that cost.
1221 if (sourceMap->isSameAs (*rowMap)) {
1222 // We can reuse the matrix's Export object, if there is one.
1223 exporter = A_->getGraph ()->getExporter ();
1224 }
1225 else { // We have to make a new Export object.
1226 exporter = rcp (new export_type (sourceMap, rangeMap));
1227 }
1228
1229 if (exporter.is_null ()) {
1230 return D; // Row Map and range Map are the same; no need to Export.
1231 }
1232 else { // Row Map and range Map are _not_ the same; must Export.
1233 RCP<V> D_out = rcp (new V (*D, Teuchos::Copy));
1234 D_out->doExport (*D, *exporter, Tpetra::ADD);
1235 return Teuchos::rcp_const_cast<const V> (D_out);
1236 }
1237 }
1238}
1239
1240
1241template<class ScalarType, class MV>
1242Teuchos::RCP<typename Chebyshev<ScalarType, MV>::V>
1243Chebyshev<ScalarType, MV>::
1244makeRangeMapVector (const Teuchos::RCP<V>& D) const
1245{
1246 using Teuchos::rcp_const_cast;
1247 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1248}
1249
1250
1251template<class ScalarType, class MV>
1252void
1253Chebyshev<ScalarType, MV>::
1254textbookApplyImpl (const op_type& A,
1255 const MV& B,
1256 MV& X,
1257 const int numIters,
1258 const ST lambdaMax,
1259 const ST lambdaMin,
1260 const ST eigRatio,
1261 const V& D_inv) const
1262{
1263 (void) lambdaMin; // Forestall compiler warning.
1264 const ST myLambdaMin = lambdaMax / eigRatio;
1265
1266 const ST zero = Teuchos::as<ST> (0);
1267 const ST one = Teuchos::as<ST> (1);
1268 const ST two = Teuchos::as<ST> (2);
1269 const ST d = (lambdaMax + myLambdaMin) / two; // Ifpack2 calls this theta
1270 const ST c = (lambdaMax - myLambdaMin) / two; // Ifpack2 calls this 1/delta
1271
1272 if (zeroStartingSolution_ && numIters > 0) {
1273 // If zero iterations, then input X is output X.
1274 X.putScalar (zero);
1275 }
1276 MV R (B.getMap (), B.getNumVectors (), false);
1277 MV P (B.getMap (), B.getNumVectors (), false);
1278 MV Z (B.getMap (), B.getNumVectors (), false);
1279 ST alpha, beta;
1280 for (int i = 0; i < numIters; ++i) {
1281 computeResidual (R, B, A, X); // R = B - A*X
1282 solve (Z, D_inv, R); // z = D_inv * R, that is, D \ R.
1283 if (i == 0) {
1284 P = Z;
1285 alpha = two / d;
1286 } else {
1287 //beta = (c * alpha / two)^2;
1288 //const ST sqrtBeta = c * alpha / two;
1289 //beta = sqrtBeta * sqrtBeta;
1290 beta = alpha * (c/two) * (c/two);
1291 alpha = one / (d - beta);
1292 P.update (one, Z, beta); // P = Z + beta*P
1293 }
1294 X.update (alpha, P, one); // X = X + alpha*P
1295 // If we compute the residual here, we could either do R = B -
1296 // A*X, or R = R - alpha*A*P. Since we choose the former, we
1297 // can move the computeResidual call to the top of the loop.
1298 }
1299}
1300
1301template<class ScalarType, class MV>
1303Chebyshev<ScalarType, MV>::maxNormInf (const MV& X) {
1304 Teuchos::Array<MT> norms (X.getNumVectors ());
1305 X.normInf (norms());
1306 return *std::max_element (norms.begin (), norms.end ());
1307}
1308
1309template<class ScalarType, class MV>
1310void
1311Chebyshev<ScalarType, MV>::
1312ifpackApplyImpl (const op_type& A,
1313 const MV& B,
1314 MV& X,
1315 const int numIters,
1316 const ST lambdaMax,
1317 const ST lambdaMin,
1318 const ST eigRatio,
1319 const V& D_inv)
1320{
1321 using std::endl;
1322#ifdef HAVE_IFPACK2_DEBUG
1323 const bool debug = debug_;
1324#else
1325 const bool debug = false;
1326#endif
1327
1328 if (debug) {
1329 *out_ << " \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1330 *out_ << " \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1331 }
1332
1333 if (numIters <= 0) {
1334 return;
1335 }
1336 const ST zero = static_cast<ST> (0.0);
1337 const ST one = static_cast<ST> (1.0);
1338 const ST two = static_cast<ST> (2.0);
1339
1340 // Quick solve when the matrix A is the identity.
1341 if (lambdaMin == one && lambdaMax == lambdaMin) {
1342 solve (X, D_inv, B);
1343 return;
1344 }
1345
1346 // Initialize coefficients
1347 const ST alpha = lambdaMax / eigRatio;
1348 const ST beta = boostFactor_ * lambdaMax;
1349 const ST delta = two / (beta - alpha);
1350 const ST theta = (beta + alpha) / two;
1351 const ST s1 = theta * delta;
1352
1353 if (debug) {
1354 *out_ << " alpha = " << alpha << endl
1355 << " beta = " << beta << endl
1356 << " delta = " << delta << endl
1357 << " theta = " << theta << endl
1358 << " s1 = " << s1 << endl;
1359 }
1360
1361 // Fetch cached temporary (multi)vector.
1362 Teuchos::RCP<MV> W_ptr = makeTempMultiVector (B);
1363 MV& W = *W_ptr;
1364
1365 if (debug) {
1366 *out_ << " Iteration " << 1 << ":" << endl
1367 << " - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1368 }
1369
1370 // Special case for the first iteration.
1371 if (! zeroStartingSolution_) {
1372 // mfh 22 May 2019: Tests don't actually exercise this path.
1373
1374 if (ck_.is_null ()) {
1375 Teuchos::RCP<const op_type> A_op = A_;
1376 ck_ = Teuchos::rcp (new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1377 }
1378 // W := (1/theta)*D_inv*(B-A*X) and X := X + W.
1379 // X := X + W
1380 ck_->compute (W, one/theta, const_cast<V&> (D_inv),
1381 const_cast<MV&> (B), X, zero);
1382 }
1383 else {
1384 // W := (1/theta)*D_inv*B and X := 0 + W.
1385 firstIterationWithZeroStartingSolution (W, one/theta, D_inv, B, X);
1386 }
1387
1388 if (debug) {
1389 *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1390 << " - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1391 }
1392
1393 if (numIters > 1 && ck_.is_null ()) {
1394 Teuchos::RCP<const op_type> A_op = A_;
1395 ck_ = Teuchos::rcp (new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1396 }
1397
1398 // The rest of the iterations.
1399 ST rhok = one / s1;
1400 ST rhokp1, dtemp1, dtemp2;
1401 for (int deg = 1; deg < numIters; ++deg) {
1402 if (debug) {
1403 *out_ << " Iteration " << deg+1 << ":" << endl
1404 << " - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1405 << " - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1406 << " - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm ()
1407 << endl << " - rhok = " << rhok << endl;
1408 }
1409
1410 rhokp1 = one / (two * s1 - rhok);
1411 dtemp1 = rhokp1 * rhok;
1412 dtemp2 = two * rhokp1 * delta;
1413 rhok = rhokp1;
1414
1415 if (debug) {
1416 *out_ << " - dtemp1 = " << dtemp1 << endl
1417 << " - dtemp2 = " << dtemp2 << endl;
1418 }
1419
1420 // W := dtemp2*D_inv*(B - A*X) + dtemp1*W.
1421 // X := X + W
1422 ck_->compute (W, dtemp2, const_cast<V&> (D_inv),
1423 const_cast<MV&> (B), (X), dtemp1);
1424
1425 if (debug) {
1426 *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1427 << " - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1428 }
1429 }
1430}
1431
1432
1433template<class ScalarType, class MV>
1435Chebyshev<ScalarType, MV>::
1436cgMethodWithInitGuess (const op_type& A,
1437 const V& D_inv,
1438 const int numIters,
1439 V& r)
1440{
1441 using std::endl;
1442 using MagnitudeType = typename STS::magnitudeType;
1443 if (debug_) {
1444 *out_ << " cgMethodWithInitGuess:" << endl;
1445 }
1446
1447 const ST one = STS::one();
1448 ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1449 // ST lambdaMin;
1450 Teuchos::ArrayRCP<MagnitudeType> diag, offdiag;
1451 Teuchos::RCP<V> p, z, Ap;
1452 diag.resize(numIters);
1453 offdiag.resize(numIters-1);
1454
1455 p = rcp(new V(A.getRangeMap ()));
1456 z = rcp(new V(A.getRangeMap ()));
1457 Ap = rcp(new V(A.getRangeMap ()));
1458
1459 // Tpetra::Details::residual (A, x, *b, *r);
1460 solve (*p, D_inv, r);
1461 rHz = r.dot (*p);
1462
1463 for (int iter = 0; iter < numIters; ++iter) {
1464 if (debug_) {
1465 *out_ << " Iteration " << iter << endl;
1466 }
1467 A.apply (*p, *Ap);
1468 pAp = p->dot (*Ap);
1469 alpha = rHz/pAp;
1470 r.update (-alpha, *Ap, one);
1471 rHzOld = rHz;
1472 solve (*z, D_inv, r);
1473 rHz = r.dot (*z);
1474 beta = rHz / rHzOld;
1475 p->update (one, *z, beta);
1476 if (iter>0) {
1477 diag[iter] = STS::real((betaOld*betaOld * pApOld + pAp) / rHzOld);
1478 offdiag[iter-1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1479 if (debug_) {
1480 *out_ << " diag[" << iter << "] = " << diag[iter] << endl;
1481 *out_ << " offdiag["<< iter-1 << "] = " << offdiag[iter-1] << endl;
1482 }
1483 } else {
1484 diag[iter] = STS::real(pAp/rHzOld);
1485 if (debug_) {
1486 *out_ << " diag[" << iter << "] = " << diag[iter] << endl;
1487 }
1488 }
1489 rHzOld2 = rHzOld;
1490 betaOld = beta;
1491 pApOld = pAp;
1492 }
1493
1494 lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1495
1496 return lambdaMax;
1497}
1498
1499
1500template<class ScalarType, class MV>
1502Chebyshev<ScalarType, MV>::
1503cgMethod (const op_type& A, const V& D_inv, const int numIters)
1504{
1505 using std::endl;
1506 if (debug_) {
1507 *out_ << "cgMethod:" << endl;
1508 }
1509
1510 Teuchos::RCP<V> r;
1511 if (eigVector_.is_null()) {
1512 r = rcp(new V(A.getDomainMap ()));
1513 if (eigKeepVectors_)
1514 eigVector_ = r;
1515 // For the first pass, just let the pseudorandom number generator
1516 // fill x with whatever values it wants; don't try to make its
1517 // entries nonnegative.
1518 PowerMethod::computeInitialGuessForPowerMethod (*r, false);
1519 } else
1520 r = eigVector_;
1521
1522 ST lambdaMax = cgMethodWithInitGuess (A, D_inv, numIters, *r);
1523
1524 return lambdaMax;
1525}
1526
1527template<class ScalarType, class MV>
1528Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
1530 return A_;
1531}
1532
1533template<class ScalarType, class MV>
1534bool
1536hasTransposeApply () const {
1537 // Technically, this is true, because the matrix must be symmetric.
1538 return true;
1539}
1540
1541template<class ScalarType, class MV>
1542Teuchos::RCP<MV>
1544makeTempMultiVector (const MV& B)
1545{
1546 // ETP 02/08/17: We must check not only if the temporary vectors are
1547 // null, but also if the number of columns match, since some multi-RHS
1548 // solvers (e.g., Belos) may call apply() with different numbers of columns.
1549
1550 const size_t B_numVecs = B.getNumVectors ();
1551 if (W_.is_null () || W_->getNumVectors () != B_numVecs) {
1552 W_ = Teuchos::rcp (new MV (B.getMap (), B_numVecs, false));
1553 }
1554 return W_;
1555}
1556
1557template<class ScalarType, class MV>
1558std::string
1560description () const {
1561 std::ostringstream oss;
1562 // YAML requires quoting the key in this case, to distinguish
1563 // key's colons from the colon that separates key from value.
1564 oss << "\"Ifpack2::Details::Chebyshev\":"
1565 << "{"
1566 << "degree: " << numIters_
1567 << ", lambdaMax: " << lambdaMaxForApply_
1568 << ", alpha: " << eigRatioForApply_
1569 << ", lambdaMin: " << lambdaMinForApply_
1570 << ", boost factor: " << boostFactor_;
1571 if (!userInvDiag_.is_null())
1572 oss << ", diagonal: user-supplied";
1573 oss << "}";
1574 return oss.str();
1575}
1576
1577template<class ScalarType, class MV>
1578void
1580describe (Teuchos::FancyOStream& out,
1581 const Teuchos::EVerbosityLevel verbLevel) const
1582{
1583 using Teuchos::TypeNameTraits;
1584 using std::endl;
1585
1586 const Teuchos::EVerbosityLevel vl =
1587 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1588 if (vl == Teuchos::VERB_NONE) {
1589 return; // print NOTHING
1590 }
1591
1592 // By convention, describe() starts with a tab.
1593 //
1594 // This does affect all processes on which it's valid to print to
1595 // 'out'. However, it does not actually print spaces to 'out'
1596 // unless operator<< gets called, so it's safe to use on all
1597 // processes.
1598 Teuchos::OSTab tab0 (out);
1599
1600 // We only print on Process 0 of the matrix's communicator. If
1601 // the matrix isn't set, we don't have a communicator, so we have
1602 // to assume that every process can print.
1603 int myRank = -1;
1604 if (A_.is_null () || A_->getComm ().is_null ()) {
1605 myRank = 0;
1606 }
1607 else {
1608 myRank = A_->getComm ()->getRank ();
1609 }
1610 if (myRank == 0) {
1611 // YAML requires quoting the key in this case, to distinguish
1612 // key's colons from the colon that separates key from value.
1613 out << "\"Ifpack2::Details::Chebyshev\":" << endl;
1614 }
1615 Teuchos::OSTab tab1 (out);
1616
1617 if (vl == Teuchos::VERB_LOW) {
1618 if (myRank == 0) {
1619 out << "degree: " << numIters_ << endl
1620 << "lambdaMax: " << lambdaMaxForApply_ << endl
1621 << "alpha: " << eigRatioForApply_ << endl
1622 << "lambdaMin: " << lambdaMinForApply_ << endl
1623 << "boost factor: " << boostFactor_ << endl;
1624 }
1625 return;
1626 }
1627
1628 // vl > Teuchos::VERB_LOW
1629
1630 if (myRank == 0) {
1631 out << "Template parameters:" << endl;
1632 {
1633 Teuchos::OSTab tab2 (out);
1634 out << "ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1635 << "MV: " << TypeNameTraits<MV>::name () << endl;
1636 }
1637
1638 // "Computed parameters" literally means "parameters whose
1639 // values were computed by compute()."
1640 if (myRank == 0) {
1641 out << "Computed parameters:" << endl;
1642 }
1643 }
1644
1645 // Print computed parameters
1646 {
1647 Teuchos::OSTab tab2 (out);
1648 // Users might want to see the values in the computed inverse
1649 // diagonal, so we print them out at the highest verbosity.
1650 if (myRank == 0) {
1651 out << "D_: ";
1652 }
1653 if (D_.is_null ()) {
1654 if (myRank == 0) {
1655 out << "unset" << endl;
1656 }
1657 }
1658 else if (vl <= Teuchos::VERB_HIGH) {
1659 if (myRank == 0) {
1660 out << "set" << endl;
1661 }
1662 }
1663 else { // D_ not null and vl > Teuchos::VERB_HIGH
1664 if (myRank == 0) {
1665 out << endl;
1666 }
1667 // By convention, describe() first indents, then prints.
1668 // We can rely on other describe() implementations to do that.
1669 D_->describe (out, vl);
1670 }
1671 if (myRank == 0) {
1672 // W_ is scratch space; its values are irrelevant.
1673 // All that matters is whether or not they have been set.
1674 out << "W_: " << (W_.is_null () ? "unset" : "set") << endl
1675 << "computedLambdaMax_: " << computedLambdaMax_ << endl
1676 << "computedLambdaMin_: " << computedLambdaMin_ << endl
1677 << "lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1678 << "lambdaMinForApply_: " << lambdaMinForApply_ << endl
1679 << "eigRatioForApply_: " << eigRatioForApply_ << endl;
1680 }
1681 } // print computed parameters
1682
1683 if (myRank == 0) {
1684 out << "User parameters:" << endl;
1685 }
1686
1687 // Print user parameters
1688 {
1689 Teuchos::OSTab tab2 (out);
1690 out << "userInvDiag_: ";
1691 if (userInvDiag_.is_null ()) {
1692 out << "unset" << endl;
1693 }
1694 else if (vl <= Teuchos::VERB_HIGH) {
1695 out << "set" << endl;
1696 }
1697 else { // userInvDiag_ not null and vl > Teuchos::VERB_HIGH
1698 if (myRank == 0) {
1699 out << endl;
1700 }
1701 userInvDiag_->describe (out, vl);
1702 }
1703 if (myRank == 0) {
1704 out << "userLambdaMax_: " << userLambdaMax_ << endl
1705 << "userLambdaMin_: " << userLambdaMin_ << endl
1706 << "userEigRatio_: " << userEigRatio_ << endl
1707 << "numIters_: " << numIters_ << endl
1708 << "eigMaxIters_: " << eigMaxIters_ << endl
1709 << "eigRelTolerance_: " << eigRelTolerance_ << endl
1710 << "eigNormalizationFreq_: " << eigNormalizationFreq_ << endl
1711 << "zeroStartingSolution_: " << zeroStartingSolution_ << endl
1712 << "assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1713 << "textbookAlgorithm_: " << textbookAlgorithm_ << endl
1714 << "computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1715 }
1716 } // print user parameters
1717}
1718
1719} // namespace Details
1720} // namespace Ifpack2
1721
1722#define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1723 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
1724
1725#endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
Declaration of Chebyshev implementation.
Definition of power methods.
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:208
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< scalar_type > 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
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:287
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
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