GRASS Programmer's Manual
6.4.2(2012)
|
00001 /* taken from i.pca */ 00002 00003 #include <stdlib.h> 00004 #include <grass/gmath.h> 00005 #include <grass/gis.h> 00006 00007 00008 static int egcmp(const void *pa, const void *pb); 00009 00010 00026 int eigen(double **M, /* Input matrix */ 00027 double **Vectors, /* eigen vector matrix -output */ 00028 double *lambda, /* Output eigenvalues */ 00029 int n /* Input matrix dimension */ 00030 ) 00031 { 00032 int i, j; 00033 double **a, *e; 00034 00035 a = G_alloc_matrix(n, n); 00036 e = G_alloc_vector(n); 00037 00038 for (i = 0; i < n; i++) 00039 for (j = 0; j < n; j++) 00040 a[i][j] = M[i][j]; 00041 00042 G_tred2(a, n, lambda, e); 00043 G_tqli(lambda, e, n, a); 00044 00045 /* Returns eigenvectors */ 00046 if (Vectors) 00047 for (i = 0; i < n; i++) 00048 for (j = 0; j < n; j++) 00049 Vectors[i][j] = a[i][j]; 00050 00051 G_free_matrix(a); 00052 G_free_vector(e); 00053 00054 return 0; 00055 } 00056 00057 00071 int egvorder2(double *d, double **z, long bands) 00072 { 00073 double *buff; 00074 double **tmp; 00075 int i, j; 00076 00077 /* allocate temporary matrix */ 00078 buff = (double *)G_malloc(bands * (bands + 1) * sizeof(double)); 00079 tmp = (double **)G_malloc(bands * sizeof(double *)); 00080 for (i = 0; i < bands; i++) 00081 tmp[i] = &buff[i * (bands + 1)]; 00082 00083 /* concatenate (vertically) z and d into tmp */ 00084 for (i = 0; i < bands; i++) { 00085 for (j = 0; j < bands; j++) 00086 tmp[i][j + 1] = z[j][i]; 00087 tmp[i][0] = d[i]; 00088 } 00089 00090 /* sort the combined matrix */ 00091 qsort(tmp, bands, sizeof(double *), egcmp); 00092 00093 /* split tmp into z and d */ 00094 for (i = 0; i < bands; i++) { 00095 for (j = 0; j < bands; j++) 00096 z[j][i] = tmp[i][j + 1]; 00097 d[i] = tmp[i][0]; 00098 } 00099 00100 /* free temporary matrix */ 00101 G_free(tmp); 00102 G_free(buff); 00103 00104 return 0; 00105 } 00106 00107 00120 int transpose2(double **eigmat, long bands) 00121 { 00122 int i, j; 00123 00124 for (i = 0; i < bands; i++) 00125 for (j = 0; j < i; j++) { 00126 double tmp = eigmat[i][j]; 00127 00128 eigmat[i][j] = eigmat[j][i]; 00129 eigmat[j][i] = tmp; 00130 } 00131 00132 return 0; 00133 } 00134 00135 00136 static int egcmp(const void *pa, const void *pb) 00137 { 00138 const double *a = *(const double *const *)pa; 00139 const double *b = *(const double *const *)pb; 00140 00141 if (*a > *b) 00142 return -1; 00143 if (*a < *b) 00144 return 1; 00145 00146 return 0; 00147 }