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 : /** Integer type GEneral Matrix-Matrix product */
35 : #include <cmath>
36 : #include <limits>
37 : #include <openGPMP/linalg/_igemm.hpp>
38 : #include <stddef.h>
39 : #include <stdio.h>
40 : #include <stdlib.h>
41 : #include <time.h>
42 :
43 : // MATRIX BUFFERS
44 : int gpmp::linalg::IGEMM::IGEMM_BUFF_A[BLOCK_SZ_M * BLOCK_SZ_K];
45 : int gpmp::linalg::IGEMM::IGEMM_BUFF_B[BLOCK_SZ_K * BLOCK_SZ_N];
46 : int gpmp::linalg::IGEMM::IGEMM_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::IGEMM::pack_micro_A(int k,
51 : const int *A,
52 : int incRowA,
53 : int incColA,
54 : int *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::IGEMM::pack_buffer_A(int mc,
68 : int kc,
69 : const int *A,
70 : int incRowA,
71 : int incColA,
72 : int *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;
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::IGEMM::pack_micro_B(int k,
99 : const int *B,
100 : int incRowB,
101 : int incColB,
102 : int *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::IGEMM::pack_buffer_B(int kc,
116 : int nc,
117 : const int *B,
118 : int incRowB,
119 : int incColB,
120 : int *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;
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::IGEMM::igemm_micro_kernel(int kc,
147 : int alpha,
148 : const int *A,
149 : const int *B,
150 : int beta,
151 : int *C,
152 : int incRowC,
153 : int incColC) {
154 : int 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 (beta == 0) {
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;
177 : }
178 : }
179 131072 : } else if (beta != 1) {
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 was already treated
188 : // in
189 : // the above layer igemm_nn)
190 196608 : if (alpha == 1) {
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 (int precision AX + Y)
208 0 : void gpmp::linalg::IGEMM::igeaxpy(int m,
209 : int n,
210 : int alpha,
211 : const int *X,
212 : int incRowX,
213 : int incColX,
214 : int *Y,
215 : int incRowY,
216 : int incColY) {
217 : int i, j;
218 :
219 0 : if (alpha != 1) {
220 0 : for (j = 0; j < n; ++j) {
221 0 : for (i = 0; i < m; ++i) {
222 0 : Y[i * incRowY + j * incColY] +=
223 0 : alpha * X[i * incRowX + j * incColX];
224 : }
225 : }
226 : }
227 :
228 : else {
229 0 : for (j = 0; j < n; ++j) {
230 0 : for (i = 0; i < m; ++i) {
231 0 : Y[i * incRowY + j * incColY] += X[i * incRowX + j * incColX];
232 : }
233 : }
234 : }
235 0 : }
236 :
237 : // Compute X *= alpha (scale elements)
238 0 : void gpmp::linalg::IGEMM::igescal(int m,
239 : int n,
240 : int alpha,
241 : int *X,
242 : int incRowX,
243 : int incColX) {
244 : int i, j;
245 :
246 0 : if (alpha != 0) {
247 0 : for (j = 0; j < n; ++j) {
248 0 : for (i = 0; i < m; ++i) {
249 0 : X[i * incRowX + j * incColX] *= alpha;
250 : }
251 : }
252 : }
253 :
254 : else {
255 0 : for (j = 0; j < n; ++j) {
256 0 : for (i = 0; i < m; ++i) {
257 0 : X[i * incRowX + j * incColX] = 0;
258 : }
259 : }
260 : }
261 0 : }
262 :
263 : // Macro Kernel for the multiplication of blocks of A and B. We assume that
264 : // these blocks were previously packed to buffers IGEMM_BUFF_A and IGEMM_BUFF_B.
265 9 : void gpmp::linalg::IGEMM::igemm_macro_kernel(int mc,
266 : int nc,
267 : int kc,
268 : int alpha,
269 : int beta,
270 : int *C,
271 : int incRowC,
272 : int incColC) {
273 :
274 9 : int mp = (mc + BLOCK_SZ_MR - 1) / BLOCK_SZ_MR;
275 9 : int np = (nc + BLOCK_SZ_NR - 1) / BLOCK_SZ_NR;
276 :
277 9 : int _mr = mc % BLOCK_SZ_MR;
278 9 : int _nr = nc % BLOCK_SZ_NR;
279 :
280 : int mr, nr;
281 : int i, j;
282 :
283 2313 : for (j = 0; j < np; ++j) {
284 2304 : nr = (j != np - 1 || _nr == 0) ? BLOCK_SZ_NR : _nr;
285 :
286 198912 : for (i = 0; i < mp; ++i) {
287 196608 : mr = (i != mp - 1 || _mr == 0) ? BLOCK_SZ_MR : _mr;
288 :
289 196608 : if (mr == BLOCK_SZ_MR && nr == BLOCK_SZ_NR) {
290 196608 : igemm_micro_kernel(
291 : kc,
292 : alpha,
293 196608 : &IGEMM_BUFF_A[i * kc * BLOCK_SZ_MR],
294 196608 : &IGEMM_BUFF_B[j * kc * BLOCK_SZ_NR],
295 : beta,
296 196608 : &C[i * BLOCK_SZ_MR * incRowC + j * BLOCK_SZ_NR * incColC],
297 : incRowC,
298 : incColC);
299 : } else {
300 0 : igemm_micro_kernel(kc,
301 : alpha,
302 0 : &IGEMM_BUFF_A[i * kc * BLOCK_SZ_MR],
303 0 : &IGEMM_BUFF_B[j * kc * BLOCK_SZ_NR],
304 : 0,
305 : IGEMM_BUFF_C,
306 : 1,
307 : BLOCK_SZ_MR);
308 0 : igescal(
309 : mr,
310 : nr,
311 : beta,
312 0 : &C[i * BLOCK_SZ_MR * incRowC + j * BLOCK_SZ_NR * incColC],
313 : incRowC,
314 : incColC);
315 0 : igeaxpy(
316 : mr,
317 : nr,
318 : 1.0,
319 : IGEMM_BUFF_C,
320 : 1,
321 : BLOCK_SZ_MR,
322 0 : &C[i * BLOCK_SZ_MR * incRowC + j * BLOCK_SZ_NR * incColC],
323 : incRowC,
324 : incColC);
325 : }
326 : }
327 : }
328 9 : }
329 :
330 : // Main IGEMM entrypoint, compute C <- beta*C + alpha*A*B
331 1 : void gpmp::linalg::IGEMM::igemm_nn(int m,
332 : int n,
333 : int k,
334 : int alpha,
335 : const int *A,
336 : int incRowA,
337 : int incColA,
338 : const int *B,
339 : int incRowB,
340 : int incColB,
341 : int beta,
342 : int *C,
343 : int incRowC,
344 : int incColC) {
345 1 : int mb = (m + BLOCK_SZ_M - 1) / BLOCK_SZ_M;
346 1 : int nb = (n + BLOCK_SZ_N - 1) / BLOCK_SZ_N;
347 1 : int kb = (k + BLOCK_SZ_K - 1) / BLOCK_SZ_K;
348 :
349 1 : int _mc = m % BLOCK_SZ_M;
350 1 : int _nc = n % BLOCK_SZ_N;
351 1 : int _kc = k % BLOCK_SZ_K;
352 :
353 : int mc, nc, kc;
354 : int i, j, l;
355 :
356 : int _beta;
357 :
358 1 : if (alpha == 0 || k == 0) {
359 0 : igescal(m, n, beta, C, incRowC, incColC);
360 0 : return;
361 : }
362 :
363 2 : for (j = 0; j < nb; ++j) {
364 1 : nc = (j != nb - 1 || _nc == 0) ? BLOCK_SZ_N : _nc;
365 :
366 4 : for (l = 0; l < kb; ++l) {
367 3 : kc = (l != kb - 1 || _kc == 0) ? BLOCK_SZ_K : _kc;
368 3 : _beta = (l == 0) ? beta : 1.0;
369 :
370 3 : pack_buffer_B(
371 : kc,
372 : nc,
373 3 : &B[l * BLOCK_SZ_K * incRowB + j * BLOCK_SZ_N * incColB],
374 : incRowB,
375 : incColB,
376 : IGEMM_BUFF_B);
377 :
378 12 : for (i = 0; i < mb; ++i) {
379 9 : mc = (i != mb - 1 || _mc == 0) ? BLOCK_SZ_M : _mc;
380 :
381 9 : pack_buffer_A(
382 : mc,
383 : kc,
384 9 : &A[i * BLOCK_SZ_M * incRowA + l * BLOCK_SZ_K * incColA],
385 : incRowA,
386 : incColA,
387 : IGEMM_BUFF_A);
388 :
389 9 : igemm_macro_kernel(
390 : mc,
391 : nc,
392 : kc,
393 : alpha,
394 : _beta,
395 9 : &C[i * BLOCK_SZ_M * incRowC + j * BLOCK_SZ_N * incColC],
396 : incRowC,
397 : incColC);
398 : }
399 : }
400 : }
401 : }
|
