numpy performance reached on reasonable matrix size
This commit is contained in:
rdgarce
@@ -1,2 +1,3 @@
|
||||
.swp
|
||||
.vscode
|
||||
*.txt
|
||||
+128
-31
@@ -3,38 +3,43 @@
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <assert.h>
|
||||
#include "../src/lina.h"
|
||||
|
||||
#define A_ROWS 1000llu
|
||||
#define A_COLS 146llu
|
||||
#define A_ROWS 960llu
|
||||
#define A_COLS 960llu
|
||||
|
||||
#define B_ROWS 146llu
|
||||
#define B_COLS 1024llu
|
||||
#define B_ROWS 960llu
|
||||
#define B_COLS 960llu
|
||||
|
||||
uint64_t nanos();
|
||||
int saveMatrixToStream(FILE *fp, double *A, int width, int height, char **error);
|
||||
static uint64_t nanos();
|
||||
|
||||
int main()
|
||||
{
|
||||
uint64_t ops = A_ROWS*B_COLS*2*A_COLS;
|
||||
uint64_t start,stop,lina_dot_time, lina_dot_mod1_time, lina_dot_mod1_1_time, lina_dot_mod2_time;
|
||||
double *A = (double *)malloc(sizeof(double)*A_ROWS*A_COLS);
|
||||
double *B = (double *)malloc(sizeof(double)*B_ROWS*B_COLS);
|
||||
uint64_t start,stop,lina_dot_time, lina_dot1_time, lina_dot2_time, lina_dot3_time, lina_dot4_time;
|
||||
|
||||
double *C1 = (double *)malloc(sizeof(double)*A_ROWS*B_COLS);
|
||||
double *C2 = (double *)malloc(sizeof(double)*A_ROWS*B_COLS);
|
||||
double *C3 = (double *)malloc(sizeof(double)*A_ROWS*B_COLS);
|
||||
double *C4 = (double *)malloc(sizeof(double)*A_ROWS*B_COLS);
|
||||
double *A = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*A_COLS);
|
||||
double *B = (double *)aligned_alloc(32,sizeof(double)*B_ROWS*B_COLS);
|
||||
|
||||
double *C1 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
|
||||
double *C2 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
|
||||
double *C3 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
|
||||
double *C4 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
|
||||
double *C5 = (double *)aligned_alloc(32,sizeof(double)*A_ROWS*B_COLS);
|
||||
|
||||
for (int i = 0; i < A_ROWS*A_COLS; i++)
|
||||
A[i] = (double)(rand()%10);
|
||||
A[i] = (double)(rand()%2);
|
||||
for (int i = 0; i < B_ROWS*B_COLS; i++)
|
||||
B[i] = (double)(rand()%10);
|
||||
B[i] = (double)(rand()%2);
|
||||
for (int i = 0; i < A_ROWS*B_COLS; i++)
|
||||
{
|
||||
C1[i] = (double)(rand()%10);
|
||||
C2[i] = (double)(rand()%10);
|
||||
C3[i] = (double)(rand()%10);
|
||||
C4[i] = (double)(rand()%10);
|
||||
C1[i] = (double)(rand()%2);
|
||||
C2[i] = (double)(rand()%2);
|
||||
C3[i] = (double)(rand()%2);
|
||||
C4[i] = (double)(rand()%2);
|
||||
C5[i] = (double)(rand()%2);
|
||||
}
|
||||
|
||||
start = nanos();
|
||||
@@ -44,36 +49,78 @@ int main()
|
||||
lina_dot_time = stop-start;
|
||||
|
||||
start = nanos();
|
||||
lina_dot_mod1(A,B,C2,A_ROWS,A_COLS,B_COLS);
|
||||
lina_dot1(A,B,C2,A_ROWS,A_COLS,B_COLS);
|
||||
stop = nanos();
|
||||
|
||||
lina_dot_mod1_time = stop-start;
|
||||
lina_dot1_time = stop-start;
|
||||
|
||||
start = nanos();
|
||||
lina_dot_mod2(A,B,C3,A_ROWS,A_COLS,B_COLS);
|
||||
lina_dot2(A,B,C3,A_ROWS,A_COLS,B_COLS);
|
||||
stop = nanos();
|
||||
|
||||
lina_dot_mod2_time = stop-start;
|
||||
lina_dot2_time = stop-start;
|
||||
|
||||
start = nanos();
|
||||
lina_dot_mod2_old(A,B,C4,A_ROWS,A_COLS,B_COLS);
|
||||
lina_dot3(A,B,C4,A_ROWS,A_COLS,B_COLS);
|
||||
stop = nanos();
|
||||
|
||||
lina_dot_mod1_1_time = stop-start;
|
||||
lina_dot3_time = stop-start;
|
||||
|
||||
if(!memcmp(C1,C2,sizeof(double)*A_ROWS*B_COLS) && !memcmp(C2,C3,sizeof(double)*A_ROWS*B_COLS) && !memcmp(C3,C4,sizeof(double)*A_ROWS*B_COLS))
|
||||
start = nanos();
|
||||
lina_dot4(A,B,C5,A_ROWS,A_COLS,B_COLS);
|
||||
stop = nanos();
|
||||
|
||||
lina_dot4_time = stop-start;
|
||||
|
||||
if(!memcmp(C1,C2,sizeof(double)*A_ROWS*B_COLS)
|
||||
&& !memcmp(C2,C3,sizeof(double)*A_ROWS*B_COLS)
|
||||
&& !memcmp(C3,C4,sizeof(double)*A_ROWS*B_COLS)
|
||||
&& !memcmp(C4,C5,sizeof(double)*A_ROWS*B_COLS))
|
||||
{
|
||||
printf( "lina_dot : %f GFLOPS\n"
|
||||
"lina_dot_mod1: %f GFLOPS\n"
|
||||
"lina_dot_mod2: %f GFLOPS\n"
|
||||
"lina_dot_mod2_old: %f GFLOPS\n", (double)ops/lina_dot_time,
|
||||
(double)ops/lina_dot_mod1_time,
|
||||
(double)ops/lina_dot_mod2_time,
|
||||
(double)ops/lina_dot_mod1_1_time);
|
||||
"lina_dot1: %f GFLOPS\n"
|
||||
"lina_dot2: %f GFLOPS\n"
|
||||
"lina_dot3: %f GFLOPS\n"
|
||||
"lina_dot4: %f GFLOPS\n", (double)ops/lina_dot_time,
|
||||
(double)ops/lina_dot1_time,
|
||||
(double)ops/lina_dot2_time,
|
||||
(double)ops/lina_dot3_time,
|
||||
(double)ops/lina_dot4_time);
|
||||
FILE *fp = fopen("lina_dots_success.txt", "w");
|
||||
if (!fp)
|
||||
return -1;
|
||||
|
||||
saveMatrixToStream(fp,C1,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C2,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C3,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C4,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C5,A_ROWS,A_COLS,NULL);
|
||||
|
||||
fclose(fp);
|
||||
}
|
||||
else
|
||||
{
|
||||
printf("ERRORE: i prodotti matriciali sono diversi!\n");
|
||||
FILE *fp = fopen("lina_dots_error.txt", "w");
|
||||
if (!fp)
|
||||
return -1;
|
||||
|
||||
saveMatrixToStream(fp,C1,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C2,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C3,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C4,A_ROWS,A_COLS,NULL);
|
||||
fprintf(fp,"\nFINE\n");
|
||||
saveMatrixToStream(fp,C5,A_ROWS,A_COLS,NULL);
|
||||
|
||||
fclose(fp);
|
||||
}
|
||||
|
||||
free(A);
|
||||
free(B);
|
||||
@@ -81,4 +128,54 @@ int main()
|
||||
free(C2);
|
||||
free(C3);
|
||||
free(C4);
|
||||
free(C5);
|
||||
}
|
||||
|
||||
static uint64_t nanos()
|
||||
{
|
||||
struct timespec time;
|
||||
clock_gettime(CLOCK_MONOTONIC, &time);
|
||||
return (uint64_t)time.tv_sec*1000000000 + (uint64_t)time.tv_nsec;
|
||||
}
|
||||
|
||||
int saveMatrixToStream(FILE *fp, double *A, int width, int height, char **error)
|
||||
{
|
||||
assert(A != NULL);
|
||||
|
||||
char *dummy;
|
||||
if (error == NULL)
|
||||
error = &dummy;
|
||||
else
|
||||
*error = NULL;
|
||||
|
||||
if (width < 1) {
|
||||
*error = "The provided width is less than one";
|
||||
return -1;
|
||||
}
|
||||
|
||||
if (height < 1) {
|
||||
*error = "The provided height is less than one";
|
||||
return -1;
|
||||
}
|
||||
|
||||
if (fp == NULL)
|
||||
fp = stdout;
|
||||
|
||||
putc('[',fp);
|
||||
|
||||
for (int i = 0; i < height-1; i++) {
|
||||
for (int j = 0; j < width-1; j++)
|
||||
fprintf(fp, "%.1f ", A[i*width + j]);
|
||||
|
||||
fprintf(fp, "%.1f,\n", A[i*width + width-1]);
|
||||
}
|
||||
|
||||
for (int j = 0; j < width-1; j++)
|
||||
fprintf(fp, "%.1f ", A[(height-1)*width + j]);
|
||||
|
||||
fprintf(fp, "%.1f", A[(height-1)*width + width-1]);
|
||||
|
||||
putc(']',fp);
|
||||
|
||||
return 0;
|
||||
}
|
||||
@@ -0,0 +1,3 @@
|
||||
gcc bench_dot.c ../src/lina.c -O3 -march=native -ffast-math -funroll-loops -o bench_dot
|
||||
./bench_dot
|
||||
python3 py_dot.py
|
||||
@@ -0,0 +1,22 @@
|
||||
import os
|
||||
os.environ['OMP_NUM_THREADS'] = '1'
|
||||
|
||||
import numpy as np
|
||||
import time
|
||||
|
||||
N = 1024
|
||||
|
||||
if __name__ == "__main__":
|
||||
A = np.random.randn(N,N).astype(np.float64)
|
||||
B = np.random.randn(N,N).astype(np.float64)
|
||||
|
||||
start = time.monotonic()
|
||||
C = A @ B
|
||||
stop = time.monotonic()
|
||||
|
||||
s = stop-start
|
||||
|
||||
ops = 2*N*N*N
|
||||
|
||||
print(f"NUMPY: {ops/s * 1e-9} GFLOPS\n")
|
||||
|
||||
+202
-1
@@ -4,6 +4,11 @@
|
||||
#include <errno.h>
|
||||
#include <ctype.h>
|
||||
#include "lina.h"
|
||||
#include <immintrin.h>
|
||||
#include <stdint.h>
|
||||
|
||||
static void
|
||||
dot_kernel_6x8(double *A_sub, double *B_sub, double *C_sub, int x, int y, int c_min, int c_max, int n, int l);
|
||||
|
||||
/* Function: lina_dot
|
||||
**
|
||||
@@ -174,7 +179,7 @@ void lina_dot2(double *A, double *B, double *C, int m, int n, int l)
|
||||
|
||||
// 2. Copy block to C
|
||||
for (int i = 0; i < BLOCKSIZE; i++)
|
||||
memcpy(&block[i*BLOCKSIZE],&C[(i+br)*l + bc], sizeof(double)*BLOCKSIZE);
|
||||
memcpy(&C[(i+br)*l + bc],&block[i*BLOCKSIZE], sizeof(double)*BLOCKSIZE);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -218,6 +223,202 @@ void lina_dot2(double *A, double *B, double *C, int m, int n, int l)
|
||||
}
|
||||
}
|
||||
|
||||
/* Function: lina_dot3
|
||||
**
|
||||
** Evaluates the dot product C = A * B. The A,B
|
||||
** matrices are, respectively, mxn and nxl, which
|
||||
** means C is mxl. The resulting C matrix is stored
|
||||
** in a memory region specified by the caller.
|
||||
**
|
||||
** Variant 3 of lina_dot:
|
||||
** This include the changes of lina_dot2 but uses
|
||||
** simd instructions to compute products and sums.
|
||||
**
|
||||
** Notes:
|
||||
**
|
||||
** - A,B must be provided as contiguous memory regions
|
||||
** represented in row-major order. Also, C is stored
|
||||
** that way too.
|
||||
**
|
||||
** - The C pointer CAN'T refer to the same memory region
|
||||
** of either A or B.
|
||||
**
|
||||
** - m,n,l must be greater than 0.
|
||||
**
|
||||
** - This function can never fail.
|
||||
*/
|
||||
void lina_dot3(double *A, double *B, double *C, int m, int n, int l)
|
||||
{
|
||||
assert(m > 0 && n > 0 && l > 0);
|
||||
assert(A != NULL && B != NULL && C != NULL);
|
||||
assert(A != C && B != C);
|
||||
|
||||
// This size is based on experimental results
|
||||
#define BLOCK_ROWS 6
|
||||
#define BLOCK_COLS 8
|
||||
|
||||
const int br_max = (m & ~(BLOCK_ROWS - 1));
|
||||
const int bc_max = (l & ~(BLOCK_COLS - 1));
|
||||
|
||||
__m256d *Bm = (__m256d *)B;
|
||||
__m256d *Cm = (__m256d *)C;
|
||||
|
||||
// problema: B non è allineato a 32 byte, cosa che pare essere il problema
|
||||
|
||||
// Dealing with the squared submatrix of C
|
||||
for (int br = 0; br < br_max; br += BLOCK_ROWS)
|
||||
{
|
||||
for (int bc = 0; bc < bc_max; bc += BLOCK_COLS)
|
||||
{
|
||||
__m256d mblock[BLOCK_ROWS][BLOCK_COLS/4] = {0};
|
||||
|
||||
// 1. Compute block
|
||||
|
||||
for(int j=0; j < n; j++)
|
||||
{
|
||||
for(int i = 0; i < BLOCK_ROWS; i++)
|
||||
{
|
||||
__m256d A_brdcst = _mm256_broadcast_sd(&A[(i+br) * n + j]);
|
||||
for(int k = 0; k < BLOCK_COLS/4; k++)
|
||||
{
|
||||
mblock[i][k] = _mm256_fmadd_pd(A_brdcst, Bm[(j * l + bc)/4 + k], mblock[i][k]);
|
||||
}
|
||||
}
|
||||
// Iteration over A's rows
|
||||
}
|
||||
|
||||
// 2. Copy block to C
|
||||
for (int i = 0; i < BLOCK_ROWS; i++)
|
||||
for (int j = 0; j < BLOCK_COLS/4; j++)
|
||||
Cm[((i+br)*l + bc)/4 + j] = mblock[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
// Dealing with the last rows and cols
|
||||
//printf("br_max: %d\nbc_max: %d\n",br_max,bc_max);
|
||||
// Last rows
|
||||
// Iteration over A's rows
|
||||
for(int i = br_max; i < m; i++) {
|
||||
// Iteration over B's columns
|
||||
for(int k = 0; k < l; k++)
|
||||
C[i*l + k] = A[i * n ] * B[k];
|
||||
}
|
||||
|
||||
// Last cols
|
||||
// Iteration over A's rows
|
||||
for (int i = 0; i < br_max; i++)
|
||||
{
|
||||
// Iteration over B's columns
|
||||
for(int k = bc_max; k < l; k++)
|
||||
C[i*l + k] = A[i * n] * B[k];
|
||||
}
|
||||
|
||||
// Iteration over the single B column
|
||||
// for executing the product of sum
|
||||
for(int j=1; j < n; j++)
|
||||
{
|
||||
// Iteration over A's rows
|
||||
for(int i = br_max; i < m; i++) {
|
||||
// Iteration over B's columns
|
||||
for(int k = 0; k < l; k++)
|
||||
C[i*l + k] += A[i * n + j] * B[j * l + k];
|
||||
}
|
||||
|
||||
// Iteration over A's rows
|
||||
for (int i = 0; i < br_max; i++)
|
||||
{
|
||||
// Iteration over B's columns
|
||||
for(int k = bc_max; k < l; k++)
|
||||
C[i*l + k] += A[i * n + j] * B[j * l + k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Function: lina_dot4
|
||||
**
|
||||
** Evaluates the dot product C = A * B. The A,B
|
||||
** matrices are, respectively, mxn and nxl, which
|
||||
** means C is mxl. The resulting C matrix is stored
|
||||
** in a memory region specified by the caller.
|
||||
**
|
||||
** Variant 4 of lina_dot:
|
||||
** This include the changes of lina_dot3 but uses the
|
||||
** micro kernel subroutine
|
||||
**
|
||||
** Notes:
|
||||
**
|
||||
** - A,B must be provided as contiguous memory regions
|
||||
** represented in row-major order. Also, C is stored
|
||||
** that way too.
|
||||
**
|
||||
** - The C pointer CAN'T refer to the same memory region
|
||||
** of either A or B.
|
||||
**
|
||||
** - m,n,l must be greater than 0.
|
||||
**
|
||||
** - This function can never fail.
|
||||
*/
|
||||
void lina_dot4(double *A, double *B, double *C, int m, int n, int l)
|
||||
{
|
||||
assert(m > 0 && n > 0 && l > 0);
|
||||
assert(A != NULL && B != NULL && C != NULL);
|
||||
assert(A != C && B != C);
|
||||
// A_sub, B_sub and C_sub must be 32 byte aligned
|
||||
assert(!((uintptr_t)A & 31llu) && !((uintptr_t)B & 31llu) && !((uintptr_t)C & 31llu));
|
||||
|
||||
#define KERNEL_ROW 6
|
||||
#define KERNEL_COLS 8
|
||||
|
||||
const int br_max = (m & ~(KERNEL_ROW - 1));
|
||||
const int bc_max = (l & ~(KERNEL_COLS - 1));
|
||||
|
||||
for (int br = 0; br < br_max; br += KERNEL_ROW)
|
||||
{
|
||||
for (int bc = 0; bc < bc_max; bc += KERNEL_COLS)
|
||||
{
|
||||
dot_kernel_6x8(A, B, C, br, bc, 0, n, n, l);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
*
|
||||
* Computes C_sub += A_sub * B_sub where:
|
||||
* - C_sub = C[x:x+6][y:y+8]
|
||||
* - A_sub = A[x:x+6][c_min:c_max]
|
||||
* - B_sub = B[c_min:c_max][y:y+8]
|
||||
* - n is the number of columns of A
|
||||
* - l the number of columns of B
|
||||
*/
|
||||
static void
|
||||
dot_kernel_6x8(double *A_sub, double *B_sub, double *C_sub, int x, int y, int c_min, int c_max, int n, int l)
|
||||
{
|
||||
// A_sub, B_sub and C_sub must be 32 byte aligned
|
||||
// assert is done in the main lina_dot function
|
||||
//assert(!((uintptr_t)A_sub & 31llu) && !((uintptr_t)B_sub & 31llu) && !((uintptr_t)C_sub & 31llu));
|
||||
|
||||
// This structure should use 12 YMM registers
|
||||
|
||||
__m256d *Bm_sub = (__m256d *)B_sub;
|
||||
__m256d *Cm_sub = (__m256d *)C_sub;
|
||||
__m256d acc[6][2] = {0};
|
||||
|
||||
|
||||
for (int k = c_min; k < c_max; k++)
|
||||
{
|
||||
for (int i = 0; i < 6; i++)
|
||||
{
|
||||
__m256d A_brdcst = _mm256_broadcast_sd(&A_sub[(x + i)*n + k]);
|
||||
for (int j = 0; j < 2; j++)
|
||||
acc[i][j] = _mm256_fmadd_pd(A_brdcst,Bm_sub[(k*l + y)/4 + j],acc[i][j]);
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < 6; i++)
|
||||
for (int j = 0; j < 2; j++)
|
||||
Cm_sub[((x+i)*l + y)/4 + j] = acc[i][j];
|
||||
}
|
||||
|
||||
/* Function: lina_add
|
||||
**
|
||||
** Evaluates the matrix addition C = A + B. The result
|
||||
|
||||
@@ -4,6 +4,8 @@
|
||||
void lina_dot(double *A, double *B, double *C, int m, int n, int l);
|
||||
void lina_dot1(double *A, double *B, double *C, int m, int n, int l);
|
||||
void lina_dot2(double *A, double *B, double *C, int m, int n, int l);
|
||||
void lina_dot3(double *A, double *B, double *C, int m, int n, int l);
|
||||
void lina_dot4(double *A, double *B, double *C, int m, int n, int l);
|
||||
void lina_add(double *A, double *B, double *C, int m, int n);
|
||||
void lina_scale(double *A, double *B, double k, int m, int n);
|
||||
void lina_conv(double *A, double *B, double *C, int Aw, int Ah, int Bw, int Bh);
|
||||
|
||||
Reference in New Issue
Block a user