4
#include <grass/gmath.h>
8
static int egcmp(const void *pa, const void *pb);
12
* \fn int eigen (double **M, double **Vectors, double *lambda, int n)
14
* \brief Computes eigenvalues (and eigen vectors if desired) for
17
* Computes eigenvalues (and eigen vectors if desired) for symmetric matices.
19
* \param M Input matrix
20
* \param Vectors eigen output vector matrix
21
* \param lambda Output eigenvalues
22
* \param n Input matrix dimension
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 */
35
a = G_alloc_matrix(n, n);
36
e = G_alloc_vector(n);
38
for (i = 0; i < n; i++)
39
for (j = 0; j < n; j++)
42
G_tred2(a, n, lambda, e);
43
G_tqli(lambda, e, n, a);
45
/* Returns eigenvectors */
47
for (i = 0; i < n; i++)
48
for (j = 0; j < n; j++)
49
Vectors[i][j] = a[i][j];
59
* \fn int egvorder2 (double *d, double **z, long bands)
71
int egvorder2(double *d, double **z, long bands)
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)];
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];
90
/* sort the combined matrix */
91
qsort(tmp, bands, sizeof(double *), egcmp);
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];
100
/* free temporary matrix */
109
* \fn int transpose2 (double **eigmat, long bands)
120
int transpose2(double **eigmat, long bands)
124
for (i = 0; i < bands; i++)
125
for (j = 0; j < i; j++) {
126
double tmp = eigmat[i][j];
128
eigmat[i][j] = eigmat[j][i];
136
static int egcmp(const void *pa, const void *pb)
138
const double *a = *(const double *const *)pa;
139
const double *b = *(const double *const *)pb;