GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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 * */
48void 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 * */
80void 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 * */
112void 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 * */
142void 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 * */
173void 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 * */
237void 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 * */
264void 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 * */
296void 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 * */
328void 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 * */
361void 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 * */
392void 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 * */
455void 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 * */
482void 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 * */
514void 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 * */
546void 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 * */
576void 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 * */
607void 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 * */
670void 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_i_euclid_norm(int *x, double *value, int rows)
Compute the euclid norm of vector x
Definition: blas_level_1.c:514
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_f_asum_norm(float *x, float *value, int rows)
Compute the asum norm of vector x
Definition: blas_level_1.c:328
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_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_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_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
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_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_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_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_copy(double *x, double *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:237
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_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_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_i_copy(int *x, int *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:670
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_f_copy(float *x, float *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:455
double b
int count
#define max(a, b)
#define x