GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-6da4f45b8e
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  * */
46 void 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  * */
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
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  * */
116 void 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  * */
147 void 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 
182 void 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 
265 void 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 
342 int 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 
376 int 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.
Definition: blas_level_2.c:342
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:182
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:46
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:265
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:147
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:116
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:376
double b
Definition: r_raster.c:39
#define x