GRASS Programmer's Manual  6.5.svn(2014)-r66266
blas_level_2.c
Go to the documentation of this file.
1
2 /*****************************************************************************
3  *
4  * MODULE: Grass PDE Numerical Library
5  * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2007
6  * soerengebbert <at> gmx <dot> de
7  *
8  * PURPOSE: linear equation system solvers
9  * part of the gpde library
10  *
11  * COPYRIGHT: (C) 2007 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 #include <grass/gisdefs.h>
27
28 #define EPSILON 0.00000000000000001
29
45 void G_math_Ax_sparse(G_math_spvector ** Asp, double *x, double *y, int rows)
46 {
47  int i, j;
48
49  double tmp;
50
51 #pragma omp for schedule (static) private(i, j, tmp)
52  for (i = 0; i < rows; i++) {
53  tmp = 0;
54  for (j = 0; j < Asp[i]->cols; j++) {
55  tmp += Asp[i]->values[j] * x[Asp[i]->index[j]];
56  }
57  y[i] = tmp;
58  }
59  return;
60 }
61
79 void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
80 {
81  int i, j;
82
83  double tmp;
84
85 #pragma omp for schedule (static) private(i, j, tmp)
86  for (i = 0; i < rows; i++) {
87  tmp = 0;
88  for (j = cols - 1; j >= 0; j--) {
89  tmp += A[i][j] * x[j];
90  }
91  y[i] = tmp;
92  }
93  return;
94 }
95
113 void G_math_f_Ax(float **A, float *x, float *y, int rows, int cols)
114 {
115  int i, j;
116
117  float tmp;
118
119 #pragma omp for schedule (static) private(i, j, tmp)
120  for (i = 0; i < rows; i++) {
121  tmp = 0;
122  for (j = cols - 1; j >= 0; j--) {
123  tmp += A[i][j] * x[j];
124  }
125  y[i] = tmp;
126  }
127  return;
128 }
129
147 void G_math_d_x_dyad_y(double *x, double *y, double **A, int rows, int cols)
148 {
149  int i, j;
150
151 #pragma omp for schedule (static) private(i, j)
152  for (i = 0; i < rows; i++) {
153  for (j = cols - 1; j >= 0; j--) {
154  A[i][j] = x[i] * y[j];
155  }
156  }
157  return;
158 }
159
177 void G_math_f_x_dyad_y(float *x, float *y, float **A, int rows, int cols)
178 {
179  int i, j;
180
181 #pragma omp for schedule (static) private(i, j)
182  for (i = 0; i < rows; i++) {
183  for (j = cols - 1; j >= 0; j--) {
184  A[i][j] = x[i] * y[j];
185  }
186  }
187  return;
188 }
189
211 void G_math_d_aAx_by(double **A, double *x, double *y, double a, double b,
212  double *z, int rows, int cols)
213 {
214  int i, j;
215
216  double tmp;
217
218  /*catch specific cases */
219  if (a == b) {
220 #pragma omp for schedule (static) private(i, j, tmp)
221  for (i = 0; i < rows; i++) {
222  tmp = 0;
223  for (j = cols - 1; j >= 0; j--) {
224  tmp += A[i][j] * x[j] + y[j];
225  }
226  z[i] = a * tmp;
227  }
228  }
229  else if (b == -1.0) {
230 #pragma omp for schedule (static) private(i, j, tmp)
231  for (i = 0; i < rows; i++) {
232  tmp = 0;
233  for (j = cols - 1; j >= 0; j--) {
234  tmp += a * A[i][j] * x[j] - y[j];
235  }
236  z[i] = tmp;
237  }
238  }
239  else if (b == 0.0) {
240 #pragma omp for schedule (static) private(i, j, tmp)
241  for (i = 0; i < rows; i++) {
242  tmp = 0;
243  for (j = cols - 1; j >= 0; j--) {
244  tmp += A[i][j] * x[j];
245  }
246  z[i] = a * tmp;
247  }
248  }
249  else if (a == -1.0) {
250 #pragma omp for schedule (static) private(i, j, tmp)
251  for (i = 0; i < rows; i++) {
252  tmp = 0;
253  for (j = cols - 1; j >= 0; j--) {
254  tmp += b * y[j] - A[i][j] * x[j];
255  }
256  z[i] = tmp;
257  }
258  }
259  else {
260 #pragma omp for schedule (static) private(i, j, tmp)
261  for (i = 0; i < rows; i++) {
262  tmp = 0;
263  for (j = cols - 1; j >= 0; j--) {
264  tmp += a * A[i][j] * x[j] + b * y[j];
265  }
266  z[i] = tmp;
267  }
268  }
269  return;
270 }
271
293 void G_math_f_aAx_by(float **A, float *x, float *y, float a, float b,
294  float *z, int rows, int cols)
295 {
296  int i, j;
297
298  float tmp;
299
300  /*catch specific cases */
301  if (a == b) {
302 #pragma omp for schedule (static) private(i, j, tmp)
303  for (i = 0; i < rows; i++) {
304  tmp = 0;
305  for (j = cols - 1; j >= 0; j--) {
306  tmp += A[i][j] * x[j] + y[j];
307  }
308  z[i] = a * tmp;
309  }
310  }
311  else if (b == -1.0) {
312 #pragma omp for schedule (static) private(i, j, tmp)
313  for (i = 0; i < rows; i++) {
314  tmp = 0;
315  for (j = cols - 1; j >= 0; j--) {
316  tmp += a * A[i][j] * x[j] - y[j];
317  }
318  z[i] = tmp;
319  }
320  }
321  else if (b == 0.0) {
322 #pragma omp for schedule (static) private(i, j, tmp)
323  for (i = 0; i < rows; i++) {
324  tmp = 0;
325  for (j = cols - 1; j >= 0; j--) {
326  tmp += A[i][j] * x[j];
327  }
328  z[i] = a * tmp;
329  }
330  }
331  else if (a == -1.0) {
332 #pragma omp for schedule (static) private(i, j, tmp)
333  for (i = 0; i < rows; i++) {
334  tmp = 0;
335  for (j = cols - 1; j >= 0; j--) {
336  tmp += b * y[j] - A[i][j] * x[j];
337  }
338  z[i] = tmp;
339  }
340  }
341  else {
342 #pragma omp for schedule (static) private(i, j, tmp)
343  for (i = 0; i < rows; i++) {
344  tmp = 0;
345  for (j = cols - 1; j >= 0; j--) {
346  tmp += a * A[i][j] * x[j] + b * y[j];
347  }
348  z[i] = tmp;
349  }
350  }
351  return;
352 }
353
354
355
371 int G_math_d_A_T(double **A, int rows)
372 {
373  int i, j;
374
375  double tmp;
376
377 #pragma omp for schedule (static) private(i, j, tmp)
378  for (i = 0; i < rows; i++)
379  for (j = 0; j < i; j++) {
380  tmp = A[i][j];
381
382  A[i][j] = A[j][i];
383  A[j][i] = tmp;
384  }
385
386  return 0;
387 }
388
404 int G_math_f_A_T(float **A, int rows)
405 {
406  int i, j;
407
408  float tmp;
409
410 #pragma omp for schedule (static) private(i, j, tmp)
411  for (i = 0; i < rows; i++)
412  for (j = 0; j < i; j++) {
413  tmp = A[i][j];
414
415  A[i][j] = A[j][i];
416  A[j][i] = tmp;
417  }
418
419  return 0;
420 }
float b
Definition: named_colr.c:8
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:371
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:79
int G_math_f_A_T(float **A, int rows)
Definition: blas_level_2.c:404
int y
Definition: plot.c:34
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:147
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:211
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:113
void G_math_f_x_dyad_y(float *x, float *y, float **A, int rows, int cols)
Compute the dyadic product of twMo vectors. The result is stored in the matrix A. ...
Definition: blas_level_2.c:177
tuple cols
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:293
void G_math_Ax_sparse(G_math_spvector **Asp, double *x, double *y, int rows)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
Definition: blas_level_2.c:45