GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
blas_level_1.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3  *
4  * MODULE: Grass Gmath Library
5  * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2007
6  * soerengebbert <at> gmx <dot> de
7  *
8  * PURPOSE: blas level 1 like functions
9  * part of the gmath library
10  *
11  * COPYRIGHT: (C) 2000 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 /* **************************************************************** */
28 /* *************** D O U B L E ************************************ */
29 /* **************************************************************** */
30 
47 void 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 
79 void 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 
111 void 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 
141 void 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 
172 void G_math_d_ax_by(double *x, double *y, double *z, double a, double b,
173  int rows)
174 {
175  int i;
176 
177  /*find specific cases */
178  if (b == 0.0) {
179 #pragma omp for schedule (static)
180  for (i = rows - 1; i >= 0; i--) {
181  z[i] = a * x[i];
182  }
183  }
184  else if ((a == 1.0) && (b == 1.0)) {
185 #pragma omp for schedule (static)
186  for (i = rows - 1; i >= 0; i--) {
187  z[i] = x[i] + y[i];
188  }
189  }
190  else if ((a == 1.0) && (b == -1.0)) {
191 #pragma omp for schedule (static)
192  for (i = rows - 1; i >= 0; i--) {
193  z[i] = x[i] - y[i];
194  }
195  }
196  else if (a == b) {
197 #pragma omp for schedule (static)
198  for (i = rows - 1; i >= 0; i--) {
199  z[i] = a * (x[i] + y[i]);
200  }
201  }
202  else if (b == -1.0) {
203 #pragma omp for schedule (static)
204  for (i = rows - 1; i >= 0; i--) {
205  z[i] = a * x[i] - y[i];
206  }
207  }
208  else if (b == 1.0) {
209 #pragma omp for schedule (static)
210  for (i = rows - 1; i >= 0; i--) {
211  z[i] = a * x[i] + y[i];
212  }
213  }
214  else {
215 #pragma omp for schedule (static)
216  for (i = rows - 1; i >= 0; i--) {
217  z[i] = a * x[i] + b * y[i];
218  }
219  }
220 
221  return;
222 }
223 
236 void G_math_d_copy(double *x, double *y, int rows)
237 {
238  y = memcpy(y, x, rows * sizeof(double));
239 
240  return;
241 }
242 
243 /* **************************************************************** */
244 /* *************** F L O A T ************************************** */
245 /* **************************************************************** */
246 
263 void G_math_f_x_dot_y(float *x, float *y, float *value, int rows)
264 {
265  int i;
266 
267  float s = 0.0;
268 
269 #pragma omp parallel for schedule (static) reduction(+:s)
270  for (i = rows - 1; i >= 0; i--) {
271  s += x[i] * y[i];
272  }
273 #pragma omp single
274  {
275  *value = s;
276  }
277  return;
278 }
279 
295 void G_math_f_euclid_norm(float *x, float *value, int rows)
296 {
297  int i;
298 
299  float s = 0.0;
300 
301 #pragma omp parallel for schedule (static) reduction(+:s)
302  for (i = rows - 1; i >= 0; i--) {
303  s += x[i] * x[i];
304  }
305 #pragma omp single
306  {
307  *value = sqrt(s);
308  }
309  return;
310 }
311 
327 void G_math_f_asum_norm(float *x, float *value, int rows)
328 {
329  int i;
330 
331  int count = 0;
332 
333  float s = 0.0;
334 
335 #pragma omp parallel for schedule (static) private(i) reduction(+:s, count)
336  for (i = 0; i < rows; i++) {
337  s += fabs(x[i]);
338  count++;
339  }
340 #pragma omp single
341  {
342  *value = s;
343  }
344  return;
345 }
346 
360 void G_math_f_max_norm(float *x, float *value, int rows)
361 {
362  int i;
363 
364  float max = 0.0;
365 
366  max = fabs(x[rows - 1]);
367  for (i = rows - 2; i >= 0; i--) {
368  if (max < fabs(x[i]))
369  max = fabs(x[i]);
370  }
371  *value = max;
372  return;
373 }
374 
391 void G_math_f_ax_by(float *x, float *y, float *z, float a, float b, int rows)
392 {
393  int i;
394 
395  /*find specific cases */
396  if (b == 0.0) {
397 #pragma omp for schedule (static)
398  for (i = rows - 1; i >= 0; i--) {
399  z[i] = a * x[i];
400  }
401  }
402  else if ((a == 1.0) && (b == 1.0)) {
403 #pragma omp for schedule (static)
404  for (i = rows - 1; i >= 0; i--) {
405  z[i] = x[i] + y[i];
406  }
407  }
408  else if ((a == 1.0) && (b == -1.0)) {
409 #pragma omp for schedule (static)
410  for (i = rows - 1; i >= 0; i--) {
411  z[i] = x[i] - y[i];
412  }
413  }
414  else if (a == b) {
415 #pragma omp for schedule (static)
416  for (i = rows - 1; i >= 0; i--) {
417  z[i] = a * (x[i] + y[i]);
418  }
419  }
420  else if (b == -1.0) {
421 #pragma omp for schedule (static)
422  for (i = rows - 1; i >= 0; i--) {
423  z[i] = a * x[i] - y[i];
424  }
425  }
426  else if (b == 1.0) {
427 #pragma omp for schedule (static)
428  for (i = rows - 1; i >= 0; i--) {
429  z[i] = a * x[i] + y[i];
430  }
431  }
432  else {
433 #pragma omp for schedule (static)
434  for (i = rows - 1; i >= 0; i--) {
435  z[i] = a * x[i] + b * y[i];
436  }
437  }
438 
439  return;
440 }
441 
454 void G_math_f_copy(float *x, float *y, int rows)
455 {
456  y = memcpy(y, x, rows * sizeof(float));
457 
458  return;
459 }
460 
461 /* **************************************************************** */
462 /* *************** I N T E G E R ********************************** */
463 /* **************************************************************** */
464 
481 void G_math_i_x_dot_y(int *x, int *y, double *value, int rows)
482 {
483  int i;
484 
485  double s = 0.0;
486 
487 #pragma omp parallel for schedule (static) reduction(+:s)
488  for (i = rows - 1; i >= 0; i--) {
489  s += x[i] * y[i];
490  }
491 #pragma omp single
492  {
493  *value = s;
494  }
495  return;
496 }
497 
513 void G_math_i_euclid_norm(int *x, double *value, int rows)
514 {
515  int i;
516 
517  double s = 0.0;
518 
519 #pragma omp parallel for schedule (static) reduction(+:s)
520  for (i = rows - 1; i >= 0; i--) {
521  s += x[i] * x[i];
522  }
523 #pragma omp single
524  {
525  *value = sqrt(s);
526  }
527  return;
528 }
529 
545 void G_math_i_asum_norm(int *x, double *value, int rows)
546 {
547  int i;
548 
549  double s = 0.0;
550 
551 #pragma omp parallel for schedule (static) reduction(+:s)
552  for (i = rows - 1; i >= 0; i--) {
553  s += fabs(x[i]);
554  }
555 #pragma omp single
556  {
557  *value = s;
558  }
559  return;
560 }
561 
575 void G_math_i_max_norm(int *x, int *value, int rows)
576 {
577  int i;
578 
579  int max = 0.0;
580 
581  max = fabs(x[rows - 1]);
582  for (i = rows - 2; i >= 0; i--) {
583  if (max < fabs(x[i]))
584  max = fabs(x[i]);
585  }
586 
587  *value = max;
588 }
589 
606 void 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 
669 void 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_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:263
void G_math_i_max_norm(int *x, int *value, int rows)
Compute the maximum norm of vector x.
Definition: blas_level_1.c:575
float b
Definition: named_colr.c:8
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:391
void G_math_d_max_norm(double *x, double *value, int rows)
Compute the maximum norm of vector x.
Definition: blas_level_1.c:141
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:47
int count
void G_math_f_copy(float *x, float *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:454
int y
Definition: plot.c:34
#define max(x, y)
Definition: draw2.c:69
void G_math_f_max_norm(float *x, float *value, int rows)
Compute the maximum norm of vector x.
Definition: blas_level_1.c:360
void G_math_f_asum_norm(float *x, float *value, int rows)
Compute the asum norm of vector x.
Definition: blas_level_1.c:327
void G_math_i_copy(int *x, int *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:669
char * value
Definition: env.c:30
void G_math_d_euclid_norm(double *x, double *value, int rows)
Compute the euclid norm of vector x.
Definition: blas_level_1.c:79
void G_math_d_copy(double *x, double *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:236
void G_math_d_asum_norm(double *x, double *value, int rows)
Compute the asum norm of vector x.
Definition: blas_level_1.c:111
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:481
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:606
void G_math_f_euclid_norm(float *x, float *value, int rows)
Compute the euclid norm of vector x.
Definition: blas_level_1.c:295
void G_math_i_asum_norm(int *x, double *value, int rows)
Compute the asum norm of vector x.
Definition: blas_level_1.c:545
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:172
void G_math_i_euclid_norm(int *x, double *value, int rows)
Compute the euclid norm of vector x.
Definition: blas_level_1.c:513