GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
blas_level_2.c
Go to the documentation of this file.
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 * */
47void 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 * */
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 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 * */
115void 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 * */
145void 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
179void 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
261void 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
339int 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
372int 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}
int G_math_d_A_T(double **A, int rows)
Compute the transposition of matrix A. Matrix A will be overwritten.
Definition: blas_level_2.c:339
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.
Definition: blas_level_2.c:179
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.
Definition: blas_level_2.c:47
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.
Definition: blas_level_2.c:261
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.
Definition: blas_level_2.c:145
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.
Definition: blas_level_2.c:81
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.
Definition: blas_level_2.c:115
int G_math_f_A_T(float **A, int rows)
Compute the transposition of matrix A. Matrix A will be overwritten.
Definition: blas_level_2.c:372
double b
#define x