1. 程式人生 > >c/c++ 程式碼中使用sse指令集加速

c/c++ 程式碼中使用sse指令集加速

使用SSE指令,首先要了解這一類用於進行初始化載入資料以及將暫存器的資料儲存到記憶體相關的指令,

我們知道,大多數SSE指令是使用的xmm0到xmm8的暫存器,那麼使用之前,就需要將資料從記憶體載入到這些暫存器。

1. load系列,用於載入資料,從記憶體到暫存器

__m128 _mm_load_ss (float *p)  
__m128 _mm_load_ps (float *p)  
__m128 _mm_load1_ps (float *p)  
__m128 _mm_loadh_pi (__m128 a, __m64 *p)  
__m128 _mm_loadl_pi (__m128 a, __m64 *p)  
__m128 _mm_loadr_ps (float *p)  
__m128 _mm_loadu_ps (float *p)

上面是從手冊查詢到的load系列的函式。其中,

_mm_load_ss用於scalar的載入,所以,載入一個單精度浮點數到暫存器的低位元組,其它三個位元組清0,(r0 := *p, r1 := r2 := r3 := 0.0)。

_mm_load_ps用於packed的載入(下面的都是用於packed的),要求p的地址是16位元組對齊,否則讀取的結果會出錯,(r0 := p[0], r1 := p[1], r2 := p[2], r3 := p[3])。

_mm_load1_ps表示將p地址的值,載入到暫存器的四個位元組,需要多條指令完成,所以,從效能考慮,在內層迴圈不要使用這類指令。(r0 := r1 := r2 := r3 := *p)。

_mm_loadh_pi和_mm_loadl_pi分別用於從兩個引數高底位元組等組合載入。具體參考手冊。

_mm_loadr_ps表示以_mm_load_ps反向的順序載入,需要多條指令完成,當然,也要求地址是16位元組對齊。(r0 := p[3], r1 := p[2], r2 := p[1], r3 := p[0])。

_mm_loadu_ps和_mm_load_ps一樣的載入,但是不要求地址是16位元組對齊,對應指令為movups。

2. set系列,用於載入資料,大部分需要多條指令完成,但是可能不需要16位元組對齊。

__m128 _mm_set_ss (float w)  
__m128 _mm_set_ps (float z, float y, float x, float w)  
__m128 _mm_set1_ps (float w)  
__m128 _mm_setr_ps (float z, float y, float x, float w)  
__m128 _mm_setzero_ps ()  

這一系列函式主要是類似於load的操作,但是可能會呼叫多條指令去完成,方便的是可能不需要考慮對齊的問題。

_mm_set_ss對應於_mm_load_ss的功能,不需要位元組對齊,需要多條指令。(r0 = w, r1 = r2 = r3 = 0.0)

_mm_set_ps對應於_mm_load_ps的功能,引數是四個單獨的單精度浮點數,所以也不需要位元組對齊,需要多條指令。(r0=w, r1 = x, r2 = y, r3 = z,注意順序)

_mm_set1_ps對應於_mm_load1_ps的功能,不需要位元組對齊,需要多條指令。(r0 = r1 = r2 = r3 = w)

_mm_setzero_ps是清0操作,只需要一條指令。(r0 = r1 = r2 = r3 = 0.0)

3. store系列,用於將計算結果等SSE暫存器的資料儲存到記憶體中。

void _mm_store_ss (float *p, __m128 a)  
void _mm_store_ps (float *p, __m128 a)  
void _mm_store1_ps (float *p, __m128 a)  
void _mm_storeh_pi (__m64 *p, __m128 a)  
void _mm_storel_pi (__m64 *p, __m128 a)  
void _mm_storer_ps (float *p, __m128 a)  
void _mm_storeu_ps (float *p, __m128 a)  
void _mm_stream_ps (float *p, __m128 a)

這一系列函式和load系列函式的功能對應,基本上都是一個反向的過程。

_mm_store_ss:一條指令,*p = a0

_mm_store_ps:一條指令,p[i] = a[i]。

_mm_store1_ps:多條指令,p[i] = a0。

