+template <class Vector >
+void
+Update(::Vector<double> &x, int k, FullMatrix<double> &h,
+ ::Vector<double> &s, vector<Vector> &v)
+{
+ ::Vector<double> y(s);
+
+ // Backsolve:
+ for (int i = k; i >= 0; i--) {
+ y(i) /= h(i,i);
+ for (int j = i - 1; j >= 0; j--)
+ y(j) -= h(j,i) * y(i);
+ }
+
+ for (int j = 0; j <= k; j++)
+ x.add (y(j), v[j]);
+}
+
+
+
+
+template < class Matrix, class Vector, class Preconditioner>
+int
+solve2(const Matrix &A, Vector &x, const Vector &b,
+ const Preconditioner &M, int &m, int &max_iter,
+ double &tol)
+{
+ FullMatrix<double> H(m+1,m);
+
+ double resid;
+ int i, j = 1, k;
+ ::Vector<double> s(m+1), cs(m+1), sn(m+1);
+ Vector w(x.size());
+
+ Vector tmp (b.size());
+ M (tmp, b);
+ double normb = tmp.l2_norm();
+ Vector r(b.size());
+ A.residual (tmp, x, b);
+ M(r, tmp);
+ double beta = r.l2_norm();
+
+ if (normb == 0.0)
+ normb = 1;
+
+ if ((resid = r.l2_norm() / normb) <= tol) {
+ tol = resid;
+ max_iter = 0;
+ return 0;
+ }
+
+ vector<Vector> v(m+1);
+
+ while (j <= max_iter) {
+ v[0] = r;
+ v[0].scale (1.0 / beta); // ??? r / beta
+ s = 0.0;
+ s(0) = beta;
+
+ for (i = 0; i < m && j <= max_iter; i++, j++) {
+ A.vmult (tmp, v[i]);
+ M(w, tmp);
+ for (k = 0; k <= i; k++) {
+ H(k, i) = w * v[k];
+ w.sadd (-H(k, i), v[k]);
+ }
+ H(i+1, i) = w.l2_norm();
+ v[i+1] = w;
+ v[i+1].scale (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
+
+ for (k = 0; k < i; k++)
+ ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
+
+ GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
+ ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
+ ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
+
+ if ((resid = fabs(s(i+1)) / normb) < tol) {
+ Update(x, i, H, s, v);
+ tol = resid;
+ max_iter = j;
+ return 0;
+ }
+ }
+ Update(x, i - 1, H, s, v);
+ A.residual (tmp, x, b);
+ M(r, tmp);
+ beta = r.l2_norm();
+ if ((resid = beta / normb) < tol) {
+ tol = resid;
+ max_iter = j;
+ return 0;
+ }
+ }
+
+ tol = resid;
+ return 1;
+}
+
+
+template<class Real>
+void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
+{
+ if (dy == 0.0) {
+ cs = 1.0;
+ sn = 0.0;
+ } else if (fabs(dy) > fabs(dx)) {
+ Real temp = dx / dy;
+ sn = 1.0 / sqrt( 1.0 + temp*temp );
+ cs = temp * sn;
+ } else {
+ Real temp = dy / dx;
+ cs = 1.0 / sqrt( 1.0 + temp*temp );
+ sn = temp * cs;
+ }
+}
+
+
+template<class Real>
+void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
+{
+ Real temp = cs * dx + sn * dy;
+ dy = -sn * dx + cs * dy;
+ dx = temp;
+}
+
+
+
+
+
/* --------------------- Inline and template functions ------------------- */