GRASS Programmer's Manual
6.4.2(2012)
|
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 }