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