← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to lib/gmath/blas_level_2.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gmath/blas_level_2.c
 
1
 
 
2
/*****************************************************************************
 
3
*
 
4
* MODULE:       Grass numerical math interface
 
5
* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
 
6
*               soerengebbert <at> googlemail <dot> com
 
7
*               
 
8
* PURPOSE:      grass blas implementation
 
9
*               part of the gmath library
 
10
*               
 
11
* COPYRIGHT:    (C) 2010 by the GRASS Development Team
 
12
*
 
13
*               This program is free software under the GNU General Public
 
14
*               License (>=v2). Read the file COPYING that comes with GRASS
 
15
*               for details.
 
16
*
 
17
*****************************************************************************/
 
18
 
 
19
#include <math.h>
 
20
#include <unistd.h>
 
21
#include <stdio.h>
 
22
#include <string.h>
 
23
#include <stdlib.h>
 
24
#include <grass/gmath.h>
 
25
#include <grass/gis.h>
 
26
 
 
27
#define EPSILON 0.00000000000000001
 
28
 
 
29
 
 
30
/*!
 
31
 * \brief Compute the matrix - vector product  
 
32
 * of matrix A and vector x.
 
33
 *
 
34
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
35
 *
 
36
 * y = A * x
 
37
 *
 
38
 *
 
39
 * \param A (double ** )
 
40
 * \param x (double *)
 
41
 * \param y (double *) 
 
42
 * \param rows (int)
 
43
 * \param cols (int)
 
44
 * \return (void)
 
45
 *
 
46
 * */
 
47
void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
 
48
{
 
49
    int i, j;
 
50
 
 
51
    double tmp;
 
52
 
 
53
#pragma omp for schedule (static) private(i, j, tmp)
 
54
    for (i = 0; i < rows; i++) {
 
55
        tmp = 0;
 
56
        for (j = cols - 1; j >= 0; j--) {
 
57
            tmp += A[i][j] * x[j];
 
58
        }
 
59
        y[i] = tmp;
 
60
    }
 
61
    return;
 
62
}
 
63
 
 
64
/*!
 
65
 * \brief Compute the matrix - vector product  
 
66
 * of matrix A and vector x.
 
67
 *
 
68
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
69
 *
 
70
 * y = A * x
 
71
 *
 
72
 *
 
73
 * \param A (float ** )
 
74
 * \param x (float *)
 
75
 * \param y (float *) 
 
76
 * \param rows (int)
 
77
 * \param cols (int)
 
78
 * \return (void)
 
79
 *
 
80
 * */
 
81
void G_math_f_Ax(float **A, float *x, float *y, int rows, int cols)
 
82
{
 
83
    int i, j;
 
84
 
 
85
    float tmp;
 
86
 
 
87
#pragma omp for schedule (static) private(i, j, tmp)
 
88
    for (i = 0; i < rows; i++) {
 
89
        tmp = 0;
 
90
        for (j = cols - 1; j >= 0; j--) {
 
91
            tmp += A[i][j] * x[j];
 
92
        }
 
93
        y[i] = tmp;
 
94
    }
 
95
    return;
 
96
}
 
97
 
 
98
/*!
 
99
 * \brief Compute the dyadic product of two vectors. 
 
100
 * The result is stored in the matrix A.
 
101
 *
 
102
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
103
 *
 
104
 * A = x * y^T
 
105
 *
 
106
 *
 
107
 * \param x (double *)
 
108
 * \param y (double *) 
 
109
 * \param A (float **)  -- matrix of size rows*cols
 
110
 * \param rows (int) -- length of vector x
 
111
 * \param cols (int) -- lengt of vector y
 
112
 * \return (void)
 
113
 *
 
114
 * */
 
115
void G_math_d_x_dyad_y(double *x, double *y, double **A, int rows, int cols)
 
116
{
 
117
    int i, j;
 
118
 
 
119
#pragma omp for schedule (static) private(i, j)
 
120
    for (i = 0; i < rows; i++) {
 
121
        for (j = cols - 1; j >= 0; j--) {
 
122
            A[i][j] = x[i] * y[j];
 
123
        }
 
124
    }
 
125
    return;
 
126
}
 
127
 
 
128
/*!
 
129
 * \brief Compute the dyadic product of two vectors. 
 
130
 * The result is stored in the matrix A.
 
131
 *
 
132
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
133
 *
 
134
 * A = x * y^T
 
135
 *
 
136
 *
 
137
 * \param x (float *)
 
138
 * \param y (float *) 
 
139
 * \param A (float **=  -- matrix of size rows*cols 
 
140
 * \param rows (int) -- length of vector x
 
141
 * \param cols (int) -- lengt of vector y
 
142
 * \return (void)
 
143
 *
 
144
 * */
 
145
void G_math_f_x_dyad_y(float *x, float *y, float **A, int rows, int cols)
 
146
{
 
147
    int i, j;
 
148
 
 
149
#pragma omp for schedule (static) private(i, j)
 
150
    for (i = 0; i < rows; i++) {
 
151
        for (j = cols - 1; j >= 0; j--) {
 
152
            A[i][j] = x[i] * y[j];
 
153
        }
 
154
    }
 
155
    return;
 
156
}
 
157
 
 
158
/*!
 
159
 * \brief Compute the scaled matrix - vector product  
 
160
 * of matrix double **A and vector x and y.
 
161
 *
 
162
 * z = a * A * x + b * y
 
163
 *
 
164
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
165
 *
 
166
 *
 
167
 * \param A (double **) 
 
168
 * \param x (double *)
 
169
 * \param y (double *) 
 
170
 * \param a (double)
 
171
 * \param b (double)
 
172
 * \param z (double *) 
 
173
 * \param rows (int)
 
174
 * \param cols (int)
 
175
 * \return (void)
 
176
 *
 
177
 * */
 
178
 
 
179
void G_math_d_aAx_by(double **A, double *x, double *y, double a, double b,
 
180
                     double *z, int rows, int cols)
 
181
{
 
182
    int i, j;
 
183
 
 
184
    double tmp;
 
185
 
 
186
    /*catch specific cases */
 
187
    if (a == b) {
 
188
#pragma omp for schedule (static) private(i, j, tmp)
 
189
        for (i = 0; i < rows; i++) {
 
190
            tmp = 0;
 
191
            for (j = cols - 1; j >= 0; j--) {
 
192
                tmp += A[i][j] * x[j] + y[j];
 
193
            }
 
194
            z[i] = a * tmp;
 
195
        }
 
196
    }
 
197
    else if (b == -1.0) {
 
198
#pragma omp for schedule (static) private(i, j, tmp)
 
199
        for (i = 0; i < rows; i++) {
 
200
            tmp = 0;
 
201
            for (j = cols - 1; j >= 0; j--) {
 
202
                tmp += a * A[i][j] * x[j] - y[j];
 
203
            }
 
204
            z[i] = tmp;
 
205
        }
 
206
    }
 
207
    else if (b == 0.0) {
 
208
#pragma omp for schedule (static) private(i, j, tmp)
 
209
        for (i = 0; i < rows; i++) {
 
210
            tmp = 0;
 
211
            for (j = cols - 1; j >= 0; j--) {
 
212
                tmp += A[i][j] * x[j];
 
213
            }
 
214
            z[i] = a * tmp;
 
215
        }
 
216
    }
 
217
    else if (a == -1.0) {
 
218
#pragma omp for schedule (static) private(i, j, tmp)
 
219
        for (i = 0; i < rows; i++) {
 
220
            tmp = 0;
 
221
            for (j = cols - 1; j >= 0; j--) {
 
222
                tmp += b * y[j] - A[i][j] * x[j];
 
223
            }
 
224
            z[i] = tmp;
 
225
        }
 
226
    }
 
227
    else {
 
228
#pragma omp for schedule (static) private(i, j, tmp)
 
229
        for (i = 0; i < rows; i++) {
 
230
            tmp = 0;
 
231
            for (j = cols - 1; j >= 0; j--) {
 
232
                tmp += a * A[i][j] * x[j] + b * y[j];
 
233
            }
 
234
            z[i] = tmp;
 
235
        }
 
236
    }
 
237
    return;
 
238
}
 
239
 
 
240
/*!
 
241
 * \brief Compute the scaled matrix - vector product  
 
242
 * of matrix A and vectors x and y.
 
243
 *
 
244
 * z = a * A * x + b * y
 
245
 *
 
246
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
247
 *
 
248
 *
 
249
 * \param A (float **) 
 
250
 * \param x (float *)
 
251
 * \param y (float *) 
 
252
 * \param a (float)
 
253
 * \param b (float)
 
254
 * \param z (float *) 
 
255
 * \param rows (int)
 
256
 * \param cols (int)
 
257
 * \return (void)
 
258
 *
 
259
 * */
 
