← Back to branch summary

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

« back to all changes in this revision

Viewing changes to vector/v.generalize/matrix.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:
vector/v.generalize/matrix.c
22
22
#include <grass/glocale.h>
23
23
#include "matrix.h"
24
24
 
25
 
int matrix_init(int rows, int cols, MATRIX * res)
 
25
int matrix_init(int rows, int cols, MATRIX *res)
26
26
{
27
27
 
28
28
    int i, j;
46
46
    return 1;
47
47
}
48
48
 
49
 
void matrix_free(MATRIX m)
 
49
void matrix_free(MATRIX *m)
50
50
{
51
51
    int i;
52
52
 
53
 
    for (i = 0; i < m.rows; i++)
54
 
        G_free(m.a[i]);
55
 
    G_free(m.a);
 
53
    for (i = 0; i < m->rows; i++)
 
54
        G_free(m->a[i]);
 
55
    G_free(m->a);
56
56
    return;
57
57
}
58
58
 
59
 
int matrix_mult(MATRIX a, MATRIX b, MATRIX * res)
 
59
int matrix_mult(MATRIX *a, MATRIX *b, MATRIX *res)
60
60
{
61
 
    if (a.cols != b.rows)
 
61
    if (a->cols != b->rows)
62
62
        return 0;
63
63
 
64
64
    /*if (!matrix_init(a.rows, b.cols, res))
67
67
 
68
68
    int i, j, k;
69
69
 
70
 
    for (i = 0; i < a.rows; i++)
71
 
        for (j = 0; j < b.cols; j++) {
 
70
    for (i = 0; i < a->rows; i++)
 
71
        for (j = 0; j < b->cols; j++) {
72
72
            res->a[i][j] = 0;
73
 
            for (k = 0; k < a.cols; k++)
74
 
                res->a[i][j] += a.a[i][k] * b.a[k][j];
 
73
            for (k = 0; k < a->cols; k++)
 
74
                res->a[i][j] += a->a[i][k] * b->a[k][j];
75
75
        }
76
76
 
77
77
    return 1;
78
78
}
79
79
 
80
 
int matrix_add_identity(double s, MATRIX * m)
 
80
int matrix_add_identity(double s, MATRIX *m)
81
81
{
82
82
    if (m->rows != m->cols)
83
83
        return 0;
92
92
/* three following functions implements elementary row operations on matrices */
93
93
 
94
94
/* auxialiry function for matrix_inverse, swaps two rows of given matrix */
95
 
void matrix_swap_rows(int x, int y, MATRIX * m)
 
95
void matrix_swap_rows(int x, int y, MATRIX *m)
96
96
{
97
97
    int i;
98
98
 
108
108
 
109
109
/* auxiliary function for matrix_inverse, multiplies row of a matrix by
110
110
 * a scalar */
111
 
void matrix_row_scalar(int row, double s, MATRIX * m)
 
111
void matrix_row_scalar(int row, double s, MATRIX *m)
112
112
{
113
113
    int i;
114
114
 
121
121
 * one row to another.
122
122
 * i.e row[ra] = row[ra] + row[rb] * s;
123
123
 */
124
 
void matrix_row_add_multiple(int ra, int rb, double s, MATRIX * m)
 
124
void matrix_row_add_multiple(int ra, int rb, double s, MATRIX *m)
125
125
{
126
126
    int i;
127
127
 
131
131
}
132
132
 
133
133
/* TODO: don't test directly equality to zero */
134
 
int matrix_inverse(MATRIX a, MATRIX * res, int percents)
135
 
{;
 
134
int matrix_inverse(MATRIX *a, MATRIX *res, int percents)
 
135
{
 
136
    int i, j;
136
137
 
137
138
    /* not a square matrix */
138
 
    if (a.rows != a.cols)
 
139
    if (a->rows != a->cols)
139
140
        return 0;
140
141
 
141
 
    int i, j;
142
 
 
143
142
    /* initialize output matrix to the identity matrix */
144
 
    if (!matrix_init(a.rows, a.rows, res)) {
 
143
    if (!matrix_init(a->rows, a->rows, res)) {
145
144
        G_fatal_error(_("Out of memory"));
146
145
        return 0;
147
146
    }
148
 
    for (i = 0; i < a.rows; i++) {
149
 
        memset(res->a[i], 0, sizeof(double) * a.cols);
 
147
    for (i = 0; i < a->rows; i++) {
 
148
        memset(res->a[i], 0, sizeof(double) * a->cols);
150
149
        res->a[i][i] = 1;
151
150
    }
152
151
 
153
 
 
154
152
    /* in order to obtain the inverse of a matrix, we run
155
153
     * gauss elimination on the matrix and each time we apply
156
154
     * elementary row operation on this matrix, we apply the
159
157
     * is row equivalent to the identity matrix.
160
158
     */
161
159
 
162
 
    int n = a.rows;
 
160
    int n = a->rows;
163
161
 
164
162
    if (percents)
165
163
        G_percent_reset();
170
168
        if (percents)
171
169
            G_percent(i, n, 1);
172
170
        for (j = i; j < n; j++) {
173
 
            if (a.a[j][i] != 0) {       /* need to change this row to something */
 
171
            if (a->a[j][i] != 0) {      /* need to change this row to something */
174
172
                found = 1;      /* more sensible */
175
 
                matrix_swap_rows(i, j, &a);
 
173
                matrix_swap_rows(i, j, a);
176
174
                matrix_swap_rows(i, j, res);
177
175
                break;
178
176
            }
179
177
        }
180
178
        if (!found)
181
179
            return 0;
182
 
        double c = (double)1.0 / a.a[i][i];
 
180
        double c = (double)1.0 / a->a[i][i];
183
181
 
184
 
        matrix_row_scalar(i, c, &a);
 
182
        matrix_row_scalar(i, c, a);
185
183
        matrix_row_scalar(i, c, res);
186
184
        for (j = 0; j < n; j++) {
187
185
            if (i == j)
188
186
                continue;
189
 
            double c = -a.a[j][i];
 
187
            double c = -a->a[j][i];
190
188
 
191
189
            if (c == 0.0)
192
190
                continue;
193
 
            matrix_row_add_multiple(j, i, c, &a);
 
191
            matrix_row_add_multiple(j, i, c, a);
194
192
            matrix_row_add_multiple(j, i, c, res);
195
193
        }
196
194
    }
198
196
    return 1;
199
197
}
200
198
 
201
 
void matrix_mult_scalar(double s, MATRIX * m)
 
199
void matrix_mult_scalar(double s, MATRIX *m)
202
200
{
203
201
    int i, j;
204
202
 
207
205
            m->a[i][j] *= s;
208
206
}
209
207
 
210
 
void matrix_add(MATRIX a, MATRIX b, MATRIX * res)
 
208
void matrix_add(MATRIX *a, MATRIX *b, MATRIX *res)
211
209
{
212
210
    int i, j;
213
211
 
214
212
    for (i = 0; i < res->rows; i++)
215
213
        for (j = 0; j < res->cols; j++)
216
 
            res->a[i][j] = a.a[i][j] + b.a[i][j];
 
214
            res->a[i][j] = a->a[i][j] + b->a[i][j];
217
215
}
218
216
 
219
 
void matrix_print(MATRIX a)
 
217
void matrix_print(MATRIX *a)
220
218
{
221
219
    int i, j;
222
220
 
223
 
    for (i = 0; i < a.rows; i++) {
 
221
    for (i = 0; i < a->rows; i++) {
224
222
        double s = 0;
225
223
 
226
 
        for (j = 0; j < a.cols; j++) {
227
 
            printf("%.3lf ", a.a[i][j]);
228
 
            s += a.a[i][j];
 
224
        for (j = 0; j < a->cols; j++) {
 
225
            printf("%.3lf ", a->a[i][j]);
 
226
            s += a->a[i][j];
229
227
        }
230
228
        printf("|%.5lf\n", s);
231
229
    }