对矩阵向量乘法进行SIMD优化

使用simd和循环展开技术，对矩阵向量乘法dgemv进行改写，以实现y=Ax计算的性能优化。原始程序如下所示：

#include <iostream>
#include <ctime>
#include <omp.h>
using namespace std;

void dgemv(const int n, const int m, const double *const A,
           const double *const x, double *y)
{
    for (int i = 0; i < n; i++)
    {
        double y_tmp = 0.0;
        for (int j = 0; j < m; j++)
        {
            y_tmp += A[i * m + j] * x[j];
        }
        y[i] = y_tmp;
    }
    return;
}

int main()
{
    int nloop = 100;
    int n = 10240;
    double *mat_A = new double[n * n];
    double *vec_x = new double[n];
    double *vec_y = new double[n];

    srand(0);
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            vec_x[j] = (double)(rand() % 10) + (double)rand() / RAND_MAX * 3.14;
            mat_A[i * n + j] = (double)(rand() % 10) + (double)rand() / RAND_MAX * 3.14;
        }
        vec_y[i] = 0.0;
    }

    clock_t t = clock();
    for (int iloop = 0; iloop < nloop; iloop++)
    {
        dgemv(n, n, mat_A, vec_x, vec_y);
    }
    t = clock() - t;

    cout << "y[0]:" << vec_y[0] << endl;
    cout << "y[1]:" << vec_y[1] << endl;
    cout << "y[2]:" << vec_y[2] << endl;
    cout << "y[3]:" << vec_y[3] << endl;
    cout << "y[4]:" << vec_y[4] << endl;

    cout << "gemv ave cost time(clock):" << t / nloop << endl;
}

运行结果如下：

使用SIMD优化后，代码如下

#include <iostream>
#include <ctime>
#include <arm_neon.h>
using namespace std;

void dgemv_neon(const int n, const int m, const double *const A, const double *const x, double *y)
{
    // #pragma omp simd
    for (int i = 0; i < n; i+=2)
    {
        float64x2_t y_sum0 = vdupq_n_f64(0.0);
        float64x2_t y_sum1 = vdupq_n_f64(0.0);

        for (int j = 0; j < m; j += 2)
        {
            float64x2_t v_x = vld1q_f64(&x[j]);
            y_sum0 = vfmaq_f64(y_sum0, vld1q_f64(&A[i * m + j]), v_x);
            y_sum1 = vfmaq_f64(y_sum1, vld1q_f64(&A[(i + 1) * m + j]), v_x);
        }
        y[i] = vaddvq_f64(y_sum0);
        y[i + 1] = vaddvq_f64(y_sum1);
    }
    return;
}

int main()
{
    int nloop = 100;
    int n = 10240;
    double *mat_A = new double[n * n];
    double *vec_x = new double[n];
    double *vec_y = new double[n];

    srand(0);
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            vec_x[j] = (double)(rand() % 10) + (double)rand() / RAND_MAX * 3.14;
            mat_A[i * n + j] = (double)(rand() % 10) + (double)rand() / RAND_MAX * 3.14;
        }
        vec_y[i] = 0.0;
    }

    clock_t t = clock();
    for (int iloop = 0; iloop < nloop; iloop++)
    {
        dgemv_neon(n, n, mat_A, vec_x, vec_y);
    }
    t = clock() - t;

    cout << "y[0]:" << vec_y[0] << endl;
    cout << "y[1]:" << vec_y[1] << endl;
    cout << "y[2]:" << vec_y[2] << endl;
    cout << "y[3]:" << vec_y[3] << endl;
    cout << "y[4]:" << vec_y[4] << endl;

    cout << "gemv ave cost time(clock):" << t / nloop << endl;
}

运行结果如下：