GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <stdlib.h> 00002 #include <math.h> 00003 #include <grass/gis.h> 00004 #include <grass/gmath.h> 00005 00006 00007 /***************************************************************************/ 00008 00009 /* this does not use the Jacobi method, but it should give the same result */ 00010 00011 int jacobi(double a[MX][MX], long n, double d[MX], double v[MX][MX]) 00012 { 00013 double *aa[MX], *vv[MX], *dd; 00014 int i; 00015 00016 for (i = 0; i < n; i++) { 00017 aa[i] = &a[i + 1][1]; 00018 vv[i] = &v[i + 1][1]; 00019 } 00020 dd = &d[1]; 00021 eigen(aa, vv, dd, n); 00022 00023 return 0; 00024 } 00025 00026 /***************************************************************************/ 00027 00028 static int egcmp(const void *pa, const void *pb) 00029 { 00030 const double *a = *(const double *const *)pa; 00031 const double *b = *(const double *const *)pb; 00032 00033 if (*a > *b) 00034 return -1; 00035 if (*a < *b) 00036 return 1; 00037 00038 return 0; 00039 } 00040 00041 int egvorder(double d[MX], double z[MX][MX], long bands) 00042 { 00043 double *buff; 00044 double **tmp; 00045 int i, j; 00046 00047 /* allocate temporary matrix */ 00048 00049 buff = (double *)G_malloc(bands * (bands + 1) * sizeof(double)); 00050 tmp = (double **)G_malloc(bands * sizeof(double *)); 00051 for (i = 0; i < bands; i++) 00052 tmp[i] = &buff[i * (bands + 1)]; 00053 00054 /* concatenate (vertically) z and d into tmp */ 00055 00056 for (i = 0; i < bands; i++) { 00057 for (j = 0; j < bands; j++) 00058 tmp[i][j + 1] = z[j + 1][i + 1]; 00059 tmp[i][0] = d[i + 1]; 00060 } 00061 00062 /* sort the combined matrix */ 00063 00064 qsort(tmp, bands, sizeof(double *), egcmp); 00065 00066 /* split tmp into z and d */ 00067 00068 for (i = 0; i < bands; i++) { 00069 for (j = 0; j < bands; j++) 00070 z[j + 1][i + 1] = tmp[i][j + 1]; 00071 d[i + 1] = tmp[i][0]; 00072 } 00073 00074 /* free temporary matrix */ 00075 00076 G_free(tmp); 00077 G_free(buff); 00078 00079 return 0; 00080 } 00081 00082 /***************************************************************************/ 00083 00084 int transpose(double eigmat[MX][MX], long bands) 00085 { 00086 int i, j; 00087 00088 for (i = 1; i <= bands; i++) 00089 for (j = 1; j < i; j++) { 00090 double tmp = eigmat[i][j]; 00091 00092 eigmat[i][j] = eigmat[j][i]; 00093 eigmat[j][i] = tmp; 00094 } 00095 00096 return 0; 00097 } 00098 00099 /***************************************************************************/