58 Vector s(m+1), cs(m+1), sn(m+1), w;
60 Real normb =
norm(M.solve(b));
61 Vector r = M.solve(b - A * x);
67 if ((resid =
norm(r) / normb) <= tol) {
73 Vector *v =
new Vector[m+1];
75 while (j <= max_iter) {
76 v[0] = r * (1.0 / beta);
80 for (i = 0; i < m && j <= max_iter; i++, j++) {
81 w = M.solve(A * v[i]);
82 for (k = 0; k <= i; k++) {
83 H(k, i) =
dot(w, v[k]);
87 v[i+1] = w * (1.0 / H(i+1, i));
89 for (k = 0; k < i; k++)
96 if ((resid =
abs(s(i+1)) / normb) < tol) {
104 Update(x, m - 1, H, s, v);
105 r = M.solve(b - A * x);
107 if ((resid = beta / normb) < tol) {
void Update(Vector &x, int k, Matrix &h, Vector &s, Vector v[])
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
A::value_type dot(const A &a, const B &b)
Dot product of array or expression.
A::value_type norm(const A &a)
Euklidean norm of array or expression.