用SSE做向量內積

SSE是Streaming SIMD Extensions的縮寫，SIMD是Single Instruction Multiple Data的縮寫，即對多條數據並行執行相同的操作，以提高CPU的性能。

從SSE4.1開始支持向量的內積：__m128 _mm_dp_ps( __m128 a,__m128 b,const int mask)。更老的版本需要綜合運用多條指令才可以辦到。_m128表示128位寄存器，可以存儲4個單精度浮點數，2個雙精度浮點數，在SSE指令集中ps代表單精度接口，pd代表雙精度接口。dp是dot production的意思。

如何查看自己的機器是否支持SSE4.1 ？

$ cat /proc/cpuinfo  | grep sse

如果你能看到“sse4_1”或更高的版本就可以了。

不多解釋，直接上代碼。

#include <smmintrin.h>       //SSE
#include <iostream>
#include <ctime>


float inner_product(const float* x, const float* y, const long & len){
    float prod = 0.0f;
    long i;
    for (i=0;i<len;i++){
        prod+=x[i]*y[i];
    }
    return prod;
}

float dot_product(const float* x, const float* y, const long & len){
    float prod = 0.0f;
    const int mask = 0xff;        //把每一位都設為1。詳見https://msdn.microsoft.com/en-us/library/bb514054(v=vs.120).aspx
    __m128 X, Y;
    float tmp;


    long i;
    for (i=0;i<len;i+=4){
        X=_mm_loadu_ps(x+i); //_mm_loadu_ps把float轉為__m128
        Y=_mm_loadu_ps(y+i);
        _mm_storeu_ps(&tmp,_mm_dp_ps(X,Y,mask));//_mm_storeu_ps把__m128轉為float。_mm_dp_ps計算向量內積（向量的長度不能超過4），在https://software.intel.com/sites/default/files/m/8/b/8/D9156103.pdf第3頁最下方對_mm_dp_ps有說明 
        prod += tmp;
    }

    return prod;
}

int main(){
    const int len1 = 10;
    float arr[len1]={2.0f,5.0f,3.0f,1.0f,7.0f,9.0f,4.0f,3.0f,7.0f,5.0f};
    float brr[len1]={9.0f,0.0f,3.0f,5.0f,3.0f,7.0f,8.0f,2.0f,6.0f,1.0f};
    std::cout<<"蠻力計算內積 "<<inner_product(arr,brr,len1)<<std::endl;
    std::cout<<"使用SSE計算內積 "<<dot_product(arr,brr,len1)<<std::endl;

    const int len2 = 1000000;
    float *crr=new float[len2];
    float *drr=new float[len2];
    for (int i=0;i<len2;i++){
        int value=i%10;
        crr[i]=value;
        drr[i]=value;
    }

    float prod;
    clock_t begin;
    clock_t end;
    begin=clock();
    prod=inner_product(crr,drr,len2);
    end=clock();
    std::cout<<"蠻力計算內積 "<<prod<<"\t用時"<<(double)(end-begin)/CLOCKS_PER_SEC<<"秒"<<std::endl;
    begin=clock();
    prod=dot_product(crr,drr,len2);
    end=clock();
    std::cout<<"使用SSE計算內積 "<<prod<<"\t用時"<<(double)(end-begin)/CLOCKS_PER_SEC<<"秒"<<std::endl;
}

//g++ -m32 -msse4.1 dp.cpp -o dp

python numpy中的向量內積調用了C語言並行庫（具體是什么並行計算庫我不知道），看一下numpy.dot()性能如何。

#coding=utf-8
__author__='orisun'

import numpy as np 
import time

LEN=1000000
arr=[]
brr=[]
for i in xrange(LEN):
    value=i%10
    arr.append(value)
    brr.append(value)
array1=np.array(arr)
array2=np.array(brr)

begin=time.time()
prod=0.0
for i in xrange(LEN):
    prod+=arr[i]*brr[i]
end=time.time()
print "蠻力計算內積 %.2e\t\t用時%f秒" % (prod,end-begin)

begin=time.time()
prod=np.dot(array1,array2)
end=time.time()
print "使用numpy計算內積 %.2e\t用時%f秒" % (prod,end-begin)

Go語言以C為原型，我們來領略下Go的威力。

package main

import (
	"fmt"
	"time"
)

func InnerProduct(x []float32, y []float32) float32{
	var rect float32=0
	for i:=0;i<len(x);i++{
		rect+=x[i]*y[i] 
	}
	return rect
}

func main() {
	const len2 int=1000000
	crr:=[]float32{}
	drr:=[]float32{}
	for i:=0;i<len2;i++{
		value:=float32(i%10)
		crr=append(crr,value)
		drr=append(drr,value)
	}

	begin:=time.Now().UnixNano()
	prod:=InnerProduct(crr,drr)
	end:=time.Now().UnixNano()
	fmt.Printf("蠻力計算內積 %.0f\t\t用時%f秒\n",prod,float32(end-begin)/1e9)
}

結果

結論：

1. 用SSE計算向量內積比用蠻力計算明顯要快，但快不到4倍（4倍應該是極限值），計算量越大SSE的優勢越明顯。
2. 用蠻力計算向量內積，python的性能遠不如C++，C++比python快了50多倍。但C++的float運算損失了0.7%的精度，python里的小數都是雙精度，沒有損失精度。
3. Go秒殺C++，Go的蠻力法比C++的並行法還要快！
4. numpy.dot()比用SSE的_mm_dp_ps還在快一點。猜測：如果numpy.dot()底層也是用的SSE，那它顯然不是直接調的_mm_dp_ps，_mm_dp_ps對簡單的指令進行了封裝並且引入了mask參數自然會慢一些。

參考

https://msdn.microsoft.com/en-us/library/bb514054(v=vs.120).aspx

https://software.intel.com/sites/default/files/m/8/b/8/D9156103.pdf

https://software.intel.com/sites/landingpage/IntrinsicsGuide/#techs=SSE4_1