28 #include "sirf/common/Operator.h"
36 template<
class value_type>
41 void set_num_iterations(
int nit)
47 template<
class vector_type>
49 vector_type& x ,
int verb=0)
52 value_type a[] = { 0, 0, 0, 0 };
59 std::unique_ptr<vector_type> uptr_y = x.clone();
60 vector_type& y = *uptr_y;
61 std::unique_ptr<vector_type> uptr_z = x.clone();
62 vector_type& z = *uptr_z;
63 std::unique_ptr<vector_type> uptr_w = x.clone();
64 vector_type& w = *uptr_w;
65 std::unique_ptr<vector_type> uptr_Az = x.clone();
66 vector_type& Az = *uptr_Az;
70 std::shared_ptr<vector_type> sptr_Ax = A(x);
71 vector_type& Ax = *sptr_Ax;
73 for (
int it = 0; it < nit_; it++) {
75 y.axpby(1.0, Ax, -lmd, x);
78 std::cout <<
'\n' << std::flush;
79 std::cout <<
"CG iteration " << it
80 <<
": largest eigenvalue " << std::real(lmd) <<
'\n';
84 w.axpby(1.0, Az, -t, z);
85 t = w.dot(y) / (lmd - t);
86 y.axpby(1.0, y, t, z);
95 y.axpby(1.0, y, -t, x);
102 std::shared_ptr<vector_type> sptr_Ay = A(y);
103 vector_type& Ay = *sptr_Ay;
110 z.axpby(u[0], x, u[1], y);
111 x.axpby(v[0], x, v[1], y);
115 Az.axpby(u[0], Ax, u[1], Ay);
116 Ax.axpby(v[0], Ax, v[1], Ay);
122 return std::real(lmd);
129 void eigh2_(
const value_type a[4], value_type lmd[2], value_type x[2],
130 value_type y[2])
const
134 value_type c = abs(a[1]);
135 value_type d = abs(a[0] - a[3]);
136 d = d * d + (c + c) * (c + c);
137 value_type s = sqrt(d);
138 value_type t = a[0] + a[3];
139 if (s == zero && t == zero) {
148 lmd[0] = (t - s) / (one + one);
149 lmd[1] = (t + s) / (one + one);
151 value_type q = lmd[0] - a[0];
152 value_type pc = (p == zero ? 0 : abs(p*p) / p);
153 value_type qc = (q == zero ? 0 : abs(q*q) / q);
156 s = sqrt(abs(pc*x[0] + qc*x[1]));
165 template<
typename vector_type,
typename value_type>
168 Wrapped_sptr(std::shared_ptr<vector_type> sptr) : sptr_(sptr) {}
169 std::shared_ptr<vector_type> sptr()
const
173 std::unique_ptr<Wrapped_sptr<vector_type, value_type> > clone()
const
175 std::shared_ptr<vector_type> sptr(sptr_->clone());
176 return std::unique_ptr<Wrapped_sptr>(
new Wrapped_sptr(sptr));
180 return sptr_->norm();
189 void* ptr = (
void*)&s;
190 sptr_->dot(*y.sptr(), ptr);
195 void* ptr_a = (
void*)&a;
196 void* ptr_b = (
void*)&b;
197 sptr_->axpby(ptr_a, *x.sptr(), ptr_b, *y.sptr());
200 std::shared_ptr<vector_type> sptr_;
Definition: JacobiCG.h:37
Definition: JacobiCG.h:166
Abstract data container.
Definition: GeometricalInfo.cpp:141