C++矩阵乘法

我想使用寄存器（逐行中的信息）通过矢量算法创建矩阵乘法。打开外环4次我有问题void matvec_XMM（双*a、双*x、双*y、int n、int lb）返回错误结果的函数这是我必须使用的算法：

i = 1,n,4
   r0 = r1 = r2 = r3 = 0
   j = 1,n,8
      r0 = r0 + aij * xj + ai,j+1 * xj+1 + … + ai,j+7 * xj+7
      r1 = r1 + ai+1,j * xj + ai+1,j+1 * xj+1 + … + ai+1,j+7 * xj+7
      r2 = r2 + ai+2,j * xj + ai+2,j+1 * xj+1 + … + ai+2,j+7 * xj+7
      r3 = r3 + ai+3,j * xj + ai+3,j+1 * xj+1 + … + ai+3,j+7 * xj+7
   end j
   yi = r0; yi+1 = r1; yi+2 = r2; yi+3 = r3; 
end i

它是ma代码：

#include "stdafx.h"
#include <iostream>
#include "mvec.h"
#include <emmintrin.h>
using namespace std;
void mult_naive(double *a, double *x, double *y, int n)
{
    int i, j, ij;
    double register reg;
    for(i=0, ij=0; i<n; ++i)
    {
        reg = 0;
        for(j=0; j<n; ++j, ++ij)
        {
            reg += a[ij]*x[j];
        }
        y[i] = reg;
    }
}
void matvec_XMM(double* a, double* x, double* y, int n, int lb)
{
int i, j;
memset((void *)y, 0, n*sizeof(double));
double res0[2];
double res1[2];
double res2[2];
double res3[2];
__m128d ry0, ry1, ry2, ry3, ra0, rx0;
double *ptr_a, *ptr_x, *ptr_y;
const int nr = 4;
ptr_a = a;
for (i = 0; i < n; i+=nr)   
{
    ry0 = _mm_setzero_pd();
    ry1 = _mm_setzero_pd();
    ry2 = _mm_setzero_pd();
    ry3 = _mm_setzero_pd(); 
    ptr_y = &y[i];
    ptr_x = x;
    for (j = 0; j<n; j+=lb)
    {
        _mm_prefetch((const char *)(ptr_a + lb*nr), _MM_HINT_NTA); 
        _mm_prefetch((const char *)(ptr_x + lb), _MM_HINT_T0);
        //----1
        rx0 = _mm_load_pd(ptr_x);       
        ra0 = _mm_load_pd(ptr_a);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry0 = _mm_add_pd(ry0, ra0);
        ra0 = _mm_load_pd(ptr_a + 2); 
        ra0 = _mm_mul_pd(ra0, rx0);
        ry1 = _mm_add_pd(ry1, ra0);
        ra0 = _mm_load_pd(ptr_a + 4);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry2 = _mm_add_pd(ry2, ra0);
        ra0 = _mm_load_pd(ptr_a + 6);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry3 = _mm_add_pd(ry3, ra0);
            //----2
        rx0 = _mm_load_pd(ptr_x + 2);   
        ra0 = _mm_load_pd(ptr_a + 8);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry0 = _mm_add_pd(ry0, ra0);
        ra0 = _mm_load_pd(ptr_a + 10);  
        ra0 = _mm_mul_pd(ra0, rx0);
        ry1 = _mm_add_pd(ry1, ra0);
        ra0 = _mm_load_pd(ptr_a + 12);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry2 = _mm_add_pd(ry2, ra0);
        ra0 = _mm_load_pd(ptr_a + 14);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry3 = _mm_add_pd(ry3, ra0);
        //----3
        rx0 = _mm_load_pd(ptr_x + 4);       
        ra0 = _mm_load_pd(ptr_a + 16);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry0 = _mm_add_pd(ry0, ra0);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry1 = _mm_add_pd(ry1, ra0);
        ra0 = _mm_load_pd(ptr_a + 20);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry2 = _mm_add_pd(ry2, ra0);
        ra0 = _mm_load_pd(ptr_a + 22);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry3 = _mm_add_pd(ry3, ra0);
        //----4
        rx0 = _mm_load_pd(ptr_x + 6);       
        ra0 = _mm_load_pd(ptr_a + 24);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry0 = _mm_add_pd(ry0, ra0);
        ra0 = _mm_load_pd(ptr_a + 26);  
        ra0 = _mm_mul_pd(ra0, rx0);
        ry1 = _mm_add_pd(ry1, ra0);
        ra0 = _mm_load_pd(ptr_a + 28);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry2 = _mm_add_pd(ry2, ra0);
        ra0 = _mm_load_pd(ptr_a + 30);
        ra0 = _mm_mul_pd(ra0, rx0);
        ry3 = _mm_add_pd(ry3, ra0);
        ptr_a += lb*nr;
        ptr_x += lb;
    }
    _mm_store_pd(res0, ry0);
    *ptr_y = res0[0] + res0[1];
    _mm_store_pd(res1, ry1);
    *(ptr_y + 1) = res1[0] + res1[1];
    _mm_store_pd(res2, ry2);
    *(ptr_y + 2)= res2[0] + res2[1];
    _mm_store_pd(res3, ry3);
    *(ptr_y + 3) = res3[0] + res3[1];
}
}

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include "windows.h"
#include "mvec.h"
using namespace std;

int main(int argc, char* argv[])
{
    double *a, *x, *y, *z;
    int n;
    DWORD tstart;
    const int lb = 8;
    double elaps_time;
    cout << "Program Mat_Vect: performance y = y +A*xn";
#ifdef _DEBUG
    cout << "DEBUG versionn";
#else
    cout << "RELEASE versionn";
#endif
    cout << "Input dimensionn";
    cin >> n;
    n = n/lb;
    n = lb*n;
    try
    {
        a = new double [n*n];
        x = new double [n+1];
        y = new double [n];
        z = new double [n];
    }
    catch(bad_alloc aa)
    {
        cout << "memory allocation error" << endl;
        system("pause");
        exit(1);
    }
    memset((void *)a, 0, _msize((void *)a));
    memset((void *)x, 0, _msize((void *)x));
    memset((void *)y, 0, _msize((void *)y));
    cout << "startn";
    prepare(a, x, n);
    //-------------------------naive algorithm-----------------------//
    cout << "naive algorithm: n";
    tstart = GetTickCount();
    mult_naive(a, x, z, n);
    elaps_time = (double)(GetTickCount()-tstart)/1000.0;
    cout << "naive algorithm: " << elaps_time << " sec" << endl;
    //-------------------------algorithm which uses XMM registers-----------------------//
    delete [] a;
    delete [] x;
    a = (double *)_aligned_malloc(n*n*sizeof(double), 16);
    x = (double *)_aligned_malloc(n*sizeof(double), 16);
    if(!a || !x)
    {
        cout << "memory allocation error" << endl;
        system("pause");
        exit(1);
    }
    cout << "algorithm which uses XMM: n";
    prepare(a, x, n);
    tstart = GetTickCount();
    matvec_XMM(a, x, y, n, lb);
    elaps_time = (double)(GetTickCount()-tstart)/1000.0;
    check(y, z, n);
    cout << "algorithm which uses XMM: " << elaps_time << " sec" << endl;
    delete [] y;
    delete [] z;
    _aligned_free(a);
    _aligned_free(x);
    system("pause");
    return 0;
}

void check(double *y, double *z, int n)
{
    int i;
    for(i=0; i<n; i++)
    {
        if(fabs(z[i] - y[i]) > 1.0e-9)
        {
            cout << "errorn";
            return;
        }
    }
    cout << "OKn";
}
void prepare(double *a, double *x, int n)
{
    int i, j, ij;
    for(i=0, ij=0; i<n; i++)
    {
        for(j=0; j<n; j++, ij++)
        {
            if(i == j)
                a[ij] = 10.0;
            else
                a[ij] = (double)(i+1);
        }
        x[i] = 1.0;
    }
}