12345678910111213141516171819202122232425262728293031323334353637383940414243444546 |
- #pragma once
- #include "ra/ra.hh"
- template <class A, class B, class X, class W>
- inline int
- cghs(A & a, B const & b, X & x, W & work, double eps)
- {
- using ra::sqr;
- work = 0.;
- auto g = work(0);
- auto r = work(1);
- auto p = work(2);
- double err = sqr(eps)*reduce_sqrm(b);
- mult(a, x, g);
- g = b-g;
- r = g;
- int i;
- for (i=0; reduce_sqrm(g)>err; ++i) {
- mult(a, r, p);
- double rho = reduce_sqrm(p);
- double sig = dot(r, p);
- double tau = dot(g, r);
- double t = tau/sig;
- double gam = (sqr(t)*rho-tau)/tau;
- x += t*r;
- g -= t*p;
- r = r*gam + g;
- }
- return i;
- }
|