# 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].

 [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, double * x, int n)
{
double v[n], c[n], s[n];
for (int i = 0; i < n; i++)
v[i] = H[n*(i + 1) - 1];
for (int k = n - 1; k >= 1; k--) {
double a = H[k*n + 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*n + k - 1];
x[l] -= ykck*v[l] - yksk*a;
v[l] = c[k]*a + s[k]*v[l];
}
a = H[(k - 1)*n + 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/v;
for (int k = 1; k < n; k++) {
double tau2 = x[k];
x[k-1] = c[k]*tau1 - s[k]*tau2;
tau1 = c[k]*tau2 + s[k]*tau1;
}
x[n-1] = tau1;
}