Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_PowerMethod.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_POWERMETHOD_HPP
45#define IFPACK2_POWERMETHOD_HPP
46
53
54#include "Kokkos_ArithTraits.hpp"
55#include "Teuchos_FancyOStream.hpp"
56#include "Teuchos_oblackholestream.hpp"
57#include "Tpetra_Details_residual.hpp"
58#include <cmath>
59#include <iostream>
60
61namespace Ifpack2 {
62namespace PowerMethod {
63
64namespace { // (anonymous)
65
66// Functor for making sure the real parts of all entries of a vector
67// are nonnegative. We use this in computeInitialGuessForPowerMethod
68// below.
69template<class OneDViewType,
70 class LocalOrdinal = typename OneDViewType::size_type>
71class PositivizeVector {
72 static_assert (Kokkos::is_view<OneDViewType>::value,
73 "OneDViewType must be a 1-D Kokkos::View.");
74 static_assert (static_cast<int> (OneDViewType::rank) == 1,
75 "This functor only works with a 1-D View.");
76 static_assert (std::is_integral<LocalOrdinal>::value,
77 "The second template parameter, LocalOrdinal, "
78 "must be an integer type.");
79public:
80 PositivizeVector (const OneDViewType& x) : x_ (x) {}
81
82 KOKKOS_INLINE_FUNCTION void
83 operator () (const LocalOrdinal& i) const
84 {
85 typedef typename OneDViewType::non_const_value_type IST;
86 typedef Kokkos::Details::ArithTraits<IST> STS;
87 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
88
89 if(STS::real (x_(i)) < STM::zero ()) {
90 x_(i) = -x_(i);
91 }
92 }
93
94private:
95 OneDViewType x_;
96};
97
98} // namespace (anonymous)
99
100
101
121template<class OperatorType, class V>
122typename V::scalar_type
123powerMethodWithInitGuess(const OperatorType& A,
124 const V& D_inv,
125 const int numIters,
126 Teuchos::RCP<V> x,
127 Teuchos::RCP<V> y,
128 const typename Teuchos::ScalarTraits<typename V::scalar_type>::magnitudeType tolerance = 1e-7,
129 const int eigNormalizationFreq = 1,
130 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::null,
131 const bool computeSpectralRadius = true)
132{
133 typedef typename V::scalar_type ST;
134 typedef Teuchos::ScalarTraits<typename V::scalar_type> STS;
135 typedef typename Teuchos::ScalarTraits<typename V::scalar_type>::magnitudeType MT;
136
137 bool verbose = (out != Teuchos::null);
138
139 using std::endl;
140 if(verbose) {
141 *out << " powerMethodWithInitGuess:" << endl;
142 }
143
144 const ST zero = static_cast<ST> (0.0);
145 const ST one = static_cast<ST> (1.0);
146 ST lambdaMax = zero;
147 ST lambdaMaxOld = one;
148 ST norm;
149
150 norm = x->norm2();
151 TEUCHOS_TEST_FOR_EXCEPTION
152 (norm == zero, std::runtime_error,
153 "Ifpack2::PowerMethod::powerMethodWithInitGuess: The initial guess "
154 "has zero norm. This could be either because Tpetra::Vector::"
155 "randomize filled the vector with zeros (if that was used to "
156 "compute the initial guess), or because the norm2 method has a "
157 "bug. The first is not impossible, but unlikely.");
158
159 if(verbose) {
160 *out << " Original norm1(x): " << x->norm1()
161 << ", norm2(x): " << norm << endl;
162 }
163
164 x->scale(one / norm);
165
166 if(verbose) {
167 *out << " norm1(x.scale(one/norm)): " << x->norm1() << endl;
168 }
169
170 if(y.is_null())
171 y = Teuchos::rcp(new V(A.getRangeMap()));
172
173 // GH 04 Nov 2022: The previous implementation of the power method was inconsistent with itself.
174 // It computed x->norm2() inside the loop, but returned x^T*Ax.
175 // This returned horribly incorrect estimates for Chebyshev spectral radius in certain
176 // cases, such as the Cheby_belos test, which has two dominant eigenvalues of 12.5839, -10.6639.
177 // In such case, the power method iteration history would appear to converge to something
178 // between 10 and 12, but the resulting spectral radius estimate was sometimes negative.
179 // These have now been split into two routes which behave consistently with themselves.
180 if(computeSpectralRadius) {
181 // In this route, the update is as follows:
182 // y=A*x, x = Dinv*y, lambda = norm(x), x = x/lambda
183 if(verbose) {
184 *out << " PowerMethod using spectral radius route" << endl;
185 }
186 for(int iter = 0; iter < numIters-1; ++iter) {
187 if(verbose) {
188 *out << " Iteration " << iter << endl;
189 }
190 A.apply(*x, *y);
191 x->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
192
193 if(((iter+1) % eigNormalizationFreq == 0) && (iter < numIters-2)) {
194 norm = x->norm2();
195 if(norm == zero) { // Return something reasonable.
196 if(verbose) {
197 *out << " norm is zero; returning zero" << endl;
198 *out << " Power method terminated after "<< iter << " iterations." << endl;
199 }
200 return zero;
201 } else {
202 lambdaMaxOld = lambdaMax;
203 lambdaMax = pow(norm, Teuchos::ScalarTraits<MT>::one() / eigNormalizationFreq);
204 if(Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld) < tolerance * Teuchos::ScalarTraits<ST>::magnitude(lambdaMax)) {
205 if(verbose) {
206 *out << " lambdaMax: " << lambdaMax << endl;
207 *out << " Power method terminated after "<< iter << " iterations." << endl;
208 }
209 return lambdaMax;
210 } else if(verbose) {
211 *out << " lambdaMaxOld: " << lambdaMaxOld << endl;
212 *out << " lambdaMax: " << lambdaMax << endl;
213 *out << " |lambdaMax-lambdaMaxOld|/|lambdaMax|: " << Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld)/Teuchos::ScalarTraits<ST>::magnitude(lambdaMax) << endl;
214 }
215 }
216 x->scale(one / norm);
217 }
218 }
219 if(verbose) {
220 *out << " lambdaMax: " << lambdaMax << endl;
221 }
222
223 norm = x->norm2();
224 if(norm == zero) { // Return something reasonable.
225 if(verbose) {
226 *out << " norm is zero; returning zero" << endl;
227 *out << " Power method terminated after "<< numIters << " iterations." << endl;
228 }
229 return zero;
230 }
231 x->scale(one / norm);
232 A.apply(*x, *y);
233 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
234 lambdaMax = y->norm2();
235 } else {
236 // In this route, the update is as follows:
237 // y=A*x, y = Dinv*y, lambda = x->dot(y), x = y/lambda
238 ST xDinvAx = norm;
239 if(verbose) {
240 *out << " PowerMethod using largest eigenvalue route" << endl;
241 }
242 for (int iter = 0; iter < numIters-1; ++iter) {
243 if(verbose) {
244 *out << " Iteration " << iter << endl;
245 }
246 A.apply(*x, *y);
247 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
248
249 if(((iter+1) % eigNormalizationFreq == 0) && (iter < numIters-2)) {
250 xDinvAx = x->dot(*y);
251 if(xDinvAx == zero) { // Return something reasonable.
252 if(verbose) {
253 *out << " xDinvAx is zero; returning zero" << endl;
254 *out << " Power method terminated after "<< iter << " iterations." << endl;
255 }
256 return zero;
257 } else {
258 lambdaMaxOld = lambdaMax;
259 lambdaMax = pow(xDinvAx, Teuchos::ScalarTraits<MT>::one() / eigNormalizationFreq);
260 if(Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld) < tolerance * Teuchos::ScalarTraits<ST>::magnitude(lambdaMax)) {
261 if(verbose) {
262 *out << " lambdaMax: " << lambdaMax << endl;
263 *out << " Power method terminated after "<< iter << " iterations." << endl;
264 }
265 return lambdaMax;
266 } else if(verbose) {
267 *out << " lambdaMaxOld: " << lambdaMaxOld << endl;
268 *out << " lambdaMax: " << lambdaMax << endl;
269 *out << " |lambdaMax-lambdaMaxOld|/|lambdaMax|: " << Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld)/Teuchos::ScalarTraits<ST>::magnitude(lambdaMax) << endl;
270 }
271 }
272 x->swap(*y);
273 norm = x->norm2();
274 x->scale(one / norm);
275 }
276 }
277 if(verbose) {
278 *out << " lambdaMax: " << lambdaMax << endl;
279 }
280
281 norm = x->norm2();
282 if(norm == zero) { // Return something reasonable.
283 if(verbose) {
284 *out << " norm is zero; returning zero" << endl;
285 *out << " Power method terminated after "<< numIters << " iterations." << endl;
286 }
287 return zero;
288 }
289 x->scale(one / norm);
290 A.apply(*x, *y);
291 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
292 lambdaMax = y->dot(*x);
293 }
294
295 if(verbose) {
296 *out << " lambdaMax: " << lambdaMax << endl;
297 *out << " Power method terminated after "<< numIters << " iterations." << endl;
298 }
299
300 return lambdaMax;
301}
302
303
315template<class V>
316void
317computeInitialGuessForPowerMethod (V& x, const bool nonnegativeRealParts)
318{
319 typedef typename V::device_type::execution_space dev_execution_space;
320 typedef typename V::local_ordinal_type LO;
321
322 x.randomize ();
323
324 if(nonnegativeRealParts) {
325 auto x_lcl = x.getLocalViewDevice (Tpetra::Access::ReadWrite);
326 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
327
328 const LO lclNumRows = static_cast<LO> (x.getLocalLength ());
329 Kokkos::RangePolicy<dev_execution_space, LO> range (0, lclNumRows);
330 PositivizeVector<decltype (x_lcl_1d), LO> functor (x_lcl_1d);
331 Kokkos::parallel_for (range, functor);
332 }
333}
334
335
352template<class OperatorType, class V>
353typename V::scalar_type
354powerMethod(const OperatorType& A,
355 const V& D_inv,
356 const int numIters,
357 Teuchos::RCP<V> y,
358 const typename Teuchos::ScalarTraits<typename V::scalar_type>::magnitudeType tolerance = 1e-7,
359 const int eigNormalizationFreq = 1,
360 Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::null,
361 const bool computeSpectralRadius = true)
362{
363 typedef typename V::scalar_type ST;
364 typedef Teuchos::ScalarTraits<typename V::scalar_type> STS;
365
366 bool verbose = (out != Teuchos::null);
367
368 if(verbose) {
369 *out << "powerMethod:" << std::endl;
370 }
371
372 const ST zero = static_cast<ST> (0.0);
373
374 Teuchos::RCP<V> x = Teuchos::rcp(new V(A.getDomainMap ()));
375 // For the first pass, just let the pseudorandom number generator
376 // fill x with whatever values it wants; don't try to make its
377 // entries nonnegative.
379
380 ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
381
382 // mfh 07 Jan 2015: Taking the real part here is only a concession
383 // to the compiler, so that this can build with ScalarType =
384 // std::complex<T>. Our Chebyshev implementation only works with
385 // real, symmetric positive definite matrices. The right thing to
386 // do would be what Belos does, which is provide a partial
387 // specialization for ScalarType = std::complex<T> with a stub
388 // implementation (that builds, but whose constructor throws).
389 if(STS::real (lambdaMax) < STS::real (zero)) {
390 if(verbose) {
391 *out << "real(lambdaMax) = " << STS::real (lambdaMax) << " < 0; "
392 "try again with a different random initial guess" << std::endl;
393 }
394 // Max eigenvalue estimate was negative. Perhaps we got unlucky
395 // with the random initial guess. Try again with a different (but
396 // still random) initial guess. Only try again once, so that the
397 // run time is bounded.
398
399 // For the second pass, make all the entries of the initial guess
400 // vector have nonnegative real parts.
402 lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
403 }
404 return lambdaMax;
405}
406
407} // namespace PowerMethod
408} // namespace Ifpack2
409
410#endif // IFPACK2_POWERMETHOD_HPP
V::scalar_type powerMethodWithInitGuess(const OperatorType &A, const V &D_inv, const int numIters, Teuchos::RCP< V > x, Teuchos::RCP< V > y, const typename Teuchos::ScalarTraits< typename V::scalar_type >::magnitudeType tolerance=1e-7, const int eigNormalizationFreq=1, Teuchos::RCP< Teuchos::FancyOStream > out=Teuchos::null, const bool computeSpectralRadius=true)
Use the power method to estimate the maximum eigenvalue of A*D_inv, given an initial guess vector x.
Definition: Ifpack2_PowerMethod.hpp:123
V::scalar_type powerMethod(const OperatorType &A, const V &D_inv, const int numIters, Teuchos::RCP< V > y, const typename Teuchos::ScalarTraits< typename V::scalar_type >::magnitudeType tolerance=1e-7, const int eigNormalizationFreq=1, Teuchos::RCP< Teuchos::FancyOStream > out=Teuchos::null, const bool computeSpectralRadius=true)
Use the power method to estimate the maximum eigenvalue of A*D_inv.
Definition: Ifpack2_PowerMethod.hpp:354
void computeInitialGuessForPowerMethod(V &x, const bool nonnegativeRealParts)
Fill x with random initial guess for power method.
Definition: Ifpack2_PowerMethod.hpp:317
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74