memset((void*)y, 0, n * sizeof(double));
    double* ptr_a = a;
    double* ptr_x = NULL;
    const int vectorSize = 4;
    const int outerStep = 4;
    const int innerStep = 16;
    __m256d ra0, ra1, ra2, ra3;
    __m256d rx0, rx1, rx2, rx3;
    __m256d ry0, ry1, ry2, ry3;
    __declspec(align(32)) double buf0[vectorSize];
    __declspec(align(32)) double buf1[vectorSize];
    __declspec(align(32)) double buf2[vectorSize];
    __declspec(align(32)) double buf3[vectorSize];

    for (int i = 0; i < n; i += outerStep, ptr_a += (outerStep - 1) * n) {
        ry0 = ry1 = ry2 = ry3 = _mm256_setzero_pd();
        ptr_x = x;
        for (int j = 0; j < n; j += innerStep, ptr_a += innerStep, ptr_x += innerStep) {
            _mm_prefetch((char*)(ptr_x + innerStep), _MM_HINT_T0);
            _mm_prefetch((char*)(ptr_x + innerStep + 8), _MM_HINT_T0);
            _mm_prefetch((char*)(ptr_a + innerStep), _MM_HINT_NTA);
            _mm_prefetch((char*)(ptr_a + innerStep + 8), _MM_HINT_NTA);
            _mm_prefetch((char*)(ptr_a + n + innerStep), _MM_HINT_NTA);
            _mm_prefetch((char*)(ptr_a + n + innerStep + 8), _MM_HINT_NTA);
            _mm_prefetch((char*)(ptr_a + 2 * n + innerStep), _MM_HINT_NTA);
            _mm_prefetch((char*)(ptr_a + 2 * n + innerStep + 8), _MM_HINT_NTA);
            _mm_prefetch((char*)(ptr_a + 3 * n + innerStep), _MM_HINT_NTA);
            _mm_prefetch((char*)(ptr_a + 3 * n + innerStep + 8), _MM_HINT_NTA);
            rx0 = _mm256_load_pd(ptr_x);
            rx1 = _mm256_load_pd(ptr_x + vectorSize);
            rx2 = _mm256_load_pd(ptr_x + 2 * vectorSize);
            rx3 = _mm256_load_pd(ptr_x + 3 * vectorSize);
            ra0 = _mm256_load_pd(ptr_a);
            ra1 = _mm256_load_pd(ptr_a + n);
            ra2 = _mm256_load_pd(ptr_a + 2 * n);
            ra3 = _mm256_load_pd(ptr_a + 3 * n);
            ry0 = _mm256_fmadd_pd(ra0, rx0, ry0);
            ry1 = _mm256_fmadd_pd(ra1, rx0, ry1);
            ry2 = _mm256_fmadd_pd(ra2, rx0, ry2);
            ry3 = _mm256_fmadd_pd(ra3, rx0, ry3);
            ra0 = _mm256_load_pd(ptr_a + vectorSize);
            ra1 = _mm256_load_pd(ptr_a + n + vectorSize);
            ra2 = _mm256_load_pd(ptr_a + 2 * n + vectorSize);
            ra3 = _mm256_load_pd(ptr_a + 3 * n + vectorSize);
            ry0 = _mm256_fmadd_pd(ra0, rx1, ry0);
            ry1 = _mm256_fmadd_pd(ra1, rx1, ry1);
            ry2 = _mm256_fmadd_pd(ra2, rx1, ry2);
            ry3 = _mm256_fmadd_pd(ra3, rx1, ry3);
            ra0 = _mm256_load_pd(ptr_a + 2 * vectorSize);
            ra1 = _mm256_load_pd(ptr_a + n + 2 * vectorSize);
            ra2 = _mm256_load_pd(ptr_a + 2 * n + 2 * vectorSize);
            ra3 = _mm256_load_pd(ptr_a + 3 * n + 2 * vectorSize);
            ry0 = _mm256_fmadd_pd(ra0, rx2, ry0);
            ry1 = _mm256_fmadd_pd(ra1, rx2, ry1);
            ry2 = _mm256_fmadd_pd(ra2, rx2, ry2);
            ry3 = _mm256_fmadd_pd(ra3, rx2, ry3);
            ra0 = _mm256_load_pd(ptr_a + 3 * vectorSize);
            ra1 = _mm256_load_pd(ptr_a + n + 3 * vectorSize);
            ra2 = _mm256_load_pd(ptr_a + 2 * n + 3 * vectorSize);
            ra3 = _mm256_load_pd(ptr_a + 3 * n + 3 * vectorSize);
            ry0 = _mm256_fmadd_pd(ra0, rx3, ry0);
            ry1 = _mm256_fmadd_pd(ra1, rx3, ry1);
            ry2 = _mm256_fmadd_pd(ra2, rx3, ry2);
            ry3 = _mm256_fmadd_pd(ra3, rx3, ry3);
        }
        _mm256_store_pd(buf0, ry0);
        _mm256_store_pd(buf1, ry1);
        _mm256_store_pd(buf2, ry2);
        _mm256_store_pd(buf3, ry3);
        y[i] = buf0[0] + buf0[1] + buf0[2] + buf0[3];
        y[i + 1] = buf1[0] + buf1[1] + buf1[2] + buf1[3];
        y[i + 2] = buf2[0] + buf2[1] + buf2[2] + buf2[3];
        y[i + 3] = buf3[0] + buf3[1] + buf3[2] + buf3[3];
    }
    ptr_a = ptr_x = NULL;