Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_PowerMethod.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4 //
5 // Copyright 2009 NTESS and the Ifpack2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef IFPACK2_POWERMETHOD_HPP
11 #define IFPACK2_POWERMETHOD_HPP
12 
19 
20 #if KOKKOS_VERSION >= 40799
21 #include "KokkosKernels_ArithTraits.hpp"
22 #else
23 #include "Kokkos_ArithTraits.hpp"
24 #endif
25 #include "Teuchos_FancyOStream.hpp"
26 #include "Teuchos_oblackholestream.hpp"
27 #include "Tpetra_Details_residual.hpp"
28 #include <cmath>
29 #include <iostream>
30 
31 namespace Ifpack2 {
32 namespace PowerMethod {
33 
34 namespace { // (anonymous)
35 
36 // Functor for making sure the real parts of all entries of a vector
37 // are nonnegative. We use this in computeInitialGuessForPowerMethod
38 // below.
39 template <class OneDViewType,
40  class LocalOrdinal = typename OneDViewType::size_type>
41 class PositivizeVector {
42  static_assert(Kokkos::is_view<OneDViewType>::value,
43  "OneDViewType must be a 1-D Kokkos::View.");
44  static_assert(static_cast<int>(OneDViewType::rank) == 1,
45  "This functor only works with a 1-D View.");
46  static_assert(std::is_integral<LocalOrdinal>::value,
47  "The second template parameter, LocalOrdinal, "
48  "must be an integer type.");
49 
50  public:
51  PositivizeVector(const OneDViewType& x)
52  : x_(x) {}
53 
54  KOKKOS_INLINE_FUNCTION void
55  operator()(const LocalOrdinal& i) const {
56  typedef typename OneDViewType::non_const_value_type IST;
57 #if KOKKOS_VERSION >= 40799
58  typedef KokkosKernels::ArithTraits<IST> STS;
59 #else
60  typedef Kokkos::ArithTraits<IST> STS;
61 #endif
62 #if KOKKOS_VERSION >= 40799
63  typedef KokkosKernels::ArithTraits<typename STS::mag_type> STM;
64 #else
65  typedef Kokkos::ArithTraits<typename STS::mag_type> STM;
66 #endif
67 
68  if (STS::real(x_(i)) < STM::zero()) {
69  x_(i) = -x_(i);
70  }
71  }
72 
73  private:
74  OneDViewType x_;
75 };
76 
77 } // namespace
78 
98 template <class OperatorType, class V>
99 typename V::scalar_type
100 powerMethodWithInitGuess(const OperatorType& A,
101  const V& D_inv,
102  const int numIters,
103  Teuchos::RCP<V> x,
104  Teuchos::RCP<V> y,
106  const int eigNormalizationFreq = 1,
107  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::null,
108  const bool computeSpectralRadius = true) {
109  typedef typename V::scalar_type ST;
112 
113  bool verbose = (out != Teuchos::null);
114 
115  using std::endl;
116  if (verbose) {
117  *out << " powerMethodWithInitGuess:" << endl;
118  }
119 
120  const ST zero = static_cast<ST>(0.0);
121  const ST one = static_cast<ST>(1.0);
122  ST lambdaMax = zero;
123  ST lambdaMaxOld = one;
124  ST norm;
125 
126  norm = x->norm2();
127  TEUCHOS_TEST_FOR_EXCEPTION(norm == zero, std::runtime_error,
128  "Ifpack2::PowerMethod::powerMethodWithInitGuess: The initial guess "
129  "has zero norm. This could be either because Tpetra::Vector::"
130  "randomize filled the vector with zeros (if that was used to "
131  "compute the initial guess), or because the norm2 method has a "
132  "bug. The first is not impossible, but unlikely.");
133 
134  if (verbose) {
135  *out << " Original norm1(x): " << x->norm1()
136  << ", norm2(x): " << norm << endl;
137  }
138 
139  x->scale(one / norm);
140 
141  if (verbose) {
142  *out << " norm1(x.scale(one/norm)): " << x->norm1() << endl;
143  }
144 
145  if (y.is_null())
146  y = Teuchos::rcp(new V(A.getRangeMap()));
147 
148  // GH 04 Nov 2022: The previous implementation of the power method was inconsistent with itself.
149  // It computed x->norm2() inside the loop, but returned x^T*Ax.
150  // This returned horribly incorrect estimates for Chebyshev spectral radius in certain
151  // cases, such as the Cheby_belos test, which has two dominant eigenvalues of 12.5839, -10.6639.
152  // In such case, the power method iteration history would appear to converge to something
153  // between 10 and 12, but the resulting spectral radius estimate was sometimes negative.
154  // These have now been split into two routes which behave consistently with themselves.
155  if (computeSpectralRadius) {
156  // In this route, the update is as follows:
157  // y=A*x, x = Dinv*y, lambda = norm(x), x = x/lambda
158  if (verbose) {
159  *out << " PowerMethod using spectral radius route" << endl;
160  }
161  for (int iter = 0; iter < numIters - 1; ++iter) {
162  if (verbose) {
163  *out << " Iteration " << iter << endl;
164  }
165  A.apply(*x, *y);
166  x->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
167 
168  if (((iter + 1) % eigNormalizationFreq == 0) && (iter < numIters - 2)) {
169  norm = x->norm2();
170  if (norm == zero) { // Return something reasonable.
171  if (verbose) {
172  *out << " norm is zero; returning zero" << endl;
173  *out << " Power method terminated after " << iter << " iterations." << endl;
174  }
175  return zero;
176  } else {
177  lambdaMaxOld = lambdaMax;
178  lambdaMax = pow(norm, Teuchos::ScalarTraits<MT>::one() / eigNormalizationFreq);
179  if (Teuchos::ScalarTraits<ST>::magnitude(lambdaMax - lambdaMaxOld) < tolerance * Teuchos::ScalarTraits<ST>::magnitude(lambdaMax)) {
180  if (verbose) {
181  *out << " lambdaMax: " << lambdaMax << endl;
182  *out << " Power method terminated after " << iter << " iterations." << endl;
183  }
184  return lambdaMax;
185  } else if (verbose) {
186  *out << " lambdaMaxOld: " << lambdaMaxOld << endl;
187  *out << " lambdaMax: " << lambdaMax << endl;
188  *out << " |lambdaMax-lambdaMaxOld|/|lambdaMax|: " << Teuchos::ScalarTraits<ST>::magnitude(lambdaMax - lambdaMaxOld) / Teuchos::ScalarTraits<ST>::magnitude(lambdaMax) << endl;
189  }
190  }
191  x->scale(one / norm);
192  }
193  }
194  if (verbose) {
195  *out << " lambdaMax: " << lambdaMax << endl;
196  }
197 
198  norm = x->norm2();
199  if (norm == zero) { // Return something reasonable.
200  if (verbose) {
201  *out << " norm is zero; returning zero" << endl;
202  *out << " Power method terminated after " << numIters << " iterations." << endl;
203  }
204  return zero;
205  }
206  x->scale(one / norm);
207  A.apply(*x, *y);
208  y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
209  lambdaMax = y->norm2();
210  } else {
211  // In this route, the update is as follows:
212  // y=A*x, y = Dinv*y, lambda = x->dot(y), x = y/lambda
213  ST xDinvAx = norm;
214  if (verbose) {
215  *out << " PowerMethod using largest eigenvalue route" << endl;
216  }
217  for (int iter = 0; iter < numIters - 1; ++iter) {
218  if (verbose) {
219  *out << " Iteration " << iter << endl;
220  }
221  A.apply(*x, *y);
222  y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
223 
224  if (((iter + 1) % eigNormalizationFreq == 0) && (iter < numIters - 2)) {
225  xDinvAx = x->dot(*y);
226  if (xDinvAx == zero) { // Return something reasonable.
227  if (verbose) {
228  *out << " xDinvAx is zero; returning zero" << endl;
229  *out << " Power method terminated after " << iter << " iterations." << endl;
230  }
231  return zero;
232  } else {
233  lambdaMaxOld = lambdaMax;
234  lambdaMax = pow(xDinvAx, Teuchos::ScalarTraits<MT>::one() / eigNormalizationFreq);
235  if (Teuchos::ScalarTraits<ST>::magnitude(lambdaMax - lambdaMaxOld) < tolerance * Teuchos::ScalarTraits<ST>::magnitude(lambdaMax)) {
236  if (verbose) {
237  *out << " lambdaMax: " << lambdaMax << endl;
238  *out << " Power method terminated after " << iter << " iterations." << endl;
239  }
240  return lambdaMax;
241  } else if (verbose) {
242  *out << " lambdaMaxOld: " << lambdaMaxOld << endl;
243  *out << " lambdaMax: " << lambdaMax << endl;
244  *out << " |lambdaMax-lambdaMaxOld|/|lambdaMax|: " << Teuchos::ScalarTraits<ST>::magnitude(lambdaMax - lambdaMaxOld) / Teuchos::ScalarTraits<ST>::magnitude(lambdaMax) << endl;
245  }
246  }
247  x->swap(*y);
248  norm = x->norm2();
249  x->scale(one / norm);
250  }
251  }
252  if (verbose) {
253  *out << " lambdaMax: " << lambdaMax << endl;
254  }
255 
256  norm = x->norm2();
257  if (norm == zero) { // Return something reasonable.
258  if (verbose) {
259  *out << " norm is zero; returning zero" << endl;
260  *out << " Power method terminated after " << numIters << " iterations." << endl;
261  }
262  return zero;
263  }
264  x->scale(one / norm);
265  A.apply(*x, *y);
266  y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
267  lambdaMax = y->dot(*x);
268  }
269 
270  if (verbose) {
271  *out << " lambdaMax: " << lambdaMax << endl;
272  *out << " Power method terminated after " << numIters << " iterations." << endl;
273  }
274 
275  return lambdaMax;
276 }
277 
289 template <class V>
290 void computeInitialGuessForPowerMethod(V& x, const bool nonnegativeRealParts) {
291  typedef typename V::device_type::execution_space dev_execution_space;
292  typedef typename V::local_ordinal_type LO;
293 
294  x.randomize();
295 
296  if (nonnegativeRealParts) {
297  auto x_lcl = x.getLocalViewDevice(Tpetra::Access::ReadWrite);
298  auto x_lcl_1d = Kokkos::subview(x_lcl, Kokkos::ALL(), 0);
299 
300  const LO lclNumRows = static_cast<LO>(x.getLocalLength());
301  Kokkos::RangePolicy<dev_execution_space, LO> range(0, lclNumRows);
302  PositivizeVector<decltype(x_lcl_1d), LO> functor(x_lcl_1d);
303  Kokkos::parallel_for(range, functor);
304  }
305 }
306 
323 template <class OperatorType, class V>
324 typename V::scalar_type
325 powerMethod(const OperatorType& A,
326  const V& D_inv,
327  const int numIters,
328  Teuchos::RCP<V> y,
330  const int eigNormalizationFreq = 1,
331  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::null,
332  const bool computeSpectralRadius = true) {
333  typedef typename V::scalar_type ST;
335 
336  bool verbose = (out != Teuchos::null);
337 
338  if (verbose) {
339  *out << "powerMethod:" << std::endl;
340  }
341 
342  const ST zero = static_cast<ST>(0.0);
343 
344  Teuchos::RCP<V> x = Teuchos::rcp(new V(A.getDomainMap()));
345  // For the first pass, just let the pseudorandom number generator
346  // fill x with whatever values it wants; don't try to make its
347  // entries nonnegative.
348  computeInitialGuessForPowerMethod(*x, false);
349 
350  ST lambdaMax = powerMethodWithInitGuess(A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
351 
352  // mfh 07 Jan 2015: Taking the real part here is only a concession
353  // to the compiler, so that this can build with ScalarType =
354  // std::complex<T>. Our Chebyshev implementation only works with
355  // real, symmetric positive definite matrices. The right thing to
356  // do would be what Belos does, which is provide a partial
357  // specialization for ScalarType = std::complex<T> with a stub
358  // implementation (that builds, but whose constructor throws).
359  if (STS::real(lambdaMax) < STS::real(zero)) {
360  if (verbose) {
361  *out << "real(lambdaMax) = " << STS::real(lambdaMax) << " < 0; "
362  "try again with a different random initial guess"
363  << std::endl;
364  }
365  // Max eigenvalue estimate was negative. Perhaps we got unlucky
366  // with the random initial guess. Try again with a different (but
367  // still random) initial guess. Only try again once, so that the
368  // run time is bounded.
369 
370  // For the second pass, make all the entries of the initial guess
371  // vector have nonnegative real parts.
372  computeInitialGuessForPowerMethod(*x, true);
373  lambdaMax = powerMethodWithInitGuess(A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
374  }
375  return lambdaMax;
376 }
377 
378 } // namespace PowerMethod
379 } // namespace Ifpack2
380 
381 #endif // IFPACK2_POWERMETHOD_HPP
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void swap(RCP< T > &r_ptr)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static magnitudeType magnitude(T a)
bool is_null() const