_mm_storeh_pi,_mm_storel_pi:值儲存其高位或低位。

_mm_storer_ps:反向,多條指令。

_mm_storeu_ps:一條指令,p[i] = a[i],不要求16位元組對齊。

_mm_stream_ps:直接寫入記憶體,不改變cache的資料。

(2)算術指令

SSE提供了大量的浮點運算指令,包括加法、減法、乘法、除法、開方、最大值、最小值、近似求倒數、求開方的倒數等等,可見SSE指令的強大之處。那麼在瞭解了上面的資料載入和資料儲存的指令之後,使用這些算術指令就很容易了,下面以加法為例。

SSE中浮點加法的指令有:

__m128 _mm_add_ss (__m128 a, __m128 b)  
__m128 _mm_add_ps (__m128 a, __m128 b) 

其中,_mm_add_ss表示scalar執行模式,_mm_add_ps表示packed執行模式。

一般而言,使用SSE指令寫程式碼,步驟為:使用load/set函式將資料從記憶體載入到SSE暫存器;使用相關SSE指令完成計算等;使用store系列函式將結果從暫存器儲存到記憶體,供後面使用。

下面是一個完成加法的例子:

#include <intrin.h>  
  
int main(int argc, char* argv[])  
{  
    float op1[4] = {1.0, 2.0, 3.0, 4.0};  
    float op2[4] = {1.0, 2.0, 3.0, 4.0};  
    float result[4];  
  
    __m128  a;  
    __m128  b;  
    __m128  c;  
  
    // Load  
    a = _mm_loadu_ps(op1);  
    b = _mm_loadu_ps(op2);  
  
    // Calculate  
    c = _mm_add_ps(a, b);   // c = a + b  
  
    // Store  
    _mm_storeu_ps(result, c);  
  
    /*      // Using the __m128 union to get the result. 
    printf("0: %lf\n", c.m128_f32[0]); 
    printf("1: %lf\n", c.m128_f32[1]); 
    printf("2: %lf\n", c.m128_f32[2]); 
    printf("3: %lf\n", c.m128_f32[3]); 
    */  
    printf("0: %lf\n", result[0]);  
    printf("1: %lf\n", result[1]);  
    printf("2: %lf\n", result[2]);  
    printf("3: %lf\n", result[3]);  
  
    return 0;  
}  

這個例子,前面已經寫過類似的加法例子,註釋中的printf部分是利用__m128這個資料型別來獲取相關的值,

這個型別是一個union型別,具體定義可以參考相關標頭檔案,但是,對於實際使用,有時候這個值是一箇中間值,需要後面計算使用,就得使用store了,效率更高。

上面使用的是_mm_loadu_ps和_mm_storeu_ps,不要求位元組對齊,如果使用_mm_load_ps和_mm_store_ps,會發現程式會崩潰或得不到正確結果。下面是指定位元組對齊後的一種實現方法:

#include <intrin.h>  
  
int main(int argc, char* argv[])  
{  
    __declspec(align(16)) float op1[4] = {1.0, 2.0, 3.0, 4.0};  
    __declspec(align(16)) float op2[4] = {1.0, 2.0, 3.0, 4.0};  
    _MM_ALIGN16 float result[4];        // A macro, same as "__declspec(align(16))"   
  
    __m128  a;  
    __m128  b;  
    __m128  c;  
  
    // Load  
    a = _mm_load_ps(op1);  
    b = _mm_load_ps(op2);  
  
    // Calculate  
    c = _mm_add_ps(a, b);   // c = a + b  
  
    // Store  
    _mm_store_ps(result, c);  
  
    /*      // Using the __m128 union to get the result. 
    printf("0: %lf\n", c.m128_f32[0]); 
    printf("1: %lf\n", c.m128_f32[1]); 
    printf("2: %lf\n", c.m128_f32[2]); 
    printf("3: %lf\n", c.m128_f32[3]); 
    */  
    printf("0: %lf\n", result[0]);  
    printf("1: %lf\n", result[1]);  
    printf("2: %lf\n", result[2]);  
    printf("3: %lf\n", result[3]);  
  
    return 0;  
}  

