GRASS Programmer's Manual  6.4.2(2012)
eigen.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines