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