Geografic-Information-System/lib/gmath/eigen_tools.c

#include <grass/gis.h>
#include <math.h>


#define MAX_ITERS 30
#define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))


int G_tqli(double d[], double e[], int n, double **z)
{
    int m, l, iter, i, k;
    double s, r, p, g, f, dd, c, b;

    for (i = 1; i < n; i++)
	e[i - 1] = e[i];
    e[n - 1] = 0.0;
    for (l = 0; l < n; l++) {
	iter = 0;

	do {
	    for (m = l; m < n - 1; m++) {
		dd = fabs(d[m]) + fabs(d[m + 1]);
		if (fabs(e[m]) + dd == dd)
		    break;
	    }

	    if (m != l) {
		if (iter++ == MAX_ITERS)
		    return 0;	/* Too many iterations in TQLI */
		g = (d[l + 1] - d[l]) / (2.0 * e[l]);
		r = sqrt((g * g) + 1.0);
		g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
		s = c = 1.0;
		p = 0.0;

		for (i = m - 1; i >= l; i--) {
		    f = s * e[i];
		    b = c * e[i];

		    if (fabs(f) >= fabs(g)) {
			c = g / f;
			r = sqrt((c * c) + 1.0);
			e[i + 1] = f * r;
			c *= (s = 1.0 / r);
		    }
		    else {
			s = f / g;
			r = sqrt((s * s) + 1.0);
			e[i + 1] = g * r;
			s *= (c = 1.0 / r);
		    }

		    g = d[i + 1] - p;
		    r = (d[i] - g) * s + 2.0 * c * b;
		    p = s * r;
		    d[i + 1] = g + p;
		    g = c * r - b;

		    /* Next loop can be omitted if eigenvectors not wanted */
		    for (k = 0; k < n; k++) {
			f = z[k][i + 1];
			z[k][i + 1] = s * z[k][i] + c * f;
			z[k][i] = c * z[k][i] - s * f;
		    }
		}
		d[l] = d[l] - p;
		e[l] = g;
		e[m] = 0.0;
	    }
	} while (m != l);
    }

    return 1;
}


void G_tred2(double **a, int n, double d[], double e[])
{
    int l, k, j, i;
    double scale, hh, h, g, f;

    for (i = n - 1; i >= 1; i--) {
	l = i - 1;
	h = scale = 0.0;

	if (l > 0) {
	    for (k = 0; k <= l; k++)
		scale += fabs(a[i][k]);

	    if (scale == 0.0)
		e[i] = a[i][l];
	    else {
		for (k = 0; k <= l; k++) {
		    a[i][k] /= scale;
		    h += a[i][k] * a[i][k];
		}

		f = a[i][l];
		g = f > 0 ? -sqrt(h) : sqrt(h);
		e[i] = scale * g;
		h -= f * g;
		a[i][l] = f - g;
		f = 0.0;

		for (j = 0; j <= l; j++) {
		    /* Next statement can be omitted if eigenvectors not wanted */
		    a[j][i] = a[i][j] / h;
		    g = 0.0;
		    for (k = 0; k <= j; k++)
			g += a[j][k] * a[i][k];
		    for (k = j + 1; k <= l; k++)
			g += a[k][j] * a[i][k];
		    e[j] = g / h;
		    f += e[j] * a[i][j];
		}

		hh = f / (h + h);
		for (j = 0; j <= l; j++) {
		    f = a[i][j];
		    e[j] = g = e[j] - hh * f;

		    for (k = 0; k <= j; k++)
			a[j][k] -= (f * e[k] + g * a[i][k]);
		}
	    }
	}
	else
	    e[i] = a[i][l];
	d[i] = h;
    }

    /* Next statement can be omitted if eigenvectors not wanted */
    d[0] = 0.0;
    e[0] = 0.0;

    /* Contents of this loop can be omitted if eigenvectors not
       wanted except for statement d[i]=a[i][i]; */
    for (i = 0; i < n; i++) {
	l = i - 1;

	if (d[i]) {
	    for (j = 0; j <= l; j++) {
		g = 0.0;
		for (k = 0; k <= l; k++)
		    g += a[i][k] * a[k][j];
		for (k = 0; k <= l; k++)
		    a[k][j] -= g * a[k][i];
	    }
	}

	d[i] = a[i][i];
	a[i][i] = 1.0;
	for (j = 0; j <= l; j++)
	    a[j][i] = a[i][j] = 0.0;
    }
}