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)
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;
137 bool verbose = (out != Teuchos::null);
141 *out <<
" powerMethodWithInitGuess:" << endl;
144 const ST zero =
static_cast<ST
> (0.0);
145 const ST one =
static_cast<ST
> (1.0);
147 ST lambdaMaxOld = one;
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.");
160 *out <<
" Original norm1(x): " << x->norm1()
161 <<
", norm2(x): " << norm << endl;
164 x->scale(one / norm);
167 *out <<
" norm1(x.scale(one/norm)): " << x->norm1() << endl;
171 y = Teuchos::rcp(
new V(A.getRangeMap()));
180 if(computeSpectralRadius) {
184 *out <<
" PowerMethod using spectral radius route" << endl;
186 for(
int iter = 0; iter < numIters-1; ++iter) {
188 *out <<
" Iteration " << iter << endl;
191 x->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
193 if(((iter+1) % eigNormalizationFreq == 0) && (iter < numIters-2)) {
197 *out <<
" norm is zero; returning zero" << endl;
198 *out <<
" Power method terminated after "<< iter <<
" iterations." << endl;
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)) {
206 *out <<
" lambdaMax: " << lambdaMax << endl;
207 *out <<
" Power method terminated after "<< iter <<
" iterations." << endl;
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;
216 x->scale(one / norm);
220 *out <<
" lambdaMax: " << lambdaMax << endl;
226 *out <<
" norm is zero; returning zero" << endl;
227 *out <<
" Power method terminated after "<< numIters <<
" iterations." << endl;
231 x->scale(one / norm);
233 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
234 lambdaMax = y->norm2();
240 *out <<
" PowerMethod using largest eigenvalue route" << endl;
242 for (
int iter = 0; iter < numIters-1; ++iter) {
244 *out <<
" Iteration " << iter << endl;
247 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
249 if(((iter+1) % eigNormalizationFreq == 0) && (iter < numIters-2)) {
250 xDinvAx = x->dot(*y);
251 if(xDinvAx == zero) {
253 *out <<
" xDinvAx is zero; returning zero" << endl;
254 *out <<
" Power method terminated after "<< iter <<
" iterations." << endl;
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)) {
262 *out <<
" lambdaMax: " << lambdaMax << endl;
263 *out <<
" Power method terminated after "<< iter <<
" iterations." << endl;
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;
274 x->scale(one / norm);
278 *out <<
" lambdaMax: " << lambdaMax << endl;
284 *out <<
" norm is zero; returning zero" << endl;
285 *out <<
" Power method terminated after "<< numIters <<
" iterations." << endl;
289 x->scale(one / norm);
291 y->elementWiseMultiply(STS::one(), D_inv, *y, STS::zero());
292 lambdaMax = y->dot(*x);
296 *out <<
" lambdaMax: " << lambdaMax << endl;
297 *out <<
" Power method terminated after "<< numIters <<
" iterations." << endl;
319 typedef typename V::device_type::execution_space dev_execution_space;
320 typedef typename V::local_ordinal_type LO;
324 if(nonnegativeRealParts) {
325 auto x_lcl = x.getLocalViewDevice (Tpetra::Access::ReadWrite);
326 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
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);
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)
363 typedef typename V::scalar_type ST;
364 typedef Teuchos::ScalarTraits<typename V::scalar_type> STS;
366 bool verbose = (out != Teuchos::null);
369 *out <<
"powerMethod:" << std::endl;
372 const ST zero =
static_cast<ST
> (0.0);
374 Teuchos::RCP<V> x = Teuchos::rcp(
new V(A.getDomainMap ()));
380 ST lambdaMax =
powerMethodWithInitGuess (A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
389 if(STS::real (lambdaMax) < STS::real (zero)) {
391 *out <<
"real(lambdaMax) = " << STS::real (lambdaMax) <<
" < 0; "
392 "try again with a different random initial guess" << std::endl;
402 lambdaMax =
powerMethodWithInitGuess (A, D_inv, numIters, x, y, tolerance, eigNormalizationFreq, out, computeSpectralRadius);
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