GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
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 * */
50void 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 * */
113void 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 * */
174void 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 * */
215void 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
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
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
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
#define NULL
Definition: ccmath.h:32