GRASS GIS 7 Programmer's Manual  7.0.3(2016)-r00000
blas_level_2.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 #include <grass/gmath.h>
25 #include <grass/gis.h>
26 
27 #define EPSILON 0.00000000000000001
28 
29 
47 void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
48 {
49  int i, j;
50 
51  double tmp;
52 
53 #pragma omp for schedule (static) private(i, j, tmp)
54  for (i = 0; i < rows; i++) {
55  tmp = 0;
56  for (j = cols - 1; j >= 0; j--) {
57  tmp += A[i][j] * x[j];
58  }
59  y[i] = tmp;
60  }
61  return;
62 }
63 
81 void G_math_f_Ax(float **A, float *x, float *y, int rows, int cols)
82 {
83  int i, j;
84 
85  float tmp;
86 
87 #pragma omp for schedule (static) private(i, j, tmp)
88  for (i = 0; i < rows; i++) {
89  tmp = 0;
90  for (j = cols - 1; j >= 0; j--) {
91  tmp += A[i][j] * x[j];
92  }
93  y[i] = tmp;
94  }
95  return;
96 }
97 
115 void G_math_d_x_dyad_y(double *x, double *y, double **A, int rows, int cols)
116 {
117  int i, j;
118 
119 #pragma omp for schedule (static) private(i, j)
120  for (i = 0; i < rows; i++) {
121  for (j = cols - 1; j >= 0; j--) {
122  A[i][j] = x[i] * y[j];
123  }
124  }
125  return;
126 }
127 
145 void G_math_f_x_dyad_y(float *x, float *y, float **A, int rows, int cols)
146 {
147  int i, j;
148 
149 #pragma omp for schedule (static) private(i, j)
150  for (i = 0; i < rows; i++) {
151  for (j = cols - 1; j >= 0; j--) {
152  A[i][j] = x[i] * y[j];
153  }
154  }
155  return;
156 }
157 
179 void G_math_d_aAx_by(double **A, double *x, double *y, double a, double b,
180  double *z, int rows, int cols)
181 {
182  int i, j;
183 
184  double tmp;
185 
186  /*catch specific cases */
187  if (a == b) {
188 #pragma omp for schedule (static) private(i, j, tmp)
189  for (i = 0; i < rows; i++) {
190  tmp = 0;
191  for (j = cols - 1; j >= 0; j--) {
192  tmp += A[i][j] * x[j] + y[j];
193  }
194  z[i] = a * tmp;
195  }
196  }
197  else if (b == -1.0) {
198 #pragma omp for schedule (static) private(i, j, tmp)
199  for (i = 0; i < rows; i++) {
200  tmp = 0;
201  for (j = cols - 1; j >= 0; j--) {
202  tmp += a * A[i][j] * x[j] - y[j];
203  }
204  z[i] = tmp;
205  }
206  }
207  else if (b == 0.0) {
208 #pragma omp for schedule (static) private(i, j, tmp)
209  for (i = 0; i < rows; i++) {
210  tmp = 0;
211  for (j = cols - 1; j >= 0; j--) {
212  tmp += A[i][j] * x[j];
213  }
214  z[i] = a * tmp;
215  }
216  }
217  else if (a == -1.0) {
218 #pragma omp for schedule (static) private(i, j, tmp)
219  for (i = 0; i < rows; i++) {
220  tmp = 0;
221  for (j = cols - 1; j >= 0; j--) {
222  tmp += b * y[j] - A[i][j] * x[j];
223  }
224  z[i] = tmp;
225  }
226  }
227  else {
228 #pragma omp for schedule (static) private(i, j, tmp)
229  for (i = 0; i < rows; i++) {
230  tmp = 0;
231  for (j = cols - 1; j >= 0; j--) {
232  tmp += a * A[i][j] * x[j] + b * y[j];
233  }
234  z[i] = tmp;
235  }
236  }
237  return;
238 }
239 
261 void G_math_f_aAx_by(float **A, float *x, float *y, float a, float b,
262  float *z, int rows, int cols)
263 {
264  int i, j;
265 
266  float tmp;
267 
268  /*catch specific cases */
269  if (a == b) {
270 #pragma omp for schedule (static) private(i, j, tmp)
271  for (i = 0; i < rows; i++) {
272  tmp = 0;
273  for (j = cols - 1; j >= 0; j--) {
274  tmp += A[i][j] * x[j] + y[j];
275  }
276  z[i] = a * tmp;
277  }
278  }
279  else if (b == -1.0) {
280 #pragma omp for schedule (static) private(i, j, tmp)
281  for (i = 0; i < rows; i++) {
282  tmp = 0;
283  for (j = cols - 1; j >= 0; j--) {
284  tmp += a * A[i][j] * x[j] - y[j];
285  }
286  z[i] = tmp;
287  }
288  }
289  else if (b == 0.0) {
290 #pragma omp for schedule (static) private(i, j, tmp)
291  for (i = 0; i < rows; i++) {
292  tmp = 0;
293  for (j = cols - 1; j >= 0; j--) {
294  tmp += A[i][j] * x[j];
295  }
296  z[i] = a * tmp;
297  }
298  }
299  else if (a == -1.0) {
300 #pragma omp for schedule (static) private(i, j, tmp)
301  for (i = 0; i < rows; i++) {
302  tmp = 0;
303  for (j = cols - 1; j >= 0; j--) {
304  tmp += b * y[j] - A[i][j] * x[j];
305  }
306  z[i] = tmp;
307  }
308  }
309  else {
310 #pragma omp for schedule (static) private(i, j, tmp)
311  for (i = 0; i < rows; i++) {
312  tmp = 0;
313  for (j = cols - 1; j >= 0; j--) {
314  tmp += a * A[i][j] * x[j] + b * y[j];
315  }
316  z[i] = tmp;
317  }
318  }
319  return;
320 }
321 
322 
323 
339 int G_math_d_A_T(double **A, int rows)
340 {
341  int i, j;
342 
343  double tmp;
344 
345 #pragma omp for schedule (static) private(i, j, tmp)
346  for (i = 0; i < rows; i++)
347  for (j = 0; j < i; j++) {
348  tmp = A[i][j];
349 
350  A[i][j] = A[j][i];
351  A[j][i] = tmp;
352  }
353 
354  return 0;
355 }
356 
372 int G_math_f_A_T(float **A, int rows)
373 {
374  int i, j;
375 
376  float tmp;
377 
378 #pragma omp for schedule (static) private(i, j, tmp)
379  for (i = 0; i < rows; i++)
380  for (j = 0; j < i; j++) {
381  tmp = A[i][j];
382 
383  A[i][j] = A[j][i];
384  A[j][i] = tmp;
385  }
386 
387  return 0;
388 }
int G_math_d_A_T(double **A, int rows)
Compute the transposition of matrix A. Matrix A will be overwritten.
Definition: blas_level_2.c:339
void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
Compute the matrix - vector product of matrix A and vector x.
Definition: blas_level_2.c:47
int G_math_f_A_T(float **A, int rows)
Compute the transposition of matrix A. Matrix A will be overwritten.
Definition: blas_level_2.c:372
void G_math_d_x_dyad_y(double *x, double *y, double **A, int rows, int cols)
Compute the dyadic product of two vectors. The result is stored in the matrix A.
Definition: blas_level_2.c:115
double b
void G_math_d_aAx_by(double **A, double *x, double *y, double a, double b, double *z, int rows, int cols)
Compute the scaled matrix - vector product of matrix double **A and vector x and y.
Definition: blas_level_2.c:179
void G_math_f_Ax(float **A, float *x, float *y, int rows, int cols)
Compute the matrix - vector product of matrix A and vector x.
Definition: blas_level_2.c:81
void G_math_f_x_dyad_y(float *x, float *y, float **A, int rows, int cols)
Compute the dyadic product of two vectors. The result is stored in the matrix A.
Definition: blas_level_2.c:145
void G_math_f_aAx_by(float **A, float *x, float *y, float a, float b, float *z, int rows, int cols)
Compute the scaled matrix - vector product of matrix A and vectors x and y.
Definition: blas_level_2.c:261