← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/gmath/eigen.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gmath/eigen.c
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
 
 
11
 
/*!
12
 
 * \fn int eigen (double **M, double **Vectors, double *lambda, int n)
13
 
 *
14
 
 * \brief Computes eigenvalues (and eigen vectors if desired) for
15
 
 * symmetric matices.
16
 
 *
17
 
 * Computes eigenvalues (and eigen vectors if desired) for symmetric matices.
18
 
 *
19
 
 * \param M Input matrix
20
 
 * \param Vectors eigen output vector matrix
21
 
 * \param lambda Output eigenvalues
22
 
 * \param n Input matrix dimension
23
 
 * \return int
24
 
 */
25
 
 
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
 
 
58
 
/*!
59
 
 * \fn int egvorder2 (double *d, double **z, long bands)
60
 
 *
61
 
 * \brief
62
 
 *
63
 
 * Returns 0.
64
 
 *
65
 
 * \param d
66
 
 * \param z
67
 
 * \param bands
68
 
 * \return int
69
 
 */
70
 
 
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
 
 
108
 
/*!
109
 
 * \fn int transpose2 (double **eigmat, long bands)
110
 
 *
111
 
 * \brief
112
 
 *
113
 
 * Returns 0.
114
 
 *
115
 
 * \param eigmat
116
 
 * \param bands
117
 
 * \return int
118
 
 */
119
 
 
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
 
}