GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
blas_level_3.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/gis.h>
25 #include <grass/gmath.h>
26 
27 
28 /*!
29  * \brief Add two matrices and scale matrix A with the scalar a
30  *
31  * \f[ {\bf C} = a {\bf A} + {\bf B} \f]
32  *
33  * In case B == NULL, matrix A will be scaled by scalar a. \n
34  * In case a == 1.0, a simple matrix addition is performed. \n
35  * In case a == -1.0 matrix A is subtracted from matrix B. \n
36  * The result is written into matrix C.
37  *
38  *
39  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
40  *
41  * \param A (double **)
42  * \param B (double **) if NULL, matrix A is scaled by scalar a only
43  * \param a (double)
44  * \param C (double **)
45  * \param rows (int)
46  * \param cols (int)
47  * \return (void)
48  *
49  * */
50 void G_math_d_aA_B(double **A, double **B, double a, double **C, int rows,
51  int cols)
52 {
53  int i, j;
54 
55 
56  /*If B is null, scale the matrix A with th scalar a */
57  if (B == NULL) {
58 #pragma omp for schedule (static) private(i, j)
59  for (i = rows - 1; i >= 0; i--)
60  for (j = cols - 1; j >= 0; j--)
61  C[i][j] = a * A[i][j];
62 
63  return;
64  }
65 
66  /*select special cases */
67  if (a == 1.0) {
68 #pragma omp for schedule (static) private(i, j)
69  for (i = rows - 1; i >= 0; i--)
70  for (j = cols - 1; j >= 0; j--)
71  C[i][j] = A[i][j] + B[i][j];
72  }
73  else if (a == -1.0) {
74 #pragma omp for schedule (static) private(i, j)
75  for (i = rows - 1; i >= 0; i--)
76  for (j = cols - 1; j >= 0; j--)
77  C[i][j] = B[i][j] - A[i][j];
78  }
79  else {
80 #pragma omp for schedule (static) private(i, j)
81  for (i = rows - 1; i >= 0; i--)
82  for (j = cols - 1; j >= 0; j--)
83  C[i][j] = a * A[i][j] + B[i][j];
84  }
85 
86  return;
87 }
88 
89 /*!
90  * \brief Add two matrices and scale matrix A with the scalar a
91  *
92  * \f[ {\bf C} = a {\bf A} + {\bf B} \f]
93  *
94  * In case B == NULL, matrix A will be scaled by scalar a. \n
95  * In case a == 1.0, a simple matrix addition is performed. \n
96  * In case a == -1.0 matrix A is subtracted from matrix B. \n
97  * The result is written into matrix C.
98  *
99  *
100  *
101  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
102  *
103  * \param A (float **)
104  * \param B (float **) if NULL, matrix A is scaled by scalar a only
105  * \param a (float)
106  * \param C (float **)
107  * \param rows (int)
108  * \param cols (int)
109 
110  * \return (void)
111  *
112  * */
113 void G_math_f_aA_B(float **A, float **B, float a, float **C, int rows,
114  int cols)
115 {
116  int i, j;
117 
118  /*If B is null, scale the matrix A with th scalar a */
119  if (B == NULL) {
120 #pragma omp for schedule (static) private(i, j)
121  for (i = rows - 1; i >= 0; i--)
122  for (j = cols - 1; j >= 0; j--)
123  C[i][j] = a * A[i][j];
124  return;
125  }
126 
127  /*select special cases */
128  if (a == 1.0) {
129 #pragma omp for schedule (static) private(i, j)
130  for (i = rows - 1; i >= 0; i--)
131  for (j = cols - 1; j >= 0; j--)
132  C[i][j] = A[i][j] + B[i][j];
133  }
134  else if (a == -1.0) {
135 #pragma omp for schedule (static) private(i, j)
136  for (i = rows - 1; i >= 0; i--)
137  for (j = cols - 1; j >= 0; j--)
138  C[i][j] = B[i][j] - A[i][j];
139  }
140  else {
141 #pragma omp for schedule (static) private(i, j)
142  for (i = rows - 1; i >= 0; i--)
143  for (j = cols - 1; j >= 0; j--)
144  C[i][j] = a * A[i][j] + B[i][j];
145  }
146 
147  return;
148 }
149 
150 
151 /*!
152  * \brief Matrix multiplication
153  *
154  * \f[ {\bf C} = {\bf A}{\bf B} \f]
155  *
156  * The result is written into matrix C.
157  *
158  * A must be of size rows_A * cols_A
159  * B must be of size rows_B * cols_B with rows_B == cols_A
160  * C must be of size rows_A * cols_B
161  *
162  *
163  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
164  *
165  * \param A (double **)
166  * \param B (double **)
167  * \param C (double **)
168  * \param rows_A (int)
169  * \param cols_A (int)
170  * \param cols_B (int)
171  * \return (void)
172  *
173  * */
174 void G_math_d_AB(double **A, double **B, double **C, int rows_A,
175  int cols_A, int cols_B)
176 {
177  int i, j, k;
178 
179 #pragma omp for schedule (static) private(i, j, k)
180  for (i = 0; i < rows_A; i++) {
181  for (j = 0; j < cols_B; j++) {
182  C[i][j] = 0.0;
183  for (k = cols_A - 1; k >= 0; k--) {
184  C[i][j] += A[i][k] * B[k][j];
185  }
186  }
187  }
188 
189  return;
190 }
191 
192 /*!
193  * \brief Matrix multiplication
194  *
195  * \f[ {\bf C} = {\bf A}{\bf B} \f]
196  *
197  * The result is written into matrix C.
198  *
199  * A must be of size rows_A * cols_A
200  * B must be of size rows_B * cols_B with rows_B == cols_A
201  * C must be of size rows_A * cols_B
202  *
203  *
204  * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
205  *
206  * \param A (float **)
207  * \param B (float **)
208  * \param C (float **)
209  * \param rows_A (int)
210  * \param cols_A (int)
211  * \param cols_B (int)
212  * \return (void)
213  *
214  * */
215 void G_math_f_AB(float **A, float **B, float **C, int rows_A,
216  int cols_A, int cols_B)
217 {
218  int i, j, k;
219 
220 #pragma omp for schedule (static) private(i, j, k)
221  for (i = 0; i < rows_A; i++) {
222  for (j = 0; j < cols_B; j++) {
223  C[i][j] = 0.0;
224  for (k = cols_A - 1; k >= 0; k--) {
225  C[i][j] += A[i][k] * B[k][j];
226  }
227  }
228  }
229 
230  return;
231 }
void G_math_f_aA_B(float **A, float **B, float a, float **C, int rows, int cols)
Add two matrices and scale matrix A with the scalar a.
Definition: blas_level_3.c:113
int rows
Definition: bitmap.h:19
#define NULL
Definition: ccmath.h:32
void G_math_d_AB(double **A, double **B, double **C, int rows_A, int cols_A, int cols_B)
Matrix multiplication.
Definition: blas_level_3.c:174
void G_math_f_AB(float **A, float **B, float **C, int rows_A, int cols_A, int cols_B)
Matrix multiplication.
Definition: blas_level_3.c:215
int cols
Definition: bitmap.h:20
void G_math_d_aA_B(double **A, double **B, double a, double **C, int rows, int cols)
Add two matrices and scale matrix A with the scalar a.
Definition: blas_level_3.c:50