Segment fault using Intel SIMD, even the space is very large and is of multiple of 32 bytes

Question

I keep getting segment fault when using SIMD instructions to optimize matrix multiplication.

Here is the core computing part. The matrices are stored like this: a large vector<double> buf of size (3 * 1025 * 1025) is allocated. matrix A starts from buf[0], matrix B starts from buf[1025] and C starts from buf[1025*2]. I performed various matrix multiplication with size of 4 to 1024. So they could all fit in this vector.

#include <immintrin.h>
#define BLOCK_SIZE 4
/*
 * performs 4 * 4 matrix multiplication C=A*B
 * C is 4-by-4, A is 4-by-4, and B is 4-by-4, column major matrices
 * lda is the size of the large matrix.
 */
static void do_block(int lda4, double* A, double* B, double* C) {
    int n=4;
    for(int i=0; i<n; i++){ // process i th column
      for(int j=0; j<n; j++){
        __m256d c = _mm256_load_pd(C+j*lda);
        c = _mm256_fmadd_pd(_mm256_load_pd(A+i*lda), _mm256_broadcast_sd(B+i+j*lda), c);
        _mm256_store_pd(C+j*lda, c);
      }
    }
}

/* This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values. */
void square_dgemm(int lda, double* A, double* B, double* C) {
    for (int j = 0; j < lda; j += BLOCK_SIZE) {
        // Accumulate block dgemms into block of C
        for (int k = 0; k < lda; k += BLOCK_SIZE) {
            // For each block-row of A
            for (int i = 0; i < lda; i += BLOCK_SIZE) {
                do_block(lda, A + i + k * lda, B + k + j * lda, C + i + j * lda);
            }
        }
    }
}

The weird thing is: When I change the vector from size of (3 * 1025 * 1025) to (3 * 1024 * 1024), it gives me segment fault.

My questions are:

1. I have learnt that these instructions require aligned data. Indeed replacing by unaligned versions like _mm256_loadu_pd eliminates this error. However, since size of (3 * 1024 * 1024 * sizeof(double)) % 32 bytes == 0, isn't it 32 bytes-aligned or I misunderstood the concept?
2. I have allocated very large contiguous space, and why it crashes from the beginning, when performing small mat mul (4*4)? I thought that as long as I calling _mm256_load_pd(addr) with at least 32 bytes allocated starting from addr, it won't crash, am I wrong?
3. Why it doesn't crash on buf of (3 * 1025 * 1025), but crashes on (3 * 1024 * 1024) ? It seems doesn't crashes when the size is an odd number, like 1025, 1027, 1029 and always crash when the number is even, like 1024, 1026.

The code was compiled using GCC, with -march=native and -O3. The CPU supports FMA, AVX and AVX2. The machine is Google Cloud VM, the CPU is Intel Xeon which I cannot get the exact model. Thanks for your advice!

Answer 1

Thanks for the comments and I think I have figured out what was wrong.

Answer to Q1: Yes, I misunderstood the concept of 'Alignment'. I have to ensure that the adress is the multiple of 32. So after manually (not an elegant way by doing pointer math) making 3 matrices starting from 36 bytes aligned space, it works.

Answer to Q2: Since the adress passed to _mm256_load_pd(addr) is not 32 bytes aligned, even if it's very large, an exception might be raised, according to Intel documentation

Answer to Q3: This is quite tricky: Not both of A and C are 32 bytes aligned. However, changing the total size will cause C to be 32 bytes aligned (But A still not aligned), and it doesn't crash. This works just by accident and it's not safe, as A is not aligned.

So vectorization of matrix multiplication is trciky, I need to

1. Ensure matrix A and C to be 32 bytes aligned, maybe use something like aligned_alloc mentioned in the comments.
2. pad the matrices into multiple of 4, to ensure in every step calling _mm256_load_pd(addr) the adress is properly aligned.