GRASS Programmer's Manual  6.4.4(2014)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
jacobi.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <math.h>
3 #include <grass/gis.h>
4 #include <grass/gmath.h>
5 
6 
7 /***************************************************************************/
8 
9 /* this does not use the Jacobi method, but it should give the same result */
10 
11 int jacobi(double a[MX][MX], long n, double d[MX], double v[MX][MX])
12 {
13  double *aa[MX], *vv[MX], *dd;
14  int i;
15 
16  for (i = 0; i < n; i++) {
17  aa[i] = &a[i + 1][1];
18  vv[i] = &v[i + 1][1];
19  }
20  dd = &d[1];
21  eigen(aa, vv, dd, n);
22 
23  return 0;
24 }
25 
26 /***************************************************************************/
27 
28 static int egcmp(const void *pa, const void *pb)
29 {
30  const double *a = *(const double *const *)pa;
31  const double *b = *(const double *const *)pb;
32 
33  if (*a > *b)
34  return -1;
35  if (*a < *b)
36  return 1;
37 
38  return 0;
39 }
40 
41 int egvorder(double d[MX], double z[MX][MX], long bands)
42 {
43  double *buff;
44  double **tmp;
45  int i, j;
46 
47  /* allocate temporary matrix */
48 
49  buff = (double *)G_malloc(bands * (bands + 1) * sizeof(double));
50  tmp = (double **)G_malloc(bands * sizeof(double *));
51  for (i = 0; i < bands; i++)
52  tmp[i] = &buff[i * (bands + 1)];
53 
54  /* concatenate (vertically) z and d into tmp */
55 
56  for (i = 0; i < bands; i++) {
57  for (j = 0; j < bands; j++)
58  tmp[i][j + 1] = z[j + 1][i + 1];
59  tmp[i][0] = d[i + 1];
60  }
61 
62  /* sort the combined matrix */
63 
64  qsort(tmp, bands, sizeof(double *), egcmp);
65 
66  /* split tmp into z and d */
67 
68  for (i = 0; i < bands; i++) {
69  for (j = 0; j < bands; j++)
70  z[j + 1][i + 1] = tmp[i][j + 1];
71  d[i + 1] = tmp[i][0];
72  }
73 
74  /* free temporary matrix */
75 
76  G_free(tmp);
77  G_free(buff);
78 
79  return 0;
80 }
81 
82 /***************************************************************************/
83 
84 int transpose(double eigmat[MX][MX], long bands)
85 {
86  int i, j;
87 
88  for (i = 1; i <= bands; i++)
89  for (j = 1; j < i; j++) {
90  double tmp = eigmat[i][j];
91 
92  eigmat[i][j] = eigmat[j][i];
93  eigmat[j][i] = tmp;
94  }
95 
96  return 0;
97 }
98 
99 /***************************************************************************/
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
float b
Definition: named_colr.c:8
int eigen(double **M, double **Vectors, double *lambda, int n)
Computes eigenvalues (and eigen vectors if desired) for symmetric matices.
Definition: eigen.c:26
int transpose(double eigmat[MX][MX], long bands)
Definition: jacobi.c:84
int jacobi(double a[MX][MX], long n, double d[MX], double v[MX][MX])
Definition: jacobi.c:11
int egvorder(double d[MX], double z[MX][MX], long bands)
Definition: jacobi.c:41