42#ifndef THYRA_SILLY_CG_SOLVE_HPP
43#define THYRA_SILLY_CG_SOLVE_HPP
45#include "Thyra_LinearOpBase.hpp"
46#include "Thyra_VectorStdOps.hpp"
47#include "Thyra_AssertOp.hpp"
65 const int maxNumIters,
74 const Scalar one = ST::one(), zero = ST::zero();
using Teuchos::as;
83 std::ios::fmtflags fmt(out.flags());
85 out <<
"\nStarting CG solver ...\n" << std::scientific <<
"\ndescribe A:\n"<<describe(A, vl)
86 <<
"\ndescribe b:\n"<<describe(b, vl)<<
"\ndescribe x:\n"<<describe(*x, vl)<<
"\n";
89 const RCP<const VectorSpaceBase<Scalar> > space = A.
domain();
90 const RCP<VectorBase<Scalar> > r = createMember(space);
92 V_V(r.ptr(), b); apply<Scalar>(A, NOTRANS, *x, r.ptr(), -one, one);
93 const ScalarMag r0_nrm = norm(*r);
94 if (r0_nrm==zero)
return true;
95 const RCP<VectorBase<Scalar> > p = createMember(space), q = createMember(space);
96 Scalar rho_old = -one;
99 for(
int iter = 0; iter <= maxNumIters; ++iter ) {
102 const ScalarMag r_nrm = norm(*r);
103 const bool isConverged = r_nrm/r0_nrm <= tolerance;
104 if( iter%(maxNumIters/10+1) == 0 || iter == maxNumIters || isConverged ) {
105 out <<
"Iter = " << iter <<
", ||b-A*x||/||b-A*x0|| = " << (r_nrm/r0_nrm) << std::endl;
106 if( r_nrm/r0_nrm < tolerance ) {
113 const Scalar rho = inner(*r, *r);
114 if (iter==0) V_V(p.ptr(), *r);
115 else Vp_V( p.ptr(), *r, rho/rho_old );
116 apply<Scalar>(A, NOTRANS, *p, q.ptr());
117 const Scalar alpha = rho/inner(*p, *q);
118 Vp_StV( x, +alpha, *p );
119 Vp_StV( r.ptr(), -alpha, *q );
Base class for all linear operators.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Abstract interface for finite-dimensional dense vectors.
Abstract interface for objects that represent a space for vectors.
#define THYRA_ASSERT_LINEAR_OP_VEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the vector version of...
@ NOTRANS
Use the non-transposed operator.
TypeTo as(const TypeFrom &t)