GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
blas_level_1.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
24#include <grass/gis.h>
25#include <grass/gmath.h>
26
27/* **************************************************************** */
28/* *************** D O U B L E ************************************ */
29/* **************************************************************** */
30
31/*!
32 * \brief Compute the dot product of vector x and y
33 *
34 * \f[ a = {\bf x}^T {\bf y} \f]
35 *
36 * The functions creates its own parallel OpenMP region.
37 * It can be called within a parallel OpenMP region if nested parallelism is
38 * supported by the compiler.
39 *
40 * \param x (double *)
41 * \param y (double *)
42 * \param value (double *) -- the return value
43 * \param rows (int)
44 * \return (void)
45 *
46 * */
47void G_math_d_x_dot_y(double *x, double *y, double *value, int rows)
48{
49 int i;
50
51 double s = 0.0;
52
53#pragma omp parallel for schedule(static) reduction(+ : s)
54 for (i = rows - 1; i >= 0; i--) {
55 s += x[i] * y[i];
56 }
57#pragma omp single
58 {
59 *value = s;
60 }
61 return;
62}
63
64/*!
65 * \brief Compute the euclid norm of vector x
66 *
67 * \f[ a = ||{\bf x}||_2 \f]
68 *
69 * The functions creates its own parallel OpenMP region.
70 * It can be called within a parallel OpenMP region if nested parallelism is
71 * supported by the compiler.
72 *
73 * \param x (double *) -- the vector
74 * \param value (double *) -- the return value
75 * \param rows (int)
76 * \return (void)
77 *
78 * */
79void G_math_d_euclid_norm(double *x, double *value, int rows)
80{
81 int i;
82
83 double s = 0.0;
84
85#pragma omp parallel for schedule(static) reduction(+ : s)
86 for (i = rows - 1; i >= 0; i--) {
87 s += x[i] * x[i];
88 }
89#pragma omp single
90 {
91 *value = sqrt(s);
92 }
93 return;
94}
95
96/*!
97 * \brief Compute the asum norm of vector x
98 *
99 * \f[ a = ||{\bf x}||_1 \f]
100 *
101 * The functions creates its own parallel OpenMP region.
102 * It can be called within a parallel OpenMP region if nested parallelism is
103 * supported by the compiler.
104 *
105 * \param x (double *)-- the vector
106 * \param value (double *) -- the return value
107 * \param rows (int)
108 * \return (void)
109 *
110 * */
111void G_math_d_asum_norm(double *x, double *value, int rows)
112{
113 int i = 0;
114
115 double s = 0.0;
116
117#pragma omp parallel for schedule(static) reduction(+ : s)
118 for (i = rows - 1; i >= 0; i--) {
119 s += fabs(x[i]);
120 }
121#pragma omp single
122 {
123 *value = s;
124 }
125 return;
126}
127
128/*!
129 * \brief Compute the maximum norm of vector x
130 *
131 * \f[ a = ||{\bf x}||_\infty \f]
132 *
133 * This function is not multi-threaded
134 *
135 * \param x (double *)-- the vector
136 * \param value (double *) -- the return value
137 * \param rows (int)
138 * \return (void)
139 *
140 * */
141void G_math_d_max_norm(double *x, double *value, int rows)
142{
143 int i;
144
145 double max = 0.0;
146
147 max = fabs(x[rows - 1]);
148 for (i = rows - 2; i >= 0; i--) {
149 if (max < fabs(x[i]))
150 max = fabs(x[i]);
151 }
152
153 *value = max;
154}
155
156/*!
157 * \brief Scales vectors x and y with the scalars a and b and adds them
158 *
159 * \f[ {\bf z} = a{\bf x} + b{\bf y} \f]
160 *
161 * This function is multi-threaded with OpenMP and can be called within a
162 * 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
255 * supported 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
288 * supported 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
320 * supported 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 float s = 0.0;
333
334#pragma omp parallel for schedule(static) private(i) reduction(+ : s)
335 for (i = 0; i < rows; i++) {
336 s += fabs(x[i]);
337 }
338#pragma omp single
339 {
340 *value = s;
341 }
342 return;
343}
344
345/*!
346 * \brief Compute the maximum norm of vector x
347 *
348 * \f[ a = ||{\bf x}||_\infty \f]
349 *
350 * This function is not multi-threaded
351 *
352 * \param x (float *)-- the vector
353 * \param value (float *) -- the return value
354 * \param rows (int)
355 * \return (void)
356 *
357 * */
358void G_math_f_max_norm(float *x, float *value, int rows)
359{
360 int i;
361
362 float max = 0.0;
363
364 max = fabs(x[rows - 1]);
365 for (i = rows - 2; i >= 0; i--) {
366 if (max < fabs(x[i]))
367 max = fabs(x[i]);
368 }
369 *value = max;
370 return;
371}
372
373/*!
374 * \brief Scales vectors x and y with the scalars a and b and adds them
375 *
376 * \f[ {\bf z} = a{\bf x} + b{\bf y} \f]
377 *
378 * This function is multi-threaded with OpenMP and can be called within a
379 * parallel OpenMP region.
380 *
381 * \param x (float *)
382 * \param y (float *)
383 * \param z (float *)
384 * \param a (float)
385 * \param b (float)
386 * \param rows (int)
387 * \return (void)
388 *
389 * */
390void G_math_f_ax_by(float *x, float *y, float *z, float a, float b, int rows)
391{
392 int i;
393
394 /*find specific cases */
395 if (b == 0.0) {
396#pragma omp for schedule(static)
397 for (i = rows - 1; i >= 0; i--) {
398 z[i] = a * x[i];
399 }
400 }
401 else if ((a == 1.0) && (b == 1.0)) {
402#pragma omp for schedule(static)
403 for (i = rows - 1; i >= 0; i--) {
404 z[i] = x[i] + y[i];
405 }
406 }
407 else if ((a == 1.0) && (b == -1.0)) {
408#pragma omp for schedule(static)
409 for (i = rows - 1; i >= 0; i--) {
410 z[i] = x[i] - y[i];
411 }
412 }
413 else if (a == b) {
414#pragma omp for schedule(static)
415 for (i = rows - 1; i >= 0; i--) {
416 z[i] = a * (x[i] + y[i]);
417 }
418 }
419 else if (b == -1.0) {
420#pragma omp for schedule(static)
421 for (i = rows - 1; i >= 0; i--) {
422 z[i] = a * x[i] - y[i];
423 }
424 }
425 else if (b == 1.0) {
426#pragma omp for schedule(static)
427 for (i = rows - 1; i >= 0; i--) {
428 z[i] = a * x[i] + y[i];
429 }
430 }
431 else {
432#pragma omp for schedule(static)
433 for (i = rows - 1; i >= 0; i--) {
434 z[i] = a * x[i] + b * y[i];
435 }
436 }
437
438 return;
439}
440
441/*!
442 * \brief Copy the vector x to y
443 *
444 * \f[ {\bf y} = {\bf x} \f]
445 *
446 * This function is not multi-threaded
447 *
448 * \param x (float *)
449 * \param y (float *)
450 * \param rows (int)
451 *
452 * */
453void G_math_f_copy(float *x, float *y, int rows)
454{
455 y = memcpy(y, x, rows * sizeof(float));
456
457 return;
458}
459
460/* **************************************************************** */
461/* *************** I N T E G E R ********************************** */
462/* **************************************************************** */
463
464/*!
465 * \brief Compute the dot product of vector x and y
466 *
467 * \f[ a = {\bf x}^T {\bf y} \f]
468 *
469 * The functions creates its own parallel OpenMP region.
470 * It can be called within a parallel OpenMP region if nested parallelism is
471 * supported by the compiler.
472 *
473 * \param x (int *)
474 * \param y (int *)
475 * \param value (double *) -- the return value
476 * \param rows (int)
477 * \return (void)
478 *
479 * */
480void G_math_i_x_dot_y(int *x, int *y, double *value, int rows)
481{
482 int i;
483
484 double s = 0.0;
485
486#pragma omp parallel for schedule(static) reduction(+ : s)
487 for (i = rows - 1; i >= 0; i--) {
488 s += x[i] * y[i];
489 }
490#pragma omp single
491 {
492 *value = s;
493 }
494 return;
495}
496
497/*!
498 * \brief Compute the euclid norm of vector x
499 *
500 * \f[ a = ||{\bf x}||_2 \f]
501 *
502 * The functions creates its own parallel OpenMP region.
503 * It can be called within a parallel OpenMP region if nested parallelism is
504 * supported by the compiler.
505 *
506 * \param x (int *) -- the vector
507 * \param value (double *) -- the return value
508 * \param rows (int)
509 * \return (void)
510 *
511 * */
512void G_math_i_euclid_norm(int *x, double *value, int rows)
513{
514 int i;
515
516 double s = 0.0;
517
518#pragma omp parallel for schedule(static) reduction(+ : s)
519 for (i = rows - 1; i >= 0; i--) {
520 s += x[i] * x[i];
521 }
522#pragma omp single
523 {
524 *value = sqrt(s);
525 }
526 return;
527}
528
529/*!
530 * \brief Compute the asum norm of vector x
531 *
532 * \f[ a = ||{\bf x}||_1 \f]
533 *
534 * The functions creates its own parallel OpenMP region.
535 * It can be called within a parallel OpenMP region if nested parallelism is
536 * supported by the compiler.
537 *
538 * \param x (int *)-- the vector
539 * \param value (double *) -- the return value
540 * \param rows (int)
541 * \return (void)
542 *
543 * */
544void G_math_i_asum_norm(int *x, double *value, int rows)
545{
546 int i;
547
548 double s = 0.0;
549
550#pragma omp parallel for schedule(static) reduction(+ : s)
551 for (i = rows - 1; i >= 0; i--) {
552 s += (double)abs(x[i]);
553 }
554#pragma omp single
555 {
556 *value = s;
557 }
558 return;
559}
560
561/*!
562 * \brief Compute the maximum norm of vector x
563 *
564 * \f[ a = ||{\bf x}||_\infty \f]
565 *
566 * This function is not multi-threaded
567 *
568 * \param x (int *)-- the vector
569 * \param value (int *) -- the return value
570 * \param rows (int)
571 * \return (void)
572 *
573 * */
574void G_math_i_max_norm(int *x, int *value, int rows)
575{
576 int i;
577
578 int max = 0.0;
579
580 max = abs(x[rows - 1]);
581 for (i = rows - 2; i >= 0; i--) {
582 if (max < abs(x[i]))
583 max = abs(x[i]);
584 }
585
586 *value = max;
587}
588
589/*!
590 * \brief Scales vectors x and y with the scalars a and b and adds them
591 *
592 * \f[ {\bf z} = a{\bf x} + b{\bf y} \f]
593 *
594 * This function is multi-threaded with OpenMP and can be called within a
595 * parallel OpenMP region.
596 *
597 * \param x (int *)
598 * \param y (int *)
599 * \param z (int *)
600 * \param a (int)
601 * \param b (int)
602 * \param rows (int)
603 * \return (void)
604 *
605 * */
606void G_math_i_ax_by(int *x, int *y, int *z, int a, int b, int rows)
607{
608 int i;
609
610 /*find specific cases */
611 if (b == 0.0) {
612#pragma omp for schedule(static)
613 for (i = rows - 1; i >= 0; i--) {
614 z[i] = a * x[i];
615 }
616 }
617 else if ((a == 1.0) && (b == 1.0)) {
618#pragma omp for schedule(static)
619 for (i = rows - 1; i >= 0; i--) {
620 z[i] = x[i] + y[i];
621 }
622 }
623 else if ((a == 1.0) && (b == -1.0)) {
624#pragma omp for schedule(static)
625 for (i = rows - 1; i >= 0; i--) {
626 z[i] = x[i] - y[i];
627 }
628 }
629 else if (a == b) {
630#pragma omp for schedule(static)
631 for (i = rows - 1; i >= 0; i--) {
632 z[i] = a * (x[i] + y[i]);
633 }
634 }
635 else if (b == -1.0) {
636#pragma omp for schedule(static)
637 for (i = rows - 1; i >= 0; i--) {
638 z[i] = a * x[i] - y[i];
639 }
640 }
641 else if (b == 1.0) {
642#pragma omp for schedule(static)
643 for (i = rows - 1; i >= 0; i--) {
644 z[i] = a * x[i] + y[i];
645 }
646 }
647 else {
648#pragma omp for schedule(static)
649 for (i = rows - 1; i >= 0; i--) {
650 z[i] = a * x[i] + b * y[i];
651 }
652 }
653
654 return;
655}
656
657/*!
658 * \brief Copy the vector x to y
659 *
660 * \f[ {\bf y} = {\bf x} \f]
661 *
662 * This function is not multi-threaded
663 *
664 * \param x (int *)
665 * \param y (int *)
666 * \param rows (int)
667 *
668 * */
669void G_math_i_copy(int *x, int *y, int rows)
670{
671 y = memcpy(y, x, rows * sizeof(int));
672
673 return;
674}
void G_math_i_euclid_norm(int *x, double *value, int rows)
Compute the euclid norm of vector x.
void G_math_d_asum_norm(double *x, double *value, int rows)
Compute the asum norm of vector x.
void G_math_f_asum_norm(float *x, float *value, int rows)
Compute the asum norm of vector x.
void G_math_d_max_norm(double *x, double *value, int rows)
Compute the maximum norm of vector x.
void G_math_i_asum_norm(int *x, double *value, int rows)
Compute the asum norm of vector x.
void G_math_f_max_norm(float *x, float *value, int rows)
Compute the maximum norm of vector x.
void G_math_d_x_dot_y(double *x, double *y, double *value, int rows)
Compute the dot product of vector x and y.
void G_math_f_euclid_norm(float *x, float *value, int rows)
Compute the euclid norm of vector x.
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.
void G_math_f_x_dot_y(float *x, float *y, float *value, int rows)
Compute the dot product of vector x and y.
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.
void G_math_d_copy(double *x, double *y, int rows)
Copy the vector x to y.
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.
void G_math_d_euclid_norm(double *x, double *value, int rows)
Compute the euclid norm of vector x.
void G_math_i_max_norm(int *x, int *value, int rows)
Compute the maximum norm of vector x.
void G_math_i_copy(int *x, int *y, int rows)
Copy the vector x to y.
void G_math_i_x_dot_y(int *x, int *y, double *value, int rows)
Compute the dot product of vector x and y.
void G_math_f_copy(float *x, float *y, int rows)
Copy the vector x to y.
double b
#define max(a, b)
#define x