GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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  * */
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 }
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_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
int rows
Definition: bitmap.h:19
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
#define 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.
Definition: blas_level_2.c:115
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
double b
Definition: r_raster.c:39
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_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
int cols
Definition: bitmap.h:20
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