260
 
 
261
void G_math_f_aAx_by(float **A, float *x, float *y, float a, float b,
 
262
                     float *z, int rows, int cols)
 
263
{
 
264
    int i, j;
 
265
 
 
266
    float tmp;
 
267
 
 
268
    /*catch specific cases */
 
269
    if (a == b) {
 
270
#pragma omp for schedule (static) private(i, j, tmp)
 
271
        for (i = 0; i < rows; i++) {
 
272
            tmp = 0;
 
273
            for (j = cols - 1; j >= 0; j--) {
 
274
                tmp += A[i][j] * x[j] + y[j];
 
275
            }
 
276
            z[i] = a * tmp;
 
277
        }
 
278
    }
 
279
    else if (b == -1.0) {
 
280
#pragma omp for schedule (static) private(i, j, tmp)
 
281
        for (i = 0; i < rows; i++) {
 
282
            tmp = 0;
 
283
            for (j = cols - 1; j >= 0; j--) {
 
284
                tmp += a * A[i][j] * x[j] - y[j];
 
285
            }
 
286
            z[i] = tmp;
 
287
        }
 
288
    }
 
289
    else if (b == 0.0) {
 
290
#pragma omp for schedule (static) private(i, j, tmp)
 
291
        for (i = 0; i < rows; i++) {
 
292
            tmp = 0;
 
293
            for (j = cols - 1; j >= 0; j--) {
 
294
                tmp += A[i][j] * x[j];
 
295
            }
 
296
            z[i] = a * tmp;
 
297
        }
 
298
    }
 
299
    else if (a == -1.0) {
 
300
#pragma omp for schedule (static) private(i, j, tmp)
 
301
        for (i = 0; i < rows; i++) {
 
302
            tmp = 0;
 
303
            for (j = cols - 1; j >= 0; j--) {
 
304
                tmp += b * y[j] - A[i][j] * x[j];
 
305
            }
 
306
            z[i] = tmp;
 
307
        }
 
308
    }
 
309
    else {
 
310
#pragma omp for schedule (static) private(i, j, tmp)
 
311
        for (i = 0; i < rows; i++) {
 
312
            tmp = 0;
 
313
            for (j = cols - 1; j >= 0; j--) {
 
314
                tmp += a * A[i][j] * x[j] + b * y[j];
 
315
            }
 
316
            z[i] = tmp;
 
317
        }
 
318
    }
 
319
    return;
 
320
}
 
321
 
 
322
 
 
323
 
 
324
/*!
 
325
 * \fn int G_math_d_A_T(double **A, int rows)
 
326
 *
 
327
 * \brief Compute the transposition of matrix A.
 
328
 * Matrix A will be overwritten.
 
329
 *
 
330
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
331
 *
 
332
 * Returns 0.
 
333
 *
 
334
 * \param A (double **)
 
335
 * \param rows (int)
 
336
 * \return int
 
337
 */
 
338
 
 
339
int G_math_d_A_T(double **A, int rows)
 
340
{
 
341
    int i, j;
 
342
 
 
343
    double tmp;
 
344
 
 
345
#pragma omp for schedule (static) private(i, j, tmp)
 
346
    for (i = 0; i < rows; i++)
 
347
        for (j = 0; j < i; j++) {
 
348
            tmp = A[i][j];
 
349
 
 
350
            A[i][j] = A[j][i];
 
351
            A[j][i] = tmp;
 
352
        }
 
353
 
 
354
    return 0;
 
355
}
 
356
 
 
357
/*!
 
358
 * \fn int G_math_f_A_T(float **A, int rows)
 
359
 *
 
360
 * \brief Compute the transposition of matrix A.
 
361
 * Matrix A will be overwritten.
 
362
 *
 
363
 * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
 
364
 *
 
365
 * Returns 0.
 
366
 *
 
367
 * \param A (float **)
 
368
 * \param rows (int)
 
369
 * \return int
 
370
 */
 
371
 
 
372
int G_math_f_A_T(float **A, int rows)
 
373
{
 
374
    int i, j;
 
375
 
 
376
    float tmp;
 
377
 
 
378
#pragma omp for schedule (static) private(i, j, tmp)
 
379
    for (i = 0; i < rows; i++)
 
380
        for (j = 0; j < i; j++) {
 
381
            tmp = A[i][j];
 
382
 
 
383
            A[i][j] = A[j][i];
 
384
            A[j][i] = tmp;
 
385
        }
 
386
 
 
387
    return 0;
 
388
}