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