下面是矩陣乘以向量的程式碼,以及sse版

vec3 operator*(const vec3 &v) const {
        vec3 ret;
        ret[0] = mat[0] * v[0] + mat[4] * v[1] + mat[8] * v[2] + mat[12];
        ret[1] = mat[1] * v[0] + mat[5] * v[1] + mat[9] * v[2] + mat[13];
        ret[2] = mat[2] * v[0] + mat[6] * v[1] + mat[10] * v[2] + mat[14];
        return ret;
    }

    vec3 MultVec3SSE(const vec3& v) const 
    {
        vec3 ret;
        
        
        __m128   res1;
        __m128   res2;
        __m128   res3;
        __m128   addres;

        /*
        _MM_ALIGN16 float addresf[4] = {mat[12],mat[13],mat[14],1};

        _MM_ALIGN16 float multf1[4] = {mat[0],mat[1],mat[2],1};
        _MM_ALIGN16 float multf2[4] = {mat[4],mat[5],mat[6],1};
        _MM_ALIGN16 float multf3[4] = {mat[8],mat[9],mat[10],1};

        _MM_ALIGN16 float mvecf1[4] = {v[0],v[0],v[0],1};
        _MM_ALIGN16 float mvecf2[4] = {v[1],v[1],v[1],1};
        _MM_ALIGN16 float mvecf3[4] = {v[2],v[2],v[2],1};
        
        addres = _mm_load_ps(addresf);
        */
        addres = _mm_set_ps(1,mat[14],mat[13],mat[12]);
/*
        
        res1 = _mm_mul_ps(_mm_load_ps(multf1),_mm_load_ps(mvecf1));
        res2 = _mm_mul_ps(_mm_load_ps(multf2),_mm_load_ps(mvecf2));
        res3 = _mm_mul_ps(_mm_load_ps(multf3),_mm_load_ps(mvecf3));  //_mm_set_ps
*/        

        res1 = _mm_mul_ps(_mm_set_ps(1,mat[2],mat[1],mat[0]),_mm_set_ps(1,v[0],v[0],v[0]));
        res2 = _mm_mul_ps(_mm_set_ps(1,mat[6],mat[5],mat[4]),_mm_set_ps(1,v[1],v[1],v[1]));
        res3 = _mm_mul_ps(_mm_set_ps(1,mat[10],mat[9],mat[8]),_mm_set_ps(1,v[2],v[2],v[2]));

        res1 = _mm_add_ps(res1,res2);
        res3 = _mm_add_ps(res1,res3);
        
        res3 = _mm_add_ps(res3,addres);
        float dest[4];
        _mm_storeu_ps(dest,res3);
        
        ret.x = dest[0] ;
        ret.y = dest[1] ;
        ret.z = dest[2] ;

        return ret;
    }

debug模式下fps只提高了個位數,release下  更是跑不過 非sse的。。。。。情何以堪

SSE 加速運算例子詳解:乘法、加法、平方、最小值、最大值、與操作

庫檔案說明

#ifndef __METHOD
#define __METHOD


void ScaleValue1(float *pArray, DWORD dwCount, float fScale);//乘法
void ScaleValue2(float *pArray, DWORD dwCount, float fScale);

void Add1(float *pArray, DWORD dwCount, float fScale);//加法
void Add2(float *pArray, DWORD dwCount, float fScale);

void Sqrt1(float *pArray, DWORD dwCount, float fScale);//平方
void Sqrt2(float *pArray, DWORD dwCount, float fScale);

void Min1(float *pArray, DWORD dwCount, float fScale);//最小值
void Min2(float *pArray, DWORD dwCount, float fScale);//最小值

void Max1(float *pArray, DWORD dwCount, float fScale);//最小值
void Max2(float *pArray, DWORD dwCount, float fScale);//最小值

void And1(float *pArray, DWORD dwCount, float fScale);//與操作
void And2(float *pArray, DWORD dwCount, float fScale);//與操作

#endif
#include <xmmintrin.h>
#include <Windows.h>
#include <math.h>

