42#ifndef THYRA_SILLY_POWER_METHOD_HPP
43#define THYRA_SILLY_POWER_METHOD_HPP
45#include "Thyra_LinearOpBase.hpp"
46#include "Thyra_VectorStdOps.hpp"
63 const int maxNumIters,
78 out <<
"\nStarting power method (target tolerance = "<<tolerance<<
") ...\n\n";
79 VectorPtr q = createMember(A.
domain()), z = createMember(A.
range()), r = createMember(A.
range());
80 Thyra::seed_randomize<Scalar>(0);
81 Thyra::randomize( Scalar(-one), Scalar(+one), z.ptr() );
84 for(
int iter = 0; iter < maxNumIters; ++iter ) {
85 const ScalarMag z_nrm = norm(*z);
86 V_StV( q.ptr(), Scalar(one/z_nrm), *z );
87 apply<Scalar>( A, NOTRANS , *q, z.ptr() );
88 *lambda = scalarProd(*q,*z);
89 if( iter%(maxNumIters/10) == 0 || iter+1 == maxNumIters ) {
90 V_StVpV(r.ptr(),Scalar(-*lambda),*q,*z);
91 const ScalarMag r_nrm = norm(*r);
92 out <<
"Iter = " << iter <<
", lambda = " << (*lambda)
93 <<
", ||A*q-lambda*q|| = " << r_nrm << std::endl;
94 if( r_nrm < tolerance )
99 out <<
"\nMaximum number of iterations exceeded with ||-lambda*q + z||"
100 " > tolerence = " << tolerance <<
"\n";
Base class for all linear operators.
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
@ NOTRANS
Use the non-transposed operator.