Line data Source code
1 : /*************************************************************************
2 : *
3 : * Project
4 : * _____ _____ __ __ _____
5 : * / ____| __ \| \/ | __ \
6 : * ___ _ __ ___ _ __ | | __| |__) | \ / | |__) |
7 : * / _ \| '_ \ / _ \ '_ \| | |_ | ___/| |\/| | ___/
8 : *| (_) | |_) | __/ | | | |__| | | | | | | |
9 : * \___/| .__/ \___|_| |_|\_____|_| |_| |_|_|
10 : * | |
11 : * |_|
12 : *
13 : * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
14 : *
15 : * This software is licensed as described in the file LICENSE, which
16 : * you should have received as part of this distribution. The terms
17 : * among other details are referenced in the official documentation
18 : * seen here : https://akielaries.github.io/openGPMP/ along with
19 : * important files seen in this project.
20 : *
21 : * You may opt to use, copy, modify, merge, publish, distribute
22 : * and/or sell copies of the Software, and permit persons to whom
23 : * the Software is furnished to do so, under the terms of the
24 : * LICENSE file. As this is an Open Source effort, all implementations
25 : * must be of the same methodology.
26 : *
27 : *
28 : *
29 : * This software is distributed on an AS IS basis, WITHOUT
30 : * WARRANTY OF ANY KIND, either express or implied.
31 : *
32 : ************************************************************************/
33 :
34 : /** Single precision GEneral Matrix-Matrix product */
35 : #include <cmath>
36 : #include <limits>
37 : #include <openGPMP/linalg/_sgemm.hpp>
38 : #include <stddef.h>
39 : #include <stdio.h>
40 : #include <stdlib.h>
41 : #include <time.h>
42 :
43 : // MATRIX BUFFERS
44 : float gpmp::linalg::SGEMM::SGEMM_BUFF_A[BLOCK_SZ_M * BLOCK_SZ_K];
45 : float gpmp::linalg::SGEMM::SGEMM_BUFF_B[BLOCK_SZ_K * BLOCK_SZ_N];
46 : float gpmp::linalg::SGEMM::SGEMM_BUFF_C[BLOCK_SZ_MR * BLOCK_SZ_NR];
47 :
48 : // pack micro panels of size BLOCK_SZ_MR rows by k columns from A without
49 : // padding
50 768 : void gpmp::linalg::SGEMM::pack_micro_A(int k,
51 : const float *A,
52 : int incRowA,
53 : int incColA,
54 : float *buffer) {
55 : int i, j;
56 :
57 262912 : for (j = 0; j < k; ++j) {
58 1310720 : for (i = 0; i < BLOCK_SZ_MR; ++i) {
59 1048576 : buffer[i] = A[i * incRowA];
60 : }
61 262144 : buffer += BLOCK_SZ_MR;
62 262144 : A += incColA;
63 : }
64 768 : }
65 :
66 : // packs panels from A with padding if needed
67 9 : void gpmp::linalg::SGEMM::pack_buffer_A(int mc,
68 : int kc,
69 : const float *A,
70 : int incRowA,
71 : int incColA,
72 : float *buffer) {
73 9 : int mp = mc / BLOCK_SZ_MR;
74 9 : int _mr = mc % BLOCK_SZ_MR;
75 :
76 : int i, j;
77 :
78 777 : for (i = 0; i < mp; ++i) {
79 768 : pack_micro_A(kc, A, incRowA, incColA, buffer);
80 768 : buffer += kc * BLOCK_SZ_MR;
81 768 : A += BLOCK_SZ_MR * incRowA;
82 : }
83 9 : if (_mr > 0) {
84 0 : for (j = 0; j < kc; ++j) {
85 0 : for (i = 0; i < _mr; ++i) {
86 0 : buffer[i] = A[i * incRowA];
87 : }
88 0 : for (i = _mr; i < BLOCK_SZ_MR; ++i) {
89 0 : buffer[i] = 0.0f;
90 : }
91 0 : buffer += BLOCK_SZ_MR;
92 0 : A += incColA;
93 : }
94 : }
95 9 : }
96 :
97 : // packing complete panels from B of size BLOCK_SZ_NR by k columns
98 768 : void gpmp::linalg::SGEMM::pack_micro_B(int k,
99 : const float *B,
100 : int incRowB,
101 : int incColB,
102 : float *buffer) {
103 : int i, j;
104 :
105 262912 : for (i = 0; i < k; ++i) {
106 1310720 : for (j = 0; j < BLOCK_SZ_NR; ++j) {
107 1048576 : buffer[j] = B[j * incColB];
108 : }
109 262144 : buffer += BLOCK_SZ_NR;
110 262144 : B += incRowB;
111 : }
112 768 : }
113 :
114 : // packing panels from B with padding if needed
115 3 : void gpmp::linalg::SGEMM::pack_buffer_B(int kc,
116 : int nc,
117 : const float *B,
118 : int incRowB,
119 : int incColB,
120 : float *buffer) {
121 3 : int np = nc / BLOCK_SZ_NR;
122 3 : int _nr = nc % BLOCK_SZ_NR;
123 :
124 : int i, j;
125 :
126 771 : for (j = 0; j < np; ++j) {
127 768 : pack_micro_B(kc, B, incRowB, incColB, buffer);
128 768 : buffer += kc * BLOCK_SZ_NR;
129 768 : B += BLOCK_SZ_NR * incColB;
130 : }
131 3 : if (_nr > 0) {
132 0 : for (i = 0; i < kc; ++i) {
133 0 : for (j = 0; j < _nr; ++j) {
134 0 : buffer[j] = B[j * incColB];
135 : }
136 0 : for (j = _nr; j < BLOCK_SZ_NR; ++j) {
137 0 : buffer[j] = 0.0f;
138 : }
139 0 : buffer += BLOCK_SZ_NR;
140 0 : B += incRowB;
141 : }
142 : }
143 3 : }
144 :
145 : // micro kernel that multiplies panels from A and B
146 196608 : void gpmp::linalg::SGEMM::sgemm_micro_kernel(int kc,
147 : float alpha,
148 : const float *A,
149 : const float *B,
150 : float beta,
151 : float *C,
152 : int incRowC,
153 : int incColC) {
154 : float AB[BLOCK_SZ_MR * BLOCK_SZ_NR];
155 :
156 : int i, j, l;
157 :
158 : // Compute AB = A*B
159 3342336 : for (l = 0; l < BLOCK_SZ_MR * BLOCK_SZ_NR; ++l) {
160 3145728 : AB[l] = 0;
161 : }
162 67305472 : for (l = 0; l < kc; ++l) {
163 335544320 : for (j = 0; j < BLOCK_SZ_NR; ++j) {
164 1342177280 : for (i = 0; i < BLOCK_SZ_MR; ++i) {
165 1073741824 : AB[i + j * BLOCK_SZ_MR] += A[i] * B[j];
166 : }
167 : }
168 67108864 : A += BLOCK_SZ_MR;
169 67108864 : B += BLOCK_SZ_NR;
170 : }
171 :
172 : // Update C <- beta*C
173 196608 : if (fabs(beta - 0.0f) < std::numeric_limits<float>::epsilon()) {
174 327680 : for (j = 0; j < BLOCK_SZ_NR; ++j) {
175 1310720 : for (i = 0; i < BLOCK_SZ_MR; ++i) {
176 1048576 : C[i * incRowC + j * incColC] = 0.0f;
177 : }
178 : }
179 131072 : } else if (fabs(beta - 1.0f) > std::numeric_limits<float>::epsilon()) {
180 0 : for (j = 0; j < BLOCK_SZ_NR; ++j) {
181 0 : for (i = 0; i < BLOCK_SZ_MR; ++i) {
182 0 : C[i * incRowC + j * incColC] *= beta;
183 : }
184 : }
185 : }
186 :
187 : // Update C <- C + alpha*AB (note: the case alpha==0.0f was already treated
188 : // in
189 : // the above layer sgemm_nn)
190 196608 : if (fabs(alpha - 1.0f) < std::numeric_limits<float>::epsilon()) {
191 983040 : for (j = 0; j < BLOCK_SZ_NR; ++j) {
192 3932160 : for (i = 0; i < BLOCK_SZ_MR; ++i) {
193 3145728 : C[i * incRowC + j * incColC] += AB[i + j * BLOCK_SZ_MR];
194 : }
195 : }
196 : }
197 :
198 : else {
199 0 : for (j = 0; j < BLOCK_SZ_NR; ++j) {
200 0 : for (i = 0; i < BLOCK_SZ_MR; ++i) {
201 0 : C[i * incRowC + j * incColC] += alpha * AB[i + j * BLOCK_SZ_MR];
202 : }
203 : }
204 : }
205 196608 : }
206 :
207 : // Compute Y += alpha*X (float precision AX + Y)
208 0 : void gpmp::linalg::SGEMM::sgeaxpy(int m,
209 : int n,
210 : float alpha,
211 : const float *X,
212 : int incRowX,
213 : int incColX,
214 : float *Y,
215 : int incRowY,
216 : int incColY) {
217 : int i, j;
218 :
219 0 : if (fabs(alpha - 1.0f) > std::numeric_limits<float>::epsilon()) {
220 :
221 0 : for (j = 0; j < n; ++j) {
222 0 : for (i = 0; i < m; ++i) {
223 0 : Y[i * incRowY + j * incColY] +=
224 0 : alpha * X[i * incRowX + j * incColX];
225 : }
226 : }
227 : }
228 :
229 : else {
230 0 : for (j = 0; j < n; ++j) {
231 0 : for (i = 0; i < m; ++i) {
232 0 : Y[i * incRowY + j * incColY] += X[i * incRowX + j * incColX];
233 : }
234 : }
235 : }
236 0 : }
237 :
238 : // Compute X *= alpha (scale elements)
239 0 : void gpmp::linalg::SGEMM::sgescal(int m,
240 : int n,
241 : float alpha,
242 : float *X,
243 : int incRowX,
244 : int incColX) {
245 : int i, j;
246 :
247 0 : if (fabs(alpha - 0.0f) > std::numeric_limits<float>::epsilon()) {
248 0 : for (j = 0; j < n; ++j) {
249 0 : for (i = 0; i < m; ++i) {
250 0 : X[i * incRowX + j * incColX] *= alpha;
251 : }
252 : }
253 : }
254 :
255 : else {
256 0 : for (j = 0; j < n; ++j) {
257 0 : for (i = 0; i < m; ++i) {
258 0 : X[i * incRowX + j * incColX] = 0.0f;
259 : }
260 : }
261 : }
262 0 : }
263 :
264 : // Macro Kernel for the multiplication of blocks of A and B. We assume that
265 : // these blocks were previously packed to buffers SGEMM_BUFF_A and SGEMM_BUFF_B.
266 9 : void gpmp::linalg::SGEMM::sgemm_macro_kernel(int mc,
267 : int nc,
268 : int kc,
269 : float alpha,
270 : float beta,
271 : float *C,
272 : int incRowC,
273 : int incColC) {
274 :
275 9 : int mp = (mc + BLOCK_SZ_MR - 1) / BLOCK_SZ_MR;
276 9 : int np = (nc + BLOCK_SZ_NR - 1) / BLOCK_SZ_NR;
277 :
278 9 : int _mr = mc % BLOCK_SZ_MR;
279 9 : int _nr = nc % BLOCK_SZ_NR;
280 :
281 : int mr, nr;
282 : int i, j;
283 :
284 2313 : for (j = 0; j < np; ++j) {
285 2304 : nr = (j != np - 1 || _nr == 0) ? BLOCK_SZ_NR : _nr;
286 :
287 198912 : for (i = 0; i < mp; ++i) {
288 196608 : mr = (i != mp - 1 || _mr == 0) ? BLOCK_SZ_MR : _mr;
289 :
290 196608 : if (mr == BLOCK_SZ_MR && nr == BLOCK_SZ_NR) {
291 196608 : sgemm_micro_kernel(
292 : kc,
293 : alpha,
294 196608 : &SGEMM_BUFF_A[i * kc * BLOCK_SZ_MR],
295 196608 : &SGEMM_BUFF_B[j * kc * BLOCK_SZ_NR],
296 : beta,
297 196608 : &C[i * BLOCK_SZ_MR * incRowC + j * BLOCK_SZ_NR * incColC],
298 : incRowC,
299 : incColC);
300 : } else {
301 0 : sgemm_micro_kernel(kc,
302 : alpha,
303 0 : &SGEMM_BUFF_A[i * kc * BLOCK_SZ_MR],
304 0 : &SGEMM_BUFF_B[j * kc * BLOCK_SZ_NR],
305 : 0.0f,
306 : SGEMM_BUFF_C,
307 : 1,
308 : BLOCK_SZ_MR);
309 0 : sgescal(
310 : mr,
311 : nr,
312 : beta,
313 0 : &C[i * BLOCK_SZ_MR * incRowC + j * BLOCK_SZ_NR * incColC],
314 : incRowC,
315 : incColC);
316 0 : sgeaxpy(
317 : mr,
318 : nr,
319 : 1.0f,
320 : SGEMM_BUFF_C,
321 : 1,
322 : BLOCK_SZ_MR,
323 0 : &C[i * BLOCK_SZ_MR * incRowC + j * BLOCK_SZ_NR * incColC],
324 : incRowC,
325 : incColC);
326 : }
327 : }
328 : }
329 9 : }
330 :
331 : // Main SGEMM entrypoint, compute C <- beta*C + alpha*A*B
332 1 : void gpmp::linalg::SGEMM::sgemm_nn(int m,
333 : int n,
334 : int k,
335 : float alpha,
336 : const float *A,
337 : int incRowA,
338 : int incColA,
339 : const float *B,
340 : int incRowB,
341 : int incColB,
342 : float beta,
343 : float *C,
344 : int incRowC,
345 : int incColC) {
346 1 : int mb = (m + BLOCK_SZ_M - 1) / BLOCK_SZ_M;
347 1 : int nb = (n + BLOCK_SZ_N - 1) / BLOCK_SZ_N;
348 1 : int kb = (k + BLOCK_SZ_K - 1) / BLOCK_SZ_K;
349 :
350 1 : int _mc = m % BLOCK_SZ_M;
351 1 : int _nc = n % BLOCK_SZ_N;
352 1 : int _kc = k % BLOCK_SZ_K;
353 :
354 : int mc, nc, kc;
355 : int i, j, l;
356 :
357 : float _beta;
358 :
359 1 : if (fabs(alpha) < std::numeric_limits<float>::epsilon() || k == 0) {
360 0 : sgescal(m, n, beta, C, incRowC, incColC);
361 0 : return;
362 : }
363 :
364 2 : for (j = 0; j < nb; ++j) {
365 1 : nc = (j != nb - 1 || _nc == 0) ? BLOCK_SZ_N : _nc;
366 :
367 4 : for (l = 0; l < kb; ++l) {
368 3 : kc = (l != kb - 1 || _kc == 0) ? BLOCK_SZ_K : _kc;
369 3 : _beta = (l == 0) ? beta : 1.0f;
370 :
371 3 : pack_buffer_B(
372 : kc,
373 : nc,
374 3 : &B[l * BLOCK_SZ_K * incRowB + j * BLOCK_SZ_N * incColB],
375 : incRowB,
376 : incColB,
377 : SGEMM_BUFF_B);
378 :
379 12 : for (i = 0; i < mb; ++i) {
380 9 : mc = (i != mb - 1 || _mc == 0) ? BLOCK_SZ_M : _mc;
381 :
382 9 : pack_buffer_A(
383 : mc,
384 : kc,
385 9 : &A[i * BLOCK_SZ_M * incRowA + l * BLOCK_SZ_K * incColA],
386 : incRowA,
387 : incColA,
388 : SGEMM_BUFF_A);
389 :
390 9 : sgemm_macro_kernel(
391 : mc,
392 : nc,
393 : kc,
394 : alpha,
395 : _beta,
396 9 : &C[i * BLOCK_SZ_M * incRowC + j * BLOCK_SZ_N * incColC],
397 : incRowC,
398 : incColC);
399 : }
400 : }
401 : }
402 : }
|
