4 #include <grass/gmath.h>
11 int jacobi(
double a[MX][MX],
long n,
double d[MX],
double v[MX][MX])
13 double *aa[MX], *vv[MX], *dd;
16 for (i = 0; i < n; i++) {
28 static int egcmp(
const void *pa,
const void *pb)
30 const double *a = *(
const double *
const *)pa;
31 const double *
b = *(
const double *
const *)pb;
41 int egvorder(
double d[MX],
double z[MX][MX],
long bands)
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)];
56 for (i = 0; i < bands; i++) {
57 for (j = 0; j < bands; j++)
58 tmp[i][j + 1] = z[j + 1][i + 1];
64 qsort(tmp, bands,
sizeof(
double *), egcmp);
68 for (i = 0; i < bands; i++) {
69 for (j = 0; j < bands; j++)
70 z[j + 1][i + 1] = tmp[i][j + 1];
88 for (i = 1; i <= bands; i++)
89 for (j = 1; j < i; j++) {
90 double tmp = eigmat[i][j];
92 eigmat[i][j] = eigmat[j][i];
void G_free(void *buf)
Free allocated memory.
int eigen(double **M, double **Vectors, double *lambda, int n)
Computes eigenvalues (and eigen vectors if desired) for symmetric matices.
int transpose(double eigmat[MX][MX], long bands)
int jacobi(double a[MX][MX], long n, double d[MX], double v[MX][MX])
int egvorder(double d[MX], double z[MX][MX], long bands)