GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
blas_level_1.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 
25 #include <grass/gis.h>
26 #include <grass/gmath.h>
27 
28 /* **************************************************************** */
29 /* *************** D O U B L E ************************************ */
30 /* **************************************************************** */
31 
32 /*!
33  * \brief Compute the dot product of vector x and y
34  *
35  * \f[ a = {\bf x}^T {\bf y} \f]
36  *
37  * The functions creates its own parallel OpenMP region.
38  * It can be called within a parallel OpenMP region if nested parallelism is supported
39  * by the compiler.
40  *
41  * \param x (double *)
42  * \param y (double *)
43  * \param value (double *) -- the return value
44  * \param rows (int)
45  * \return (void)
46  *
47  * */
48 void G_math_d_x_dot_y(double *x, double *y, double *value, int rows)
49 {
50  int i;
51 
52  double s = 0.0;
53 
54 #pragma omp parallel for schedule (static) reduction(+:s)
55  for (i = rows - 1; i >= 0; i--) {
56  s += x[i] * y[i];
57  }
58 #pragma omp single
59  {
60  *value = s;
61  }
62  return;
63 }
64 
65 /*!
66  * \brief Compute the euclid norm of vector x
67  *
68  * \f[ a = ||{\bf x}||_2 \f]
69  *
70  * The functions creates its own parallel OpenMP region.
71  * It can be called within a parallel OpenMP region if nested parallelism is supported
72  * by the compiler.
73  *
74  * \param x (double *) -- the vector
75  * \param value (double *) -- the return value
76  * \param rows (int)
77  * \return (void)
78  *
79  * */
80 void G_math_d_euclid_norm(double *x, double *value, int rows)
81 {
82  int i;
83 
84  double s = 0.0;
85 
86 #pragma omp parallel for schedule (static) reduction(+:s)
87  for (i = rows - 1; i >= 0; i--) {
88  s += x[i] * x[i];
89  }
90 #pragma omp single
91  {
92  *value = sqrt(s);
93  }
94  return;
95 }
96 
97 /*!
98  * \brief Compute the asum norm of vector x
99  *
100  * \f[ a = ||{\bf x}||_1 \f]
101  *
102  * The functions creates its own parallel OpenMP region.
103  * It can be called within a parallel OpenMP region if nested parallelism is supported
104  * by the compiler.
105  *
106  * \param x (double *)-- the vector
107  * \param value (double *) -- the return value
108  * \param rows (int)
109  * \return (void)
110  *
111  * */
112 void G_math_d_asum_norm(double *x, double *value, int rows)
113 {
114  int i = 0;
115 
116  double s = 0.0;
117 
118 #pragma omp parallel for schedule (static) reduction(+:s)
119  for (i = rows - 1; i >= 0; i--) {
120  s += fabs(x[i]);
121  }
122 #pragma omp single
123  {
124  *value = s;
125  }
126  return;
127 }
128 
129 /*!
130  * \brief Compute the maximum norm of vector x
131  *
132  * \f[ a = ||{\bf x}||_\infty \f]
133  *
134  * This function is not multi-threaded
135  *
136  * \param x (double *)-- the vector
137  * \param value (double *) -- the return value
138  * \param rows (int)
139  * \return (void)
140  *
141  * */
142 void G_math_d_max_norm(double *x, double *value, int rows)
143 {
144  int i;
145 
146  double max = 0.0;
147 
148  max = fabs(x[rows - 1]);
149  for (i = rows - 2; i >= 0; i--) {
150  if (max < fabs(x[i]))
151  max = fabs(x[i]);
152  }
153 
154  *value = max;
155 }
156 
157 /*!
158  * \brief Scales vectors x and y with the scalars a and b and adds them
159  *
160  * \f[ {\bf z} = a{\bf x} + b{\bf y} \f]
161  *
162  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
163  *
164  * \param x (double *)
165  * \param y (double *)
166  * \param z (double *)
167  * \param a (double)
168  * \param b (double)
169  * \param rows (int)
170  * \return (void)
171  *
172  * */
173 void G_math_d_ax_by(double *x, double *y, double *z, double a, double b,
174  int rows)
175 {
176  int i;
177 
178  /*find specific cases */
179  if (b == 0.0) {
180 #pragma omp for schedule (static)
181  for (i = rows - 1; i >= 0; i--) {
182  z[i] = a * x[i];
183  }
184  }
185  else if ((a == 1.0) && (b == 1.0)) {
186 #pragma omp for schedule (static)
187  for (i = rows - 1; i >= 0; i--) {
188  z[i] = x[i] + y[i];
189  }
190  }
191  else if ((a == 1.0) && (b == -1.0)) {
192 #pragma omp for schedule (static)
193  for (i = rows - 1; i >= 0; i--) {
194  z[i] = x[i] - y[i];
195  }
196  }
197  else if (a == b) {
198 #pragma omp for schedule (static)
199  for (i = rows - 1; i >= 0; i--) {
200  z[i] = a * (x[i] + y[i]);
201  }
202  }
203  else if (b == -1.0) {
204 #pragma omp for schedule (static)
205  for (i = rows - 1; i >= 0; i--) {
206  z[i] = a * x[i] - y[i];
207  }
208  }
209  else if (b == 1.0) {
210 #pragma omp for schedule (static)
211  for (i = rows - 1; i >= 0; i--) {
212  z[i] = a * x[i] + y[i];
213  }
214  }
215  else {
216 #pragma omp for schedule (static)
217  for (i = rows - 1; i >= 0; i--) {
218  z[i] = a * x[i] + b * y[i];
219  }
220  }
221 
222  return;
223 }
224 
225 /*!
226  * \brief Copy the vector x to y
227  *
228  * \f[ {\bf y} = {\bf x} \f]
229  *
230  * This function is not multi-threaded
231  *
232  * \param x (double *)
233  * \param y (double *)
234  * \param rows (int)
235  *
236  * */
237 void G_math_d_copy(double *x, double *y, int rows)
238 {
239  y = memcpy(y, x, rows * sizeof(double));
240 
241  return;
242 }
243 
244 /* **************************************************************** */
245 /* *************** F L O A T ************************************** */
246 /* **************************************************************** */
247 
248 /*!
249  * \brief Compute the dot product of vector x and y
250  *
251  * \f[ a = {\bf x}^T {\bf y} \f]
252  *
253  * The functions creates its own parallel OpenMP region.
254  * It can be called within a parallel OpenMP region if nested parallelism is supported
255  * by the compiler.
256  *
257  * \param x (float *)
258  * \param y (float *)
259  * \param value (float *) -- the return value
260  * \param rows (int)
261  * \return (void)
262  *
263  * */
264 void G_math_f_x_dot_y(float *x, float *y, float *value, int rows)
265 {
266  int i;
267 
268  float s = 0.0;
269 
270 #pragma omp parallel for schedule (static) reduction(+:s)
271  for (i = rows - 1; i >= 0; i--) {
272  s += x[i] * y[i];
273  }
274 #pragma omp single
275  {
276  *value = s;
277  }
278  return;
279 }
280 
281 /*!
282  * \brief Compute the euclid norm of vector x
283  *
284  * \f[ a = ||{\bf x}||_2 \f]
285  *
286  * The functions creates its own parallel OpenMP region.
287  * It can be called within a parallel OpenMP region if nested parallelism is supported
288  * by the compiler.
289  *
290  * \param x (double *) -- the vector
291  * \param value (float *) -- the return value
292  * \param rows (int)
293  * \return (void)
294  *
295  * */
296 void G_math_f_euclid_norm(float *x, float *value, int rows)
297 {
298  int i;
299 
300  float s = 0.0;
301 
302 #pragma omp parallel for schedule (static) reduction(+:s)
303  for (i = rows - 1; i >= 0; i--) {
304  s += x[i] * x[i];
305  }
306 #pragma omp single
307  {
308  *value = sqrt(s);
309  }
310  return;
311 }
312 
313 /*!
314  * \brief Compute the asum norm of vector x
315  *
316  * \f[ a = ||{\bf x}||_1 \f]
317  *
318  * The functions creates its own parallel OpenMP region.
319  * It can be called within a parallel OpenMP region if nested parallelism is supported
320  * by the compiler.
321  *
322  * \param x (float *)-- the vector
323  * \param value (float *) -- the return value
324  * \param rows (int)
325  * \return (void)
326  *
327  * */
328 void G_math_f_asum_norm(float *x, float *value, int rows)
329 {
330  int i;
331 
332  int count = 0;
333 
334  float s = 0.0;
335 
336 #pragma omp parallel for schedule (static) private(i) reduction(+:s, count)
337  for (i = 0; i < rows; i++) {
338  s += fabs(x[i]);
339  count++;
340  }
341 #pragma omp single
342  {
343  *value = s;
344  }
345  return;
346 }
347 
348 /*!
349  * \brief Compute the maximum norm of vector x
350  *
351  * \f[ a = ||{\bf x}||_\infty \f]
352  *
353  * This function is not multi-threaded
354  *
355  * \param x (float *)-- the vector
356  * \param value (float *) -- the return value
357  * \param rows (int)
358  * \return (void)
359  *
360  * */
361 void G_math_f_max_norm(float *x, float *value, int rows)
362 {
363  int i;
364 
365  float max = 0.0;
366 
367  max = fabs(x[rows - 1]);
368  for (i = rows - 2; i >= 0; i--) {
369  if (max < fabs(x[i]))
370  max = fabs(x[i]);
371  }
372  *value = max;
373  return;
374 }
375 
376 /*!
377  * \brief Scales vectors x and y with the scalars a and b and adds them
378  *
379  * \f[ {\bf z} = a{\bf x} + b{\bf y} \f]
380  *
381  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
382  *
383  * \param x (float *)
384  * \param y (float *)
385  * \param z (float *)
386  * \param a (float)
387  * \param b (float)
388  * \param rows (int)
389  * \return (void)
390  *
391  * */
392 void G_math_f_ax_by(float *x, float *y, float *z, float a, float b, int rows)
393 {
394  int i;
395 
396  /*find specific cases */
397  if (b == 0.0) {
398 #pragma omp for schedule (static)
399  for (i = rows - 1; i >= 0; i--) {
400  z[i] = a * x[i];
401  }
402  }
403  else if ((a == 1.0) && (b == 1.0)) {
404 #pragma omp for schedule (static)
405  for (i = rows - 1; i >= 0; i--) {
406  z[i] = x[i] + y[i];
407  }
408  }
409  else if ((a == 1.0) && (b == -1.0)) {
410 #pragma omp for schedule (static)
411  for (i = rows - 1; i >= 0; i--) {
412  z[i] = x[i] - y[i];
413  }
414  }
415  else if (a == b) {
416 #pragma omp for schedule (static)
417  for (i = rows - 1; i >= 0; i--) {
418  z[i] = a * (x[i] + y[i]);
419  }
420  }
421  else if (b == -1.0) {
422 #pragma omp for schedule (static)
423  for (i = rows - 1; i >= 0; i--) {
424  z[i] = a * x[i] - y[i];
425  }
426  }
427  else if (b == 1.0) {
428 #pragma omp for schedule (static)
429  for (i = rows - 1; i >= 0; i--) {
430  z[i] = a * x[i] + y[i];
431  }
432  }
433  else {
434 #pragma omp for schedule (static)
435  for (i = rows - 1; i >= 0; i--) {
436  z[i] = a * x[i] + b * y[i];
437  }
438  }
439 
440  return;
441 }
442 
443 /*!
444  * \brief Copy the vector x to y
445  *
446  * \f[ {\bf y} = {\bf x} \f]
447  *
448  * This function is not multi-threaded
449  *
450  * \param x (float *)
451  * \param y (float *)
452  * \param rows (int)
453  *
454  * */
455 void G_math_f_copy(float *x, float *y, int rows)
456 {
457  y = memcpy(y, x, rows * sizeof(float));
458 
459  return;
460 }
461 
462 /* **************************************************************** */
463 /* *************** I N T E G E R ********************************** */
464 /* **************************************************************** */
465 
466 /*!
467  * \brief Compute the dot product of vector x and y
468  *
469  * \f[ a = {\bf x}^T {\bf y} \f]
470  *
471  * The functions creates its own parallel OpenMP region.
472  * It can be called within a parallel OpenMP region if nested parallelism is supported
473  * by the compiler.
474  *
475  * \param x (int *)
476  * \param y (int *)
477  * \param value (double *) -- the return value
478  * \param rows (int)
479  * \return (void)
480  *
481  * */
482 void G_math_i_x_dot_y(int *x, int *y, double *value, int rows)
483 {
484  int i;
485 
486  double s = 0.0;
487 
488 #pragma omp parallel for schedule (static) reduction(+:s)
489  for (i = rows - 1; i >= 0; i--) {
490  s += x[i] * y[i];
491  }
492 #pragma omp single
493  {
494  *value = s;
495  }
496  return;
497 }
498 
499 /*!
500  * \brief Compute the euclid norm of vector x
501  *
502  * \f[ a = ||{\bf x}||_2 \f]
503  *
504  * The functions creates its own parallel OpenMP region.
505  * It can be called within a parallel OpenMP region if nested parallelism is supported
506  * by the compiler.
507  *
508  * \param x (int *) -- the vector
509  * \param value (double *) -- the return value
510  * \param rows (int)
511  * \return (void)
512  *
513  * */
514 void G_math_i_euclid_norm(int *x, double *value, int rows)
515 {
516  int i;
517 
518  double s = 0.0;
519 
520 #pragma omp parallel for schedule (static) reduction(+:s)
521  for (i = rows - 1; i >= 0; i--) {
522  s += x[i] * x[i];
523  }
524 #pragma omp single
525  {
526  *value = sqrt(s);
527  }
528  return;
529 }
530 
531 /*!
532  * \brief Compute the asum norm of vector x
533  *
534  * \f[ a = ||{\bf x}||_1 \f]
535  *
536  * The functions creates its own parallel OpenMP region.
537  * It can be called within a parallel OpenMP region if nested parallelism is supported
538  * by the compiler.
539  *
540  * \param x (int *)-- the vector
541  * \param value (double *) -- the return value
542  * \param rows (int)
543  * \return (void)
544  *
545  * */
546 void G_math_i_asum_norm(int *x, double *value, int rows)
547 {
548  int i;
549 
550  double s = 0.0;
551 
552 #pragma omp parallel for schedule (static) reduction(+:s)
553  for (i = rows - 1; i >= 0; i--) {
554  s += (double)abs(x[i]);
555  }
556 #pragma omp single
557  {
558  *value = s;
559  }
560  return;
561 }
562 
563 /*!
564  * \brief Compute the maximum norm of vector x
565  *
566  * \f[ a = ||{\bf x}||_\infty \f]
567  *
568  * This function is not multi-threaded
569  *
570  * \param x (int *)-- the vector
571  * \param value (int *) -- the return value
572  * \param rows (int)
573  * \return (void)
574  *
575  * */
576 void G_math_i_max_norm(int *x, int *value, int rows)
577 {
578  int i;
579 
580  int max = 0.0;
581 
582  max = abs(x[rows - 1]);
583  for (i = rows - 2; i >= 0; i--) {
584  if (max < abs(x[i]))
585  max = abs(x[i]);
586  }
587 
588  *value = max;
589 }
590 
591 /*!
592  * \brief Scales vectors x and y with the scalars a and b and adds them
593  *
594  * \f[ {\bf z} = a{\bf x} + b{\bf y} \f]
595  *
596  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
597  *
598  * \param x (int *)
599  * \param y (int *)
600  * \param z (int *)
601  * \param a (int)
602  * \param b (int)
603  * \param rows (int)
604  * \return (void)
605  *
606  * */
607 void G_math_i_ax_by(int *x, int *y, int *z, int a, int b, int rows)
608 {
609  int i;
610 
611  /*find specific cases */
612  if (b == 0.0) {
613 #pragma omp for schedule (static)
614  for (i = rows - 1; i >= 0; i--) {
615  z[i] = a * x[i];
616  }
617  }
618  else if ((a == 1.0) && (b == 1.0)) {
619 #pragma omp for schedule (static)
620  for (i = rows - 1; i >= 0; i--) {
621  z[i] = x[i] + y[i];
622  }
623  }
624  else if ((a == 1.0) && (b == -1.0)) {
625 #pragma omp for schedule (static)
626  for (i = rows - 1; i >= 0; i--) {
627  z[i] = x[i] - y[i];
628  }
629  }
630  else if (a == b) {
631 #pragma omp for schedule (static)
632  for (i = rows - 1; i >= 0; i--) {
633  z[i] = a * (x[i] + y[i]);
634  }
635  }
636  else if (b == -1.0) {
637 #pragma omp for schedule (static)
638  for (i = rows - 1; i >= 0; i--) {
639  z[i] = a * x[i] - y[i];
640  }
641  }
642  else if (b == 1.0) {
643 #pragma omp for schedule (static)
644  for (i = rows - 1; i >= 0; i--) {
645  z[i] = a * x[i] + y[i];
646  }
647  }
648  else {
649 #pragma omp for schedule (static)
650  for (i = rows - 1; i >= 0; i--) {
651  z[i] = a * x[i] + b * y[i];
652  }
653  }
654 
655  return;
656 }
657 
658 /*!
659  * \brief Copy the vector x to y
660  *
661  * \f[ {\bf y} = {\bf x} \f]
662  *
663  * This function is not multi-threaded
664  *
665  * \param x (int *)
666  * \param y (int *)
667  * \param rows (int)
668  *
669  * */
670 void G_math_i_copy(int *x, int *y, int rows)
671 {
672  y = memcpy(y, x, rows * sizeof(int));
673 
674  return;
675 }
void G_math_f_x_dot_y(float *x, float *y, float *value, int rows)
Compute the dot product of vector x and y.
Definition: blas_level_1.c:264
void G_math_i_max_norm(int *x, int *value, int rows)
Compute the maximum norm of vector x.
Definition: blas_level_1.c:576
void G_math_f_ax_by(float *x, float *y, float *z, float a, float b, int rows)
Scales vectors x and y with the scalars a and b and adds them.
Definition: blas_level_1.c:392
void G_math_d_max_norm(double *x, double *value, int rows)
Compute the maximum norm of vector x.
Definition: blas_level_1.c:142
void G_math_d_x_dot_y(double *x, double *y, double *value, int rows)
Compute the dot product of vector x and y.
Definition: blas_level_1.c:48
int count
int rows
Definition: bitmap.h:19
void G_math_f_copy(float *x, float *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:455
#define x
#define max(x, y)
Definition: draw2.c:32
void G_math_f_max_norm(float *x, float *value, int rows)
Compute the maximum norm of vector x.
Definition: blas_level_1.c:361
void G_math_f_asum_norm(float *x, float *value, int rows)
Compute the asum norm of vector x.
Definition: blas_level_1.c:328
double b
Definition: r_raster.c:39
void G_math_i_copy(int *x, int *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:670
void G_math_d_euclid_norm(double *x, double *value, int rows)
Compute the euclid norm of vector x.
Definition: blas_level_1.c:80
void G_math_d_copy(double *x, double *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:237
void G_math_d_asum_norm(double *x, double *value, int rows)
Compute the asum norm of vector x.
Definition: blas_level_1.c:112
void G_math_i_x_dot_y(int *x, int *y, double *value, int rows)
Compute the dot product of vector x and y.
Definition: blas_level_1.c:482
void G_math_i_ax_by(int *x, int *y, int *z, int a, int b, int rows)
Scales vectors x and y with the scalars a and b and adds them.
Definition: blas_level_1.c:607
void G_math_f_euclid_norm(float *x, float *value, int rows)
Compute the euclid norm of vector x.
Definition: blas_level_1.c:296
void G_math_i_asum_norm(int *x, double *value, int rows)
Compute the asum norm of vector x.
Definition: blas_level_1.c:546
void G_math_d_ax_by(double *x, double *y, double *z, double a, double b, int rows)
Scales vectors x and y with the scalars a and b and adds them.
Definition: blas_level_1.c:173
void G_math_i_euclid_norm(int *x, double *value, int rows)
Compute the euclid norm of vector x.
Definition: blas_level_1.c:514