simd и умножение матриц

#include <iostream>
#include <immintrin.h>
#include <chrono>
using namespace std;
 
const int M = 1024, N = 1024;
 
float scalar_mul(float** a, float** b, float** c);
float vector_mul(float** a, float** b, float** c);
float print_m(float** matrix);
float set_zero(float** matrix);
 
 
int main()
{
    setlocale(LC_ALL, "Russian");   
 
    float** a = new float* [M];
    float** b = new float* [M];
    float** c = new float* [M];
 
    for (int i = 0; i < M; i++) {
        a[i] = new float[N];
        b[i] = new float[N];
        c[i] = new float[N];
    }
 
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            a[i][j] = j;
            b[i][j] = j;
            c[i][j] = 0;
        }
    }
 
    printf("---Скалярное умножение матриц %dx%d---\n", M, N);
    auto start1 = chrono::high_resolution_clock::now();
    scalar_mul(a, b, c);
    auto end1 = chrono::high_resolution_clock::now();
    //print_m(c);
    chrono::duration<float> duration1 = end1 - start1;
    printf("Время выполнения: %f секунд\n\n", duration1.count());
 
    set_zero(c);
 
    printf("---Векторное умножение матриц %dx%d---\n", M, N);
    auto start2 = chrono::high_resolution_clock::now();
    vector_mul(a, b, c);
    auto end2 = chrono::high_resolution_clock::now();
    //print_m(c);
    chrono::duration<float> duration2 = end2 - start2;
    printf("Время выполнения: %f секунд\n\n", duration2.count());
 
    printf("Done!\n\n");
}
 
 
float scalar_mul(float** a, float** b, float** c) {
    /*  СКАЛЯРНОЕ УМНОЖЕНИЕ   */
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < M; k++) {
                c[i][j] = c[i][j] + a[i][k] * b[k][j];
            }
        }
    }
    return **c;
}
 
 
float vector_mul(float** a, float** b, float** c) {
    /*  ВЕКТОРНОЕ УМНОЖЕНИЕ   */
    //  транспонирование матрицы B
    float** bt = new float* [M];
 
    for (int i = 0; i < M; i++) {
        bt[i] = new float[N];
    }
 
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            bt[i][j] = b[j][i];
        }
    }
    //  умножение матриц
    __m256 c_line = _mm256_setzero_ps();
    int step = sizeof(c_line) / sizeof(a[0][0]);
 
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            __m256 c_line = _mm256_setzero_ps();
            for (int k = 0; k < N; k += step) {
 
                __m256 a_line = _mm256_load_ps(&a[i][k]);
 
                __m256 b_line = _mm256_load_ps(&bt[j][k]);
                __m256 tmp_line = _mm256_mul_ps(a_line, b_line);
                c_line = _mm256_add_ps(tmp_line, c_line);
            }
 
            __m256 t1 = _mm256_hadd_ps(c_line, c_line);
            __m256 t2 = _mm256_hadd_ps(t1, t1);
            __m128 t3 = _mm256_extractf128_ps(t2, 1);
            __m128 t4 = _mm_add_ss(_mm256_castps256_ps128(t2), t3);
 
            _mm_store_ps(&c[i][j], t4);
        }
    }
    return **c;
}
 
 
float print_m(float** matrix) {
    /*  ВЫВОД МАТРИЦЫ НА ЭКРАН   */
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
           cout << matrix[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
    return 0;
}
 
 
float set_zero(float** matrix) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            matrix[i][j] = 0;
        }
    }
    return **matrix;
}