GRASS GIS 7 Programmer's Manual  7.0.5(2016)-r00000
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 += fabs(x[i]);
555  }
556 #pragma omp single
557  {
558  *value = s;
559  }
560  return;
561 }
562 
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 = fabs(x[rows - 1]);
583  for (i = rows - 2; i >= 0; i--) {
584  if (max < fabs(x[i]))
585  max = fabs(x[i]);
586  }
587 
588  *value = max;
589 }
590 
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 
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
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
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
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
#define max(a, b)
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