phynexis-deps/OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #5972 from Orcina-Ltd/asum-alignment-determinism

Browse Source 
Title: kernel/x86_64: make AVX-512 asum/sum kernels independent of buffer alignment
This commit is contained in:
Martin Kroeker
2026-08-16 16:42:53 +02:00
committed by GitHub
parent 377753094d d793be85b4
commit b08df585f6
6 changed files with 142 additions and 180 deletions
Download Patch File Download Diff File Expand all files Collapse all files
kernel/x86_64/casum_microk_skylakex-2.c
+31 -38
View File
@@ -130,7 +130,7 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
accum_2 = _mm512_setzero_ps();
accum_3 = _mm512_setzero_ps();
// alignment has side-effect when the size of input array is not large enough
// shorter unrolled path for small inputs
if (n2 < 256) {
if (n2 >= 128) {
x00 = _mm512_loadu_ps(&x1[ 0]);
@@ -217,29 +217,22 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
sumf = _mm512_reduce_add_ps(accum_0);
}
// n2 >= 256, doing alignment
/* n2 >= 256. Deliberately no peeling to an alignment boundary:
peeling makes the grouping of the sum into accumulators - and
hence the rounding of the result - depend on the buffer address.
Unaligned loads keep the result a function of n alone, and cost
the same as aligned loads when the address happens to be
aligned. */
else {
int align_header = ((64 - ((uintptr_t)x1 & (uintptr_t)0x3f)) >> 2) & 0xf;
if (0 != align_header) {
uint16_t align_mask16 = (((uint16_t)0xffff) >> (16 - align_header));
x00 = _mm512_maskz_loadu_ps(*((__mmask16*) &align_mask16), &x1[0]);
x00 = _mm512_and_ps(x00, abs_mask);
accum_0 = _mm512_add_ps(accum_0, x00);
n2 -= align_header;
x1 += align_header;
}
x00 = _mm512_load_ps(&x1[ 0]);
x01 = _mm512_load_ps(&x1[ 16]);
x02 = _mm512_load_ps(&x1[ 32]);
x03 = _mm512_load_ps(&x1[ 48]);
x04 = _mm512_load_ps(&x1[ 64]);
x05 = _mm512_load_ps(&x1[ 80]);
x06 = _mm512_load_ps(&x1[ 96]);
x07 = _mm512_load_ps(&x1[112]);
x00 = _mm512_loadu_ps(&x1[ 0]);
x01 = _mm512_loadu_ps(&x1[ 16]);
x02 = _mm512_loadu_ps(&x1[ 32]);
x03 = _mm512_loadu_ps(&x1[ 48]);
x04 = _mm512_loadu_ps(&x1[ 64]);
x05 = _mm512_loadu_ps(&x1[ 80]);
x06 = _mm512_loadu_ps(&x1[ 96]);
x07 = _mm512_loadu_ps(&x1[112]);
n2 -= 128;
x1 += 128;
@@ -251,26 +244,26 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
x03 = _mm512_and_ps(x03, abs_mask);
accum_0 = _mm512_add_ps(accum_0, x00);
x00 = _mm512_load_ps(&x1[ 0]);
x00 = _mm512_loadu_ps(&x1[ 0]);
accum_1 = _mm512_add_ps(accum_1, x01);
x01 = _mm512_load_ps(&x1[ 16]);
x01 = _mm512_loadu_ps(&x1[ 16]);
accum_2 = _mm512_add_ps(accum_2, x02);
x02 = _mm512_load_ps(&x1[ 32]);
x02 = _mm512_loadu_ps(&x1[ 32]);
accum_3 = _mm512_add_ps(accum_3, x03);
x03 = _mm512_load_ps(&x1[ 48]);
x03 = _mm512_loadu_ps(&x1[ 48]);
x04 = _mm512_and_ps(x04, abs_mask);
x05 = _mm512_and_ps(x05, abs_mask);
x06 = _mm512_and_ps(x06, abs_mask);
x07 = _mm512_and_ps(x07, abs_mask);
accum_0 = _mm512_add_ps(accum_0, x04);
x04 = _mm512_load_ps(&x1[ 64]);
x04 = _mm512_loadu_ps(&x1[ 64]);
accum_1 = _mm512_add_ps(accum_1, x05);
x05 = _mm512_load_ps(&x1[ 80]);
x05 = _mm512_loadu_ps(&x1[ 80]);
accum_2 = _mm512_add_ps(accum_2, x06);
x06 = _mm512_load_ps(&x1[ 96]);
x06 = _mm512_loadu_ps(&x1[ 96]);
accum_3 = _mm512_add_ps(accum_3, x07);
x07 = _mm512_load_ps(&x1[112]);
x07 = _mm512_loadu_ps(&x1[112]);
n2 -= 128;
x1 += 128;
@@ -296,10 +289,10 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
accum_3 = _mm512_add_ps(accum_3, x07);
if (n2 >= 64) {
x00 = _mm512_load_ps(&x1[ 0]);
x01 = _mm512_load_ps(&x1[16]);
x02 = _mm512_load_ps(&x1[32]);
x03 = _mm512_load_ps(&x1[48]);
x00 = _mm512_loadu_ps(&x1[ 0]);
x01 = _mm512_loadu_ps(&x1[16]);
x02 = _mm512_loadu_ps(&x1[32]);
x03 = _mm512_loadu_ps(&x1[48]);
x00 = _mm512_and_ps(x00, abs_mask);
x01 = _mm512_and_ps(x01, abs_mask);
x02 = _mm512_and_ps(x02, abs_mask);
@@ -314,8 +307,8 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 32) {
x00 = _mm512_load_ps(&x1[ 0]);
x01 = _mm512_load_ps(&x1[16]);
x00 = _mm512_loadu_ps(&x1[ 0]);
x01 = _mm512_loadu_ps(&x1[16]);
x00 = _mm512_and_ps(x00, abs_mask);
x01 = _mm512_and_ps(x01, abs_mask);
accum_0 = _mm512_add_ps(accum_0, x00);
@@ -326,7 +319,7 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 16) {
x00 = _mm512_load_ps(&x1[ 0]);
x00 = _mm512_loadu_ps(&x1[ 0]);
x00 = _mm512_and_ps(x00, abs_mask);
accum_0 = _mm512_add_ps(accum_0, x00);
@@ -336,7 +329,7 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
if (n2) {
uint16_t tail_mask16 = (((uint16_t) 0xffff) >> (16 - n2));
x00 = _mm512_maskz_load_ps(*((__mmask16*) &tail_mask16), &x1[ 0]);
x00 = _mm512_maskz_loadu_ps(*((__mmask16*) &tail_mask16), &x1[ 0]);
x00 = _mm512_and_ps(x00, abs_mask);
accum_0 = _mm512_add_ps(accum_0, x00);
}
kernel/x86_64/csum_microk_skylakex-2.c
+31 -37
View File
@@ -109,7 +109,7 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
accum_2 = _mm512_setzero_ps();
accum_3 = _mm512_setzero_ps();
// alignment has side-effect when the size of input array is not large enough
// shorter unrolled path for small inputs
if (n2 < 256) {
if (n2 >= 128) {
x00 = _mm512_loadu_ps(&x1[ 0]);
@@ -178,28 +178,22 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
sumf = _mm512_reduce_add_ps(accum_0);
}
// n2 >= 256, doing alignment
/* n2 >= 256. Deliberately no peeling to an alignment boundary:
peeling makes the grouping of the sum into accumulators - and
hence the rounding of the result - depend on the buffer address.
Unaligned loads keep the result a function of n alone, and cost
the same as aligned loads when the address happens to be
aligned. */
else {
int align_header = ((64 - ((uintptr_t)x1 & (uintptr_t)0x3f)) >> 2) & 0xf;
if (0 != align_header) {
uint16_t align_mask16 = (((uint16_t)0xffff) >> (16 - align_header));
x00 = _mm512_maskz_loadu_ps(*((__mmask16*) &align_mask16), &x1[0]);
accum_0 = _mm512_add_ps(accum_0, x00);
n2 -= align_header;
x1 += align_header;
}
x00 = _mm512_load_ps(&x1[ 0]);
x01 = _mm512_load_ps(&x1[ 16]);
x02 = _mm512_load_ps(&x1[ 32]);
x03 = _mm512_load_ps(&x1[ 48]);
x04 = _mm512_load_ps(&x1[ 64]);
x05 = _mm512_load_ps(&x1[ 80]);
x06 = _mm512_load_ps(&x1[ 96]);
x07 = _mm512_load_ps(&x1[112]);
x00 = _mm512_loadu_ps(&x1[ 0]);
x01 = _mm512_loadu_ps(&x1[ 16]);
x02 = _mm512_loadu_ps(&x1[ 32]);
x03 = _mm512_loadu_ps(&x1[ 48]);
x04 = _mm512_loadu_ps(&x1[ 64]);
x05 = _mm512_loadu_ps(&x1[ 80]);
x06 = _mm512_loadu_ps(&x1[ 96]);
x07 = _mm512_loadu_ps(&x1[112]);
n2 -= 128;
x1 += 128;
@@ -207,22 +201,22 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
while (n2 >= 128) {
accum_0 = _mm512_add_ps(accum_0, x00);
x00 = _mm512_load_ps(&x1[ 0]);
x00 = _mm512_loadu_ps(&x1[ 0]);
accum_1 = _mm512_add_ps(accum_1, x01);
x01 = _mm512_load_ps(&x1[ 16]);
x01 = _mm512_loadu_ps(&x1[ 16]);
accum_2 = _mm512_add_ps(accum_2, x02);
x02 = _mm512_load_ps(&x1[ 32]);
x02 = _mm512_loadu_ps(&x1[ 32]);
accum_3 = _mm512_add_ps(accum_3, x03);
x03 = _mm512_load_ps(&x1[ 48]);
x03 = _mm512_loadu_ps(&x1[ 48]);
accum_0 = _mm512_add_ps(accum_0, x04);
x04 = _mm512_load_ps(&x1[ 64]);
x04 = _mm512_loadu_ps(&x1[ 64]);
accum_1 = _mm512_add_ps(accum_1, x05);
x05 = _mm512_load_ps(&x1[ 80]);
x05 = _mm512_loadu_ps(&x1[ 80]);
accum_2 = _mm512_add_ps(accum_2, x06);
x06 = _mm512_load_ps(&x1[ 96]);
x06 = _mm512_loadu_ps(&x1[ 96]);
accum_3 = _mm512_add_ps(accum_3, x07);
x07 = _mm512_load_ps(&x1[112]);
x07 = _mm512_loadu_ps(&x1[112]);
n2 -= 128;
x1 += 128;
@@ -239,10 +233,10 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
accum_3 = _mm512_add_ps(accum_3, x07);
if (n2 >= 64) {
x00 = _mm512_load_ps(&x1[ 0]);
x01 = _mm512_load_ps(&x1[16]);
x02 = _mm512_load_ps(&x1[32]);
x03 = _mm512_load_ps(&x1[48]);
x00 = _mm512_loadu_ps(&x1[ 0]);
x01 = _mm512_loadu_ps(&x1[16]);
x02 = _mm512_loadu_ps(&x1[32]);
x03 = _mm512_loadu_ps(&x1[48]);
accum_0 = _mm512_add_ps(accum_0, x00);
accum_1 = _mm512_add_ps(accum_1, x01);
accum_2 = _mm512_add_ps(accum_2, x02);
@@ -253,8 +247,8 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 32) {
x00 = _mm512_load_ps(&x1[ 0]);
x01 = _mm512_load_ps(&x1[16]);
x00 = _mm512_loadu_ps(&x1[ 0]);
x01 = _mm512_loadu_ps(&x1[16]);
accum_0 = _mm512_add_ps(accum_0, x00);
accum_1 = _mm512_add_ps(accum_1, x01);
@@ -263,7 +257,7 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 16) {
x00 = _mm512_load_ps(&x1[ 0]);
x00 = _mm512_loadu_ps(&x1[ 0]);
accum_0 = _mm512_add_ps(accum_0, x00);
n2 -= 16;
@@ -272,7 +266,7 @@ static FLOAT casum_kernel(BLASLONG n, FLOAT *x)
if (n2) {
uint16_t tail_mask16 = (((uint16_t) 0xffff) >> (16 - n2));
x00 = _mm512_maskz_load_ps(*((__mmask16*) &tail_mask16), &x1[ 0]);
x00 = _mm512_maskz_loadu_ps(*((__mmask16*) &tail_mask16), &x1[ 0]);
accum_0 = _mm512_add_ps(accum_0, x00);
}
kernel/x86_64/dasum_microk_skylakex-2.c
+9 -16
View File
@@ -19,21 +19,14 @@ static FLOAT dasum_kernel(BLASLONG n, FLOAT *x1)
BLASLONG i = 0;
FLOAT sumf = 0.0;
if (n >= 256) {
BLASLONG align_512 = ((64 - ((uintptr_t)x1 & (uintptr_t)0x3f)) >> 3) & 0x7;
for (i = 0; i < align_512; i++) {
sumf += ABS_K(x1[i]);
}
n -= align_512;
x1 += align_512;
}
BLASLONG tail_index_SSE = n&(~7);
BLASLONG tail_index_AVX512 = n&(~255);
//
/* Deliberately no peeling to an alignment boundary: peeling makes the
grouping of the sum into accumulators - and hence the rounding of the
result - depend on the buffer address. Unaligned loads keep the result
a function of n alone, and cost the same as aligned loads when the
address happens to be aligned. */
if ( n >= 256 ) {
__m512d accum_0, accum_1, accum_2, accum_3;
@@ -42,10 +35,10 @@ static FLOAT dasum_kernel(BLASLONG n, FLOAT *x1)
accum_2 = _mm512_setzero_pd();
accum_3 = _mm512_setzero_pd();
for (i = 0; i < tail_index_AVX512; i += 32) {
accum_0 += _mm512_abs_pd(_mm512_load_pd(&x1[i + 0]));
accum_1 += _mm512_abs_pd(_mm512_load_pd(&x1[i + 8]));
accum_2 += _mm512_abs_pd(_mm512_load_pd(&x1[i +16]));
accum_3 += _mm512_abs_pd(_mm512_load_pd(&x1[i +24]));
accum_0 += _mm512_abs_pd(_mm512_loadu_pd(&x1[i + 0]));
accum_1 += _mm512_abs_pd(_mm512_loadu_pd(&x1[i + 8]));
accum_2 += _mm512_abs_pd(_mm512_loadu_pd(&x1[i +16]));
accum_3 += _mm512_abs_pd(_mm512_loadu_pd(&x1[i +24]));
}
accum_0 = accum_0 + accum_1 + accum_2 + accum_3;
kernel/x86_64/sasum_microk_skylakex-2.c
+9 -14
View File
@@ -18,19 +18,14 @@ static FLOAT sasum_kernel(BLASLONG n, FLOAT *x1)
BLASLONG i = 0;
FLOAT sumf = 0.0;
if (n >= 256) {
BLASLONG align_512 = ((64 - ((uintptr_t)x1 & (uintptr_t)0x3f)) >> 2) & 0xf;
for (i = 0; i < align_512; i++) {
sumf += ABS_K(x1[i]);
}
n -= align_512;
x1 += align_512;
}
BLASLONG tail_index_SSE = n&(~7);
BLASLONG tail_index_AVX512 = n&(~255);
/* Deliberately no peeling to an alignment boundary: peeling makes the
grouping of the sum into accumulators - and hence the rounding of the
result - depend on the buffer address. Unaligned loads keep the result
a function of n alone, and cost the same as aligned loads when the
address happens to be aligned. */
if (n >= 256) {
__m512 accum_0, accum_1, accum_2, accum_3;
accum_0 = _mm512_setzero_ps();
@@ -39,10 +34,10 @@ static FLOAT sasum_kernel(BLASLONG n, FLOAT *x1)
accum_3 = _mm512_setzero_ps();
for (i = 0; i < tail_index_AVX512; i += 64) {
accum_0 += _mm512_abs_ps(_mm512_load_ps(&x1[i + 0]));
accum_1 += _mm512_abs_ps(_mm512_load_ps(&x1[i +16]));
accum_2 += _mm512_abs_ps(_mm512_load_ps(&x1[i +32]));
accum_3 += _mm512_abs_ps(_mm512_load_ps(&x1[i +48]));
accum_0 += _mm512_abs_ps(_mm512_loadu_ps(&x1[i + 0]));
accum_1 += _mm512_abs_ps(_mm512_loadu_ps(&x1[i +16]));
accum_2 += _mm512_abs_ps(_mm512_loadu_ps(&x1[i +32]));
accum_3 += _mm512_abs_ps(_mm512_loadu_ps(&x1[i +48]));
}
accum_0 = accum_0 + accum_1 + accum_2 + accum_3;
kernel/x86_64/zasum_microk_skylakex-2.c
+31 -38
View File
@@ -122,7 +122,7 @@ static FLOAT zasum_kernel(BLASLONG n, FLOAT *x)
accum_2 = _mm512_setzero_pd();
accum_3 = _mm512_setzero_pd();
// alignment has side-effect when the size of input array is not large enough
// shorter unrolled path for small inputs
if (n2 < 128) {
if (n2 >= 64) {
x00 = _mm512_loadu_pd(&x1[ 0]);
@@ -208,29 +208,22 @@ static FLOAT zasum_kernel(BLASLONG n, FLOAT *x)
accum_0 = _mm512_add_pd(accum_0, accum_2);
sumf = _mm512_reduce_add_pd(accum_0);
}
// n2 >= 128, doing alignment
/* n2 >= 128. Deliberately no peeling to an alignment boundary:
peeling makes the grouping of the sum into accumulators - and
hence the rounding of the result - depend on the buffer address.
Unaligned loads keep the result a function of n alone, and cost
the same as aligned loads when the address happens to be
aligned. */
else {
int align_header = ((64 - ((uintptr_t)x1 & (uintptr_t)0x3f)) >> 3) & 0x7;
if (0 != align_header) {
unsigned char align_mask8 = (((unsigned char)0xff) >> (8 - align_header));
x00 = _mm512_maskz_loadu_pd(*((__mmask8*) &align_mask8), &x1[0]);
x00 = _mm512_and_pd(x00, abs_mask);
accum_0 = _mm512_add_pd(accum_0, x00);
n2 -= align_header;
x1 += align_header;
}
x00 = _mm512_load_pd(&x1[ 0]);
x01 = _mm512_load_pd(&x1[ 8]);
x02 = _mm512_load_pd(&x1[16]);
x03 = _mm512_load_pd(&x1[24]);
x04 = _mm512_load_pd(&x1[32]);
x05 = _mm512_load_pd(&x1[40]);
x06 = _mm512_load_pd(&x1[48]);
x07 = _mm512_load_pd(&x1[56]);
x00 = _mm512_loadu_pd(&x1[ 0]);
x01 = _mm512_loadu_pd(&x1[ 8]);
x02 = _mm512_loadu_pd(&x1[16]);
x03 = _mm512_loadu_pd(&x1[24]);
x04 = _mm512_loadu_pd(&x1[32]);
x05 = _mm512_loadu_pd(&x1[40]);
x06 = _mm512_loadu_pd(&x1[48]);
x07 = _mm512_loadu_pd(&x1[56]);
n2 -= 64;
x1 += 64;
@@ -241,26 +234,26 @@ static FLOAT zasum_kernel(BLASLONG n, FLOAT *x)
x02 = _mm512_and_pd(x02, abs_mask);
x03 = _mm512_and_pd(x03, abs_mask);
accum_0 = _mm512_add_pd(accum_0, x00);
x00 = _mm512_load_pd(&x1[ 0]);
x00 = _mm512_loadu_pd(&x1[ 0]);
accum_1 = _mm512_add_pd(accum_1, x01);
x01 = _mm512_load_pd(&x1[ 8]);
x01 = _mm512_loadu_pd(&x1[ 8]);
accum_2 = _mm512_add_pd(accum_2, x02);
x02 = _mm512_load_pd(&x1[16]);
x02 = _mm512_loadu_pd(&x1[16]);
accum_3 = _mm512_add_pd(accum_3, x03);
x03 = _mm512_load_pd(&x1[24]);
x03 = _mm512_loadu_pd(&x1[24]);
x04 = _mm512_and_pd(x04, abs_mask);
x05 = _mm512_and_pd(x05, abs_mask);
x06 = _mm512_and_pd(x06, abs_mask);
x07 = _mm512_and_pd(x07, abs_mask);
accum_0 = _mm512_add_pd(accum_0, x04);
x04 = _mm512_load_pd(&x1[32]);
x04 = _mm512_loadu_pd(&x1[32]);
accum_1 = _mm512_add_pd(accum_1, x05);
x05 = _mm512_load_pd(&x1[40]);
x05 = _mm512_loadu_pd(&x1[40]);
accum_2 = _mm512_add_pd(accum_2, x06);
x06 = _mm512_load_pd(&x1[48]);
x06 = _mm512_loadu_pd(&x1[48]);
accum_3 = _mm512_add_pd(accum_3, x07);
x07 = _mm512_load_pd(&x1[56]);
x07 = _mm512_loadu_pd(&x1[56]);
n2 -= 64;
x1 += 64;
@@ -286,10 +279,10 @@ static FLOAT zasum_kernel(BLASLONG n, FLOAT *x)
accum_3 = _mm512_add_pd(accum_3, x07);
if (n2 >= 32) {
x00 = _mm512_load_pd(&x1[ 0]);
x01 = _mm512_load_pd(&x1[ 8]);
x02 = _mm512_load_pd(&x1[16]);
x03 = _mm512_load_pd(&x1[24]);
x00 = _mm512_loadu_pd(&x1[ 0]);
x01 = _mm512_loadu_pd(&x1[ 8]);
x02 = _mm512_loadu_pd(&x1[16]);
x03 = _mm512_loadu_pd(&x1[24]);
x00 = _mm512_and_pd(x00, abs_mask);
x01 = _mm512_and_pd(x01, abs_mask);
x02 = _mm512_and_pd(x02, abs_mask);
@@ -304,8 +297,8 @@ static FLOAT zasum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 16) {
x00 = _mm512_load_pd(&x1[ 0]);
x01 = _mm512_load_pd(&x1[ 8]);
x00 = _mm512_loadu_pd(&x1[ 0]);
x01 = _mm512_loadu_pd(&x1[ 8]);
x00 = _mm512_and_pd(x00, abs_mask);
x01 = _mm512_and_pd(x01, abs_mask);
accum_0 = _mm512_add_pd(accum_0, x00);
@@ -316,7 +309,7 @@ static FLOAT zasum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 8) {
x00 = _mm512_load_pd(&x1[ 0]);
x00 = _mm512_loadu_pd(&x1[ 0]);
x00 = _mm512_and_pd(x00, abs_mask);
accum_0 = _mm512_add_pd(accum_0, x00);
@@ -326,7 +319,7 @@ static FLOAT zasum_kernel(BLASLONG n, FLOAT *x)
if (n2) {
unsigned char tail_mask8 = (((unsigned char) 0xff) >> (8 - n2));
x00 = _mm512_maskz_load_pd(*((__mmask8*) &tail_mask8), &x1[ 0]);
x00 = _mm512_maskz_loadu_pd(*((__mmask8*) &tail_mask8), &x1[ 0]);
x00 = _mm512_and_pd(x00, abs_mask);
accum_0 = _mm512_add_pd(accum_0, x00);
}
kernel/x86_64/zsum_microk_skylakex-2.c
+31 -37
View File
@@ -102,7 +102,7 @@ static FLOAT zsum_kernel(BLASLONG n, FLOAT *x)
accum_2 = _mm512_setzero_pd();
accum_3 = _mm512_setzero_pd();
// alignment has side-effect when the size of input array is not large enough
// shorter unrolled path for small inputs
if (n2 < 128) {
if (n2 >= 64) {
x00 = _mm512_loadu_pd(&x1[ 0]);
@@ -170,50 +170,44 @@ static FLOAT zsum_kernel(BLASLONG n, FLOAT *x)
accum_0 = _mm512_add_pd(accum_0, accum_2);
sumf = _mm512_reduce_add_pd(accum_0);
}
// n2 >= 128, doing alignment
/* n2 >= 128. Deliberately no peeling to an alignment boundary:
peeling makes the grouping of the sum into accumulators - and
hence the rounding of the result - depend on the buffer address.
Unaligned loads keep the result a function of n alone, and cost
the same as aligned loads when the address happens to be
aligned. */
else {
int align_header = ((64 - ((uintptr_t)x1 & (uintptr_t)0x3f)) >> 3) & 0x7;
if (0 != align_header) {
unsigned char align_mask8 = (((unsigned char)0xff) >> (8 - align_header));
x00 = _mm512_maskz_loadu_pd(*((__mmask8*) &align_mask8), &x1[0]);
accum_0 = _mm512_add_pd(accum_0, x00);
n2 -= align_header;
x1 += align_header;
}
x00 = _mm512_load_pd(&x1[ 0]);
x01 = _mm512_load_pd(&x1[ 8]);
x02 = _mm512_load_pd(&x1[16]);
x03 = _mm512_load_pd(&x1[24]);
x04 = _mm512_load_pd(&x1[32]);
x05 = _mm512_load_pd(&x1[40]);
x06 = _mm512_load_pd(&x1[48]);
x07 = _mm512_load_pd(&x1[56]);
x00 = _mm512_loadu_pd(&x1[ 0]);
x01 = _mm512_loadu_pd(&x1[ 8]);
x02 = _mm512_loadu_pd(&x1[16]);
x03 = _mm512_loadu_pd(&x1[24]);
x04 = _mm512_loadu_pd(&x1[32]);
x05 = _mm512_loadu_pd(&x1[40]);
x06 = _mm512_loadu_pd(&x1[48]);
x07 = _mm512_loadu_pd(&x1[56]);
n2 -= 64;
x1 += 64;
while (n2 >= 64) {
accum_0 = _mm512_add_pd(accum_0, x00);
x00 = _mm512_load_pd(&x1[ 0]);
x00 = _mm512_loadu_pd(&x1[ 0]);
accum_1 = _mm512_add_pd(accum_1, x01);
x01 = _mm512_load_pd(&x1[ 8]);
x01 = _mm512_loadu_pd(&x1[ 8]);
accum_2 = _mm512_add_pd(accum_2, x02);
x02 = _mm512_load_pd(&x1[16]);
x02 = _mm512_loadu_pd(&x1[16]);
accum_3 = _mm512_add_pd(accum_3, x03);
x03 = _mm512_load_pd(&x1[24]);
x03 = _mm512_loadu_pd(&x1[24]);
accum_0 = _mm512_add_pd(accum_0, x04);
x04 = _mm512_load_pd(&x1[32]);
x04 = _mm512_loadu_pd(&x1[32]);
accum_1 = _mm512_add_pd(accum_1, x05);
x05 = _mm512_load_pd(&x1[40]);
x05 = _mm512_loadu_pd(&x1[40]);
accum_2 = _mm512_add_pd(accum_2, x06);
x06 = _mm512_load_pd(&x1[48]);
x06 = _mm512_loadu_pd(&x1[48]);
accum_3 = _mm512_add_pd(accum_3, x07);
x07 = _mm512_load_pd(&x1[56]);
x07 = _mm512_loadu_pd(&x1[56]);
n2 -= 64;
x1 += 64;
@@ -230,10 +224,10 @@ static FLOAT zsum_kernel(BLASLONG n, FLOAT *x)
accum_3 = _mm512_add_pd(accum_3, x07);
if (n2 >= 32) {
x00 = _mm512_load_pd(&x1[ 0]);
x01 = _mm512_load_pd(&x1[ 8]);
x02 = _mm512_load_pd(&x1[16]);
x03 = _mm512_load_pd(&x1[24]);
x00 = _mm512_loadu_pd(&x1[ 0]);
x01 = _mm512_loadu_pd(&x1[ 8]);
x02 = _mm512_loadu_pd(&x1[16]);
x03 = _mm512_loadu_pd(&x1[24]);
accum_0 = _mm512_add_pd(accum_0, x00);
accum_1 = _mm512_add_pd(accum_1, x01);
accum_2 = _mm512_add_pd(accum_2, x02);
@@ -244,8 +238,8 @@ static FLOAT zsum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 16) {
x00 = _mm512_load_pd(&x1[ 0]);
x01 = _mm512_load_pd(&x1[ 8]);
x00 = _mm512_loadu_pd(&x1[ 0]);
x01 = _mm512_loadu_pd(&x1[ 8]);
accum_0 = _mm512_add_pd(accum_0, x00);
accum_1 = _mm512_add_pd(accum_1, x01);
@@ -254,7 +248,7 @@ static FLOAT zsum_kernel(BLASLONG n, FLOAT *x)
}
if (n2 >= 8) {
x00 = _mm512_load_pd(&x1[ 0]);
x00 = _mm512_loadu_pd(&x1[ 0]);
accum_0 = _mm512_add_pd(accum_0, x00);
n2 -= 8;
@@ -263,7 +257,7 @@ static FLOAT zsum_kernel(BLASLONG n, FLOAT *x)
if (n2) {
unsigned char tail_mask8 = (((unsigned char) 0xff) >> (8 - n2));
x00 = _mm512_maskz_load_pd(*((__mmask8*) &tail_mask8), &x1[ 0]);
x00 = _mm512_maskz_loadu_pd(*((__mmask8*) &tail_mask8), &x1[ 0]);
accum_0 = _mm512_add_pd(accum_0, x00);
}