10 #ifndef IFPACK2_POWERMETHOD_HPP
11 #define IFPACK2_POWERMETHOD_HPP
20 #if KOKKOS_VERSION >= 40799
21 #include "KokkosKernels_ArithTraits.hpp"
23 #include "Kokkos_ArithTraits.hpp"
25 #include "Teuchos_FancyOStream.hpp"
26 #include "Teuchos_oblackholestream.hpp"
27 #include "Tpetra_Details_residual.hpp"
32 namespace PowerMethod {
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.");
51 PositivizeVector(
const OneDViewType& x)
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;
60 typedef Kokkos::ArithTraits<IST> STS;
62 #if KOKKOS_VERSION >= 40799
63 typedef KokkosKernels::ArithTraits<typename STS::mag_type> STM;
65 typedef Kokkos::ArithTraits<typename STS::mag_type> STM;
68 if (STS::real(x_(i)) < STM::zero()) {
98 template <
class OperatorType,
class V>
99 typename V::scalar_type
100 powerMethodWithInitGuess(
const OperatorType& A,
106 const int eigNormalizationFreq = 1,
108 const bool computeSpectralRadius =
true) {
109 typedef typename V::scalar_type ST;
113 bool verbose = (out != Teuchos::null);
117 *out <<
" powerMethodWithInitGuess:" << endl;
120 const ST zero =
static_cast<ST
>(0.0);
121 const ST one =
static_cast<ST
>(1.0);
123 ST lambdaMaxOld = one;
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.");
135 *out <<
" Original norm1(x): " << x->norm1()
136 <<
", norm2(x): " << norm << endl;
139 x->scale(one / norm);
142 *out <<
" norm1(x.scale(one/norm)): " << x->norm1() << endl;
155 if (computeSpectralRadius) {
159 *out <<
" PowerMethod using spectral radius route" << endl;
161 for (
int iter = 0; iter < numIters - 1; ++iter) {
163 *out <<
" Iteration " << iter << endl;
166 x->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
168 if (((iter + 1) % eigNormalizationFreq == 0) && (iter < numIters - 2)) {
172 *out <<
" norm is zero; returning zero" << endl;
173 *out <<
" Power method terminated after " << iter <<
" iterations." << endl;
177 lambdaMaxOld = lambdaMax;
181 *out <<
" lambdaMax: " << lambdaMax << endl;
182 *out <<
" Power method terminated after " << iter <<
" iterations." << endl;
185 }
else if (verbose) {
186 *out <<
" lambdaMaxOld: " << lambdaMaxOld << endl;
187 *out <<
" lambdaMax: " << lambdaMax << endl;
191 x->scale(one / norm);
195 *out <<
" lambdaMax: " << lambdaMax << endl;
201 *out <<
" norm is zero; returning zero" << endl;
202 *out <<
" Power method terminated after " << numIters <<
" iterations." << endl;
206 x->scale(one / norm);
208 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
209 lambdaMax = y->norm2();
215 *out <<
" PowerMethod using largest eigenvalue route" << endl;
217 for (
int iter = 0; iter < numIters - 1; ++iter) {
219 *out <<
" Iteration " << iter << endl;
222 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
224 if (((iter + 1) % eigNormalizationFreq == 0) && (iter < numIters - 2)) {
225 xDinvAx = x->dot(*y);
226 if (xDinvAx == zero) {
228 *out <<
" xDinvAx is zero; returning zero" << endl;
229 *out <<
" Power method terminated after " << iter <<
" iterations." << endl;
233 lambdaMaxOld = lambdaMax;
237 *out <<
" lambdaMax: " << lambdaMax << endl;
238 *out <<
" Power method terminated after " << iter <<
" iterations." << endl;
241 }
else if (verbose) {
242 *out <<
" lambdaMaxOld: " << lambdaMaxOld << endl;
243 *out <<
" lambdaMax: " << lambdaMax << endl;
249 x->scale(one / norm);
253 *out <<
" lambdaMax: " << lambdaMax << endl;
259 *out <<
" norm is zero; returning zero" << endl;
260 *out <<
" Power method terminated after " << numIters <<
" iterations." << endl;
264 x->scale(one / norm);
266 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
267 lambdaMax = y->dot(*x);
271 *out <<
" lambdaMax: " << lambdaMax << endl;
272 *out <<
" Power method terminated after " << numIters <<
" iterations." << endl;
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;
296 if (nonnegativeRealParts) {
297 auto x_lcl = x.getLocalViewDevice(Tpetra::Access::ReadWrite);
298 auto x_lcl_1d = Kokkos::subview(x_lcl, Kokkos::ALL(), 0);
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);
323 template <
class OperatorType,
class V>
324 typename V::scalar_type
325 powerMethod(
const OperatorType& A,
330 const int eigNormalizationFreq = 1,
332 const bool computeSpectralRadius =
true) {
333 typedef typename V::scalar_type ST;
336 bool verbose = (out != Teuchos::null);
339 *out <<
"powerMethod:" << std::endl;
342 const ST zero =
static_cast<ST
>(0.0);
348 computeInitialGuessForPowerMethod(*x,
false);
350 ST lambdaMax = powerMethodWithInitGuess(A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
359 if (STS::real(lambdaMax) < STS::real(zero)) {
361 *out <<
"real(lambdaMax) = " << STS::real(lambdaMax) <<
" < 0; "
362 "try again with a different random initial guess"
372 computeInitialGuessForPowerMethod(*x,
true);
373 lambdaMax = powerMethodWithInitGuess(A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
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)