src/hessenberg.h

    A solver for Hessenberg systems

    An Hessenberg matrix is an “almost triangular” matrix i.e. the sum of a triangular matrix and a tridiagonal matrix.

    The function below solves Hx=b where H is an upper Hessenberg matrix of rank n. The right-hand side b is given as vector x and is replaced by the solution. H is given as a one dimensional array where each matrix element is indexed as H_{ij} = H[in+j].

    References

    [henry1994]

    Greg Henry. The shifted Hessenberg system solve computation. Cornell Theory Center, Cornell University, 1994. [ .pdf ]

    		 
    static inline void givens (double x, double y, double * c, double * s)
{
#if 0
  #define sign2(a,b) ((b) > 0. ? ((a) > 0. ? (a) : -(a)) : ((a) > 0. ? -(a) : (a)))

  if (x == 0. && y == 0.)
    *c = 1., *s = 0.;
  else if (fabs(y) > fabs(x)) {
    double t = x/y;
    x = sqrt(1. + t*t);
    *s = - sign2(1./x, y);
    *c = t*(*s);
  }
  else {
    double t = y/x;
    y = sqrt2(1. + t*t);
    *c = sign2(1./y, x);
    *s = - t*(*c);
  }
#else
  double t = sqrt (sq(x) + sq(y));
  *c = x/t, *s = -y/t;
#endif
}

void solve_hessenberg (double H[nl*nl], double x[nl])
{
  double v[nl], c[nl], s[nl];
  for (int i = 0; i < nl; i++)
    v[i] = H[nl*(i + 1) - 1];
  for (int k = nl - 1; k >= 1; k--) {
    double a = H[k*nl + k - 1];
    givens (v[k], a, &c[k], &s[k]);
    x[k] /= c[k]*v[k] - s[k]*a;
    double ykck = x[k]*c[k], yksk = x[k]*s[k];
    for (int l = 0; l <= k - 2; l++) {
      a = H[l*nl + k - 1];
      x[l] -= ykck*v[l] - yksk*a;
      v[l] = c[k]*a + s[k]*v[l];
    }
    a = H[(k - 1)*nl + k - 1];
    x[k-1] -= ykck*v[k-1] - yksk*a;
    v[k-1] = c[k]*a + s[k]*v[k-1];
  }
  double tau1 = x[0]/v[0];
  for (int k = 1; k < nl; k++) {
    double tau2 = x[k];
    x[k-1] = c[k]*tau1 - s[k]*tau2;
    tau1 = c[k]*tau2 + s[k]*tau1;
  }
  x[nl-1] = tau1;
}