GRASS Programmer's Manual  6.4.2(2012)
eigen_tools.c
Go to the documentation of this file.
00001 #include <grass/gis.h>
00002 #include <math.h>
00003 
00004 
00005 #define MAX_ITERS 30
00006 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
00007 
00008 
00009 int G_tqli(double d[], double e[], int n, double **z)
00010 {
00011     int m, l, iter, i, k;
00012     double s, r, p, g, f, dd, c, b;
00013 
00014     for (i = 1; i < n; i++)
00015         e[i - 1] = e[i];
00016     e[n - 1] = 0.0;
00017     for (l = 0; l < n; l++) {
00018         iter = 0;
00019 
00020         do {
00021             for (m = l; m < n - 1; m++) {
00022                 dd = fabs(d[m]) + fabs(d[m + 1]);
00023                 if (fabs(e[m]) + dd == dd)
00024                     break;
00025             }
00026 
00027             if (m != l) {
00028                 if (iter++ == MAX_ITERS)
00029                     return 0;   /* Too many iterations in TQLI */
00030                 g = (d[l + 1] - d[l]) / (2.0 * e[l]);
00031                 r = sqrt((g * g) + 1.0);
00032                 g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
00033                 s = c = 1.0;
00034                 p = 0.0;
00035 
00036                 for (i = m - 1; i >= l; i--) {
00037                     f = s * e[i];
00038                     b = c * e[i];
00039 
00040                     if (fabs(f) >= fabs(g)) {
00041                         c = g / f;
00042                         r = sqrt((c * c) + 1.0);
00043                         e[i + 1] = f * r;
00044                         c *= (s = 1.0 / r);
00045                     }
00046                     else {
00047                         s = f / g;
00048                         r = sqrt((s * s) + 1.0);
00049                         e[i + 1] = g * r;
00050                         s *= (c = 1.0 / r);
00051                     }
00052 
00053                     g = d[i + 1] - p;
00054                     r = (d[i] - g) * s + 2.0 * c * b;
00055                     p = s * r;
00056                     d[i + 1] = g + p;
00057                     g = c * r - b;
00058 
00059                     /* Next loop can be omitted if eigenvectors not wanted */
00060                     for (k = 0; k < n; k++) {
00061                         f = z[k][i + 1];
00062                         z[k][i + 1] = s * z[k][i] + c * f;
00063                         z[k][i] = c * z[k][i] - s * f;
00064                     }
00065                 }
00066                 d[l] = d[l] - p;
00067                 e[l] = g;
00068                 e[m] = 0.0;
00069             }
00070         } while (m != l);
00071     }
00072 
00073     return 1;
00074 }
00075 
00076 
00077 void G_tred2(double **a, int n, double d[], double e[])
00078 {
00079     int l, k, j, i;
00080     double scale, hh, h, g, f;
00081 
00082     for (i = n - 1; i >= 1; i--) {
00083         l = i - 1;
00084         h = scale = 0.0;
00085 
00086         if (l > 0) {
00087             for (k = 0; k <= l; k++)
00088                 scale += fabs(a[i][k]);
00089 
00090             if (scale == 0.0)
00091                 e[i] = a[i][l];
00092             else {
00093                 for (k = 0; k <= l; k++) {
00094                     a[i][k] /= scale;
00095                     h += a[i][k] * a[i][k];
00096                 }
00097 
00098                 f = a[i][l];
00099                 g = f > 0 ? -sqrt(h) : sqrt(h);
00100                 e[i] = scale * g;
00101                 h -= f * g;
00102                 a[i][l] = f - g;
00103                 f = 0.0;
00104 
00105                 for (j = 0; j <= l; j++) {
00106                     /* Next statement can be omitted if eigenvectors not wanted */
00107                     a[j][i] = a[i][j] / h;
00108                     g = 0.0;
00109                     for (k = 0; k <= j; k++)
00110                         g += a[j][k] * a[i][k];
00111                     for (k = j + 1; k <= l; k++)
00112                         g += a[k][j] * a[i][k];
00113                     e[j] = g / h;
00114                     f += e[j] * a[i][j];
00115                 }
00116 
00117                 hh = f / (h + h);
00118                 for (j = 0; j <= l; j++) {
00119                     f = a[i][j];
00120                     e[j] = g = e[j] - hh * f;
00121 
00122                     for (k = 0; k <= j; k++)
00123                         a[j][k] -= (f * e[k] + g * a[i][k]);
00124                 }
00125             }
00126         }
00127         else
00128             e[i] = a[i][l];
00129         d[i] = h;
00130     }
00131 
00132     /* Next statement can be omitted if eigenvectors not wanted */
00133     d[0] = 0.0;
00134     e[0] = 0.0;
00135 
00136     /* Contents of this loop can be omitted if eigenvectors not
00137        wanted except for statement d[i]=a[i][i]; */
00138     for (i = 0; i < n; i++) {
00139         l = i - 1;
00140 
00141         if (d[i]) {
00142             for (j = 0; j <= l; j++) {
00143                 g = 0.0;
00144                 for (k = 0; k <= l; k++)
00145                     g += a[i][k] * a[k][j];
00146                 for (k = 0; k <= l; k++)
00147                     a[k][j] -= g * a[k][i];
00148             }
00149         }
00150 
00151         d[i] = a[i][i];
00152         a[i][i] = 1.0;
00153         for (j = 0; j <= l; j++)
00154             a[j][i] = a[i][j] = 0.0;
00155     }
00156 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines