ArhSSE

#include <iostream>
#include <ctime>
#include <xmmintrin.h>
#include <cstdlib>
using namespace std;
#define p 8000
float **matrA = (float **)_mm_malloc(p * sizeof(float), 16);
float **matrB = (float **)_mm_malloc(p * sizeof(float), 16);
float *vectX = (float *)_mm_malloc(p * sizeof(float), 16);
float *vectY = (float *)_mm_malloc(p * sizeof(float), 16);
void standart_function(float**matrA, float **matrB, float*vectX, float*vectY, float a, float b, int n);
void sse_function(float**matrA, float **matrB, float*vectX, float*vectY, float a, float b, int n);
float a, b;
int main()
{
    //srand(time(0));
    setlocale(LC_ALL, "rus");
    clock_t time;
    bool exitbool = true;
    a = (rand() % 200 - 100) + ((rand() % 100) / 100.);
    b = (rand() % 200 - 100) + ((rand() % 100) / 100.);
    int step = 1000, highborder = 8000, tempnum;
    for (int i = 0; i < p; i++)
    {
        matrA[i] = (float *)_mm_malloc(p * sizeof(float), 16);
        matrB[i] = (float *)_mm_malloc(p * sizeof(float), 16);
        vectX[i] = (rand() % 200 - 100) + ((rand() % 100) / 100.);
        vectY[i] = (rand() % 200 - 100) + ((rand() % 100) / 100.);
    }
    for (int i = 0; i < p; i++)
    {
        for (int j = 0; j < p; j++)
        {
            matrA[i][j] = (rand() % 200 - 100) + ((rand() % 100) / 100.);
            matrB[i][j] = (rand() % 200 - 100) + ((rand() % 100) / 100.);
        }
    }
    if (highborder > 0 && highborder < p)
    {
        cout << "Введите размер массива: ";
        cin >> highborder;
    }
    if (step > 0 && step < p&&step <= highborder)
    {
        cout << "Введите шаг(массив будет отсортирован от нуля с заданным шагом): ";
        cin >> step;
    }
    while (exitbool)
    {
        cout << "Выберите способ выполнения операции alpha*matrix(A)*vector(x)+beta*matrix(B)*vector(b):\n1)выпонить стандартными операциями;\n2)выполнить с использованием расширений SSE;\nИначе, введите \"0\": ";
        cin >> tempnum;
        cout << endl;

        for (int j = 0; j <= highborder; j += step)
        {
            switch (tempnum)
            {
            case 1: {
                time = clock();
                standart_function(matrA, matrB, vectX, vectY, a, b, j);
                cout << "Стандартными операциями (" << j << "): " << ((double)(clock() - time)) / CLOCKS_PER_SEC << " секунд" << endl;
                break;
            }
            case 2: {
                time = clock();
                sse_function(matrA, matrB, vectX, vectY, a, b, j);
                cout << "SSE (" << j << "): " << ((double)(clock() - time)) / CLOCKS_PER_SEC << " секунд" << endl;
                break; }
            case 0: {exitbool = false; }
            }
        }
        cout << endl;
    }
        for (int i = 0; i < p; i++)
        {
            _mm_free(matrA[i]);
            _mm_free(matrB[i]);
        }
        _mm_free(vectX);
        _mm_free(vectY);
        _mm_free(matrA);
        _mm_free(matrB);


    return 0;
}

void standart_function(float**matrA, float **matrB, float*vectX, float*vectY, float a, float b, int n)
{
    float *rdvect1, *rdvect2;

    rdvect1 = (float *)_mm_malloc(n * sizeof(float), 16);//выделяем память
    rdvect2 = (float *)_mm_malloc(n * sizeof(float), 16);
    for (int i = 0; i < n; i++)
    {
        rdvect1[i] = 0;
        for (int j = 0; j < n; j++)
        {
            rdvect1[i] += matrA[i][j] * vectX[j];
            rdvect2[i] += matrB[i][j] * vectY[j];
        }

        rdvect1[i] = a*rdvect1[i] + b*rdvect2[i];
    }
}

void sse_function(float**matrA, float **matrB, float*vectX, float*vectY, float a, float b, int n)
{
    __m128 *xx, *yy, ssrdvect1, sstemp1, ssrdvect2, sstemp2, aa, bb;
    xx = (__m128 *)vectX;
    yy = (__m128 *)vectY;

    aa = _mm_set_ps1(a);
    bb = _mm_set_ps1(b);

    float *rdvect11 = (float *)_mm_malloc(n * sizeof(float), 16);//выделяем память
    float *rdvect22 = (float *)_mm_malloc(n * sizeof(float), 16);
    float *ready_vector = (float *)_mm_malloc(n * sizeof(float), 16);
    for (int i = 0; i < n; i++)
    {

        sstemp1 = _mm_set_ps1(a);//устанавливаем 4 значения, выведенных по адресу
        sstemp2 = _mm_set_ps1(b);

        for (int j = 0; j<n / 4; j++)
        {
            ssrdvect1 = _mm_mul_ps(*((__m128 *)*(matrA + i) + j), *(xx + j));// векторное умножение четырех чисел
            sstemp1 = _mm_add_ps(sstemp1, _mm_add_ps(ssrdvect1, _mm_movehl_ps(ssrdvect1, ssrdvect1)));

            ssrdvect2 = _mm_mul_ps(*((__m128 *)*(matrB + i) + j), *(yy + j));
            ssrdvect2 = _mm_add_ps(sstemp2, _mm_add_ps(ssrdvect2, _mm_movehl_ps(ssrdvect2, ssrdvect2)));
        }
        _mm_store_ss((rdvect11 + i), _mm_add_ps(sstemp1, _mm_shuffle_ps(sstemp1, sstemp1, 1)));//запись младшего значения в память
        _mm_store_ss((rdvect22 + i), _mm_add_ps(sstemp2, _mm_shuffle_ps(sstemp2, sstemp2, 1)));
    }

    xx = (__m128 *)rdvect11;
    yy = (__m128 *)rdvect22;
    for (int i = 0; i < n / 4; i++, xx++, yy++)
        _mm_store_ps((ready_vector + i * 4), _mm_add_ps(_mm_mul_ps(aa, *xx), _mm_mul_ps(bb, *yy)));

    _mm_free(rdvect11);//освобождаем память
    _mm_free(rdvect22);
    _mm_free(ready_vector);

}