GRASS Programmer's Manual  6.4.4(2014)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
eigen.c
Go to the documentation of this file.
1 /* taken from i.pca */
2 
3 #include <stdlib.h>
4 #include <grass/gmath.h>
5 #include <grass/gis.h>
6 
7 
8 static int egcmp(const void *pa, const void *pb);
9 
10 
26 int eigen(double **M, /* Input matrix */
27  double **Vectors, /* eigen vector matrix -output */
28  double *lambda, /* Output eigenvalues */
29  int n /* Input matrix dimension */
30  )
31 {
32  int i, j;
33  double **a, *e;
34 
35  a = G_alloc_matrix(n, n);
36  e = G_alloc_vector(n);
37 
38  for (i = 0; i < n; i++)
39  for (j = 0; j < n; j++)
40  a[i][j] = M[i][j];
41 
42  G_tred2(a, n, lambda, e);
43  G_tqli(lambda, e, n, a);
44 
45  /* Returns eigenvectors */
46  if (Vectors)
47  for (i = 0; i < n; i++)
48  for (j = 0; j < n; j++)
49  Vectors[i][j] = a[i][j];
50 
51  G_free_matrix(a);
52  G_free_vector(e);
53 
54  return 0;
55 }
56 
57 
71 int egvorder2(double *d, double **z, long bands)
72 {
73  double *buff;
74  double **tmp;
75  int i, j;
76 
77  /* allocate temporary matrix */
78  buff = (double *)G_malloc(bands * (bands + 1) * sizeof(double));
79  tmp = (double **)G_malloc(bands * sizeof(double *));
80  for (i = 0; i < bands; i++)
81  tmp[i] = &buff[i * (bands + 1)];
82 
83  /* concatenate (vertically) z and d into tmp */
84  for (i = 0; i < bands; i++) {
85  for (j = 0; j < bands; j++)
86  tmp[i][j + 1] = z[j][i];
87  tmp[i][0] = d[i];
88  }
89 
90  /* sort the combined matrix */
91  qsort(tmp, bands, sizeof(double *), egcmp);
92 
93  /* split tmp into z and d */
94  for (i = 0; i < bands; i++) {
95  for (j = 0; j < bands; j++)
96  z[j][i] = tmp[i][j + 1];
97  d[i] = tmp[i][0];
98  }
99 
100  /* free temporary matrix */
101  G_free(tmp);
102  G_free(buff);
103 
104  return 0;
105 }
106 
107 
120 int transpose2(double **eigmat, long bands)
121 {
122  int i, j;
123 
124  for (i = 0; i < bands; i++)
125  for (j = 0; j < i; j++) {
126  double tmp = eigmat[i][j];
127 
128  eigmat[i][j] = eigmat[j][i];
129  eigmat[j][i] = tmp;
130  }
131 
132  return 0;
133 }
134 
135 
136 static int egcmp(const void *pa, const void *pb)
137 {
138  const double *a = *(const double *const *)pa;
139  const double *b = *(const double *const *)pb;
140 
141  if (*a > *b)
142  return -1;
143  if (*a < *b)
144  return 1;
145 
146  return 0;
147 }
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
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
int egvorder2(double *d, double **z, long bands)
Returns 0.
Definition: eigen.c:71
void G_free_vector(double *v)
Vector memory deallocation.
Definition: dalloc.c:129
void G_free_matrix(double **m)
Matrix memory deallocation.
Definition: dalloc.c:169
int G_tqli(double d[], double e[], int n, double **z)
Definition: eigen_tools.c:9
void G_tred2(double **a, int n, double d[], double e[])
Definition: eigen_tools.c:77
int transpose2(double **eigmat, long bands)
Returns 0.
Definition: eigen.c:120
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41