void ScaleValue1(float *pArray, DWORD dwCount, float fScale)//乘法
{
    DWORD dwGroupCount = dwCount/4;
    __m128 e_Scale = _mm_set_ps1(fScale);//設定所有4個值為同一值

    for (DWORD i=0; i<dwGroupCount; i++)
    {
        *(__m128*)(pArray + i*4) = _mm_mul_ps( *(__m128*)(pArray + i*4),e_Scale);
    }
}

void ScaleValue2(float *pArray, DWORD dwCount, float fScale)
{
    for (DWORD i =0; i<dwCount; i++)
    {
        pArray[i] *= fScale;
    }
}

void Add1(float *pArray, DWORD dwCount, float fScale)//加法
{
    DWORD dwGroupCount = dwCount/4;
    __m128 e_Scale = _mm_set_ps1(fScale);//設定所有4個值為同一值

    for (DWORD i=0; i<dwGroupCount; i++)
    {
        *(__m128*)(pArray + i*4) = _mm_add_ps( *(__m128*)(pArray + i*4),e_Scale);
    }
}

void Add2(float *pArray, DWORD dwCount, float fScale)
{
    for (DWORD i =0; i<dwCount; i++)
    {
        pArray[i] += fScale;
    }
}

void Sqrt1(float *pArray, DWORD dwCount, float fScale)//平方
{
    DWORD dwGroupCount = dwCount/4;
    __m128 e_Scale = _mm_set_ps1(fScale);//設定所有4個值為同一值

    for (DWORD i=0; i<dwGroupCount; i++)
    {
        *(__m128*)(pArray + i*4) = _mm_sqrt_ps(e_Scale);
    }
}

void Sqrt2(float *pArray, DWORD dwCount, float fScale)
{
    for (DWORD i =0; i<dwCount; i++)
    {
        pArray[i] = sqrt(fScale);
    }
}

void Min1(float *pArray, DWORD dwCount, float fScale)//最小值
{
    DWORD dwGroupCount = dwCount/4;
    __m128 e_Scale = _mm_set_ps1(fScale);//設定所有4個值為同一值

    for (DWORD i=0; i<dwGroupCount; i++)
    {
        *(__m128*)(pArray + i*4) = _mm_min_ps( *(__m128*)(pArray + i*4),e_Scale);
    }
}

void Min2(float *pArray, DWORD dwCount, float fScale)
{
    for (DWORD i =0; i<dwCount; i++)
    {
        pArray[i] = (pArray[i]>fScale? fScale : pArray[i]);
    }
}

void Max1(float *pArray, DWORD dwCount, float fScale)//最大值
{
    DWORD dwGroupCount = dwCount/4;
    __m128 e_Scale = _mm_set_ps1(fScale);//設定所有4個值為同一值

    for (DWORD i=0; i<dwGroupCount; i++)
    {
        *(__m128*)(pArray + i*4) = _mm_max_ps( *(__m128*)(pArray + i*4),e_Scale);
    }
}

void Max2(float *pArray, DWORD dwCount, float fScale)
{
    for (DWORD i =0; i<dwCount; i++)
    {
        pArray[i] = (pArray[i]<fScale? fScale : pArray[i]);
    }
}

void And1(float *pArray, DWORD dwCount, float fScale)//與操作
{
    DWORD dwGroupCount = dwCount/4;
    __m128 e_Scale = _mm_set_ps1(fScale);//設定所有4個值為同一值

    for (DWORD i=0; i<dwGroupCount; i++)
    {
        *(__m128*)(pArray + i*4) = _mm_and_ps( *(__m128*)(pArray + i*4),e_Scale);
    }
}

void And2(float *pArray, DWORD dwCount, float fScale)
{
    for (DWORD i =0; i<dwCount; i++)
    {
        pArray[i] = (int)(pArray[i]) & (int)(fScale);
    }
}

SSE指令集加速運算

先上程式碼:

/*g++ -msse2 main.cpp -lrt*/

#include <iostream>
#include <xmmintrin.h>//SSE指令集需包含詞標頭檔案
#include <time.h>
using namespace std;

#define N 120
int main() {

    struct timespec tpstart,tpend;
    clock_gettime(CLOCK_MONOTONIC, &tpstart);
    /////////////////do something
        {

        __m128 *p1,*p2,*p3;//__m128是一個長128位的資料型別,存放在暫存器中
        __attribute(aligned(128)) float f1[N],f2[N],f3[N];//新建一些浮點型陣列
        /*__attribute(aligned(128))強制使編譯器在給f1等分配
        記憶體空間時對齊在128位(8位元組)上*/
        //cout<<f1<<"t"<<(int)f1%128<<endl;
        /*執行此條語句可以看出地址f1總是128的倍數,在記憶體上即
        是總存在一行上,使得匯流排取值時能夠一次取完*/
        for(int i=0; i<N; i++) {
            f1[i]=i+0.12;
            f2[i]=i+0.16;
            }
        for(int time=0; time<10000; time++) {
            p1=(__m128*)f1;
            p2=(__m128*)f2;
            p3=(__m128*)f3;
            for(int i=0; i<N; i+=4) {
                *p3=_mm_mul_ps(*p1,*p2);//此函式封裝了一些彙編原語,使得可以同時計算4個值的乘法
                p3++;
                p2++;
                p1++;
                }
            /*注意N的取值,N取為4的整數倍,使得最後一次調
            用SSE指令時剛好能全部處理完。否則容易引發段錯誤。
            所以,即使不能使用完,也一定要申請4的整數倍空間*/
            }
        }
    /////////////////done
    clock_gettime(CLOCK_MONOTONIC, &tpend);
    long timedif = (tpend.tv_sec-tpstart.tv_sec)*1000*1000+(tpend.tv_nsec-tpstart.tv_nsec)/1000;
    cout<<"SSE:t"<<timedif<<endl;

/////////////////////////不採用SSE指令/////////////////////////////////////


    clock_gettime(CLOCK_MONOTONIC, &tpstart);
    /////////////////do something
        {

        float f1[N],f2[N],f3[N];
        for(int i=0; i<N; i++) {
            f1[i]=i+0.12;
            f2[i]=i+0.16;
            }

        for(int time=0; time<10000; time++) {
            for(int i=0; i<N; i++) {
                f3[i]=f1[i]*f2[i];
                }
            }
        }
    /////////////////done
    clock_gettime(CLOCK_MONOTONIC, &tpend);
    timedif = (tpend.tv_sec-tpstart.tv_sec)*1000*1000+(tpend.tv_nsec-tpstart.tv_nsec)/1000;
    cout<<"noSSE:t"<<timedif<<endl;

    return 0;
    }

注意計時函式需要包含的標頭檔案和編譯選項,clock_gettime函式能精確到納秒。

執行結果為:

SSE:    4929
noSSE:    9980

以微秒為單位,平均提高約兩倍。

下面是另一份程式碼,演示如何動態申請空間,注意記憶體對齊:

#include <xmmintrin.h>
#include <iostream>
using namespace std;
const int N=120;

int main()
{
    __m128 *p1,*p2,*p3;
    __attribute(aligned(128)) float *pf1,*pf2,*pf3;
    pf1=(float *)_mm_malloc(sizeof(float)*N,128);
    pf2=(float *)_mm_malloc(sizeof(float)*N,128);
    pf3=(float *)_mm_malloc(sizeof(float)*N,128);
    cout<<(int)pf1%128<<endl;
    cout<<(int)pf2%128<<endl;
    cout<<(int)pf3%128<<endl;

    for(int i=0;i<N;i++){
        pf1[i]=i+0.12;
        pf2[i]=i+0.16;
    };
    p1=(__m128*)pf1;
    p2=(__m128*)pf2;
    p3=(__m128*)pf3;
    for(int i=0;i<N;i+=4){
        *p3=_mm_add_ps(*p1,*p2);
        p3++;
        p2++;
        p1++;
    }
    _mm_free(pf1);
    _mm_free(pf2);
    _mm_free(pf3);
    return 0;
}