图片 1

SSE图像算法优化系列十一:使用FFT变换实现图像卷积。

在GPU上执行能够加快我的应用程序吗?

  

     
本文重点主要不在于FFT的SSE优化,而在于使用FFT实现快速卷积的相关技巧和过程。

GPU能够对符合以下标准的应用程序进行加速:

  这个月6号开始,着手解决一个具有实际意义的计算任务。任务数据有9879896条,每条包含30个整数,任务是计算每两条数据之间的斯皮尔相关系数及其P值。原始数据只有500+MB,因此我并不认为这是个多么大的计算任务。随后稍加计算,我还是很惊呆的,要计算(9879896×9879895)÷2≈4.88亿亿组数据,但此时这还只是个数字概念,我也没意识到时间复杂度和空间复杂度的问题。

      关于FFT变换,有很多参考的代码,特别是对于长度为2的整数次幂的序列,实现起来也是非常简易的,而对于非2次幂的序列,就稍微有点麻烦了,matlab中是可以实现任意长度FFT的,FFTW也是可以的,而Opencv则有选择性的实现了某些长度序列的变换,查看Opencv的代码,可以发现其只有对是4的整数次幂的数据部分采用了SSE优化,比如4、16、64、256、1024这样的序列部分,因此基4的FFT是最快的,而剩余的部分则依旧是普通的C语言实现,我这里主要是把Opencv的代码抠出来了。

大规模并行—计算能够被分割成上百个或上千个独立的工作单元。

 

     Opencv关于FFT实现的代码在Opencv
3.0\opencv\sources\modules\core\src\dxt.cpp中,代码写的特别复杂,扣取工作也做的相当艰苦,其基4的SSE优化核心代码如下所示:

计算密集型—计算消耗的时间显著超过了花费转移数据到GPU内存以及从GPU内存转移出数据的时间。

1. 计算规模初体验

// optimized radix-4 transform
template<> struct DFT_VecR4<float>
{
    int operator()(Complex<float>* dst, int N, int n0, int& _dw0, const Complex<float>* wave) const
    {
        int n = 1, i, j, nx, dw, dw0 = _dw0;
        __m128 z = _mm_setzero_ps(), x02=z, x13=z, w01=z, w23=z, y01, y23, t0, t1;
        Cv32suf t; t.i = 0x80000000;
        __m128 neg0_mask = _mm_load_ss(&t.f);
        __m128 neg3_mask = _mm_shuffle_ps(neg0_mask, neg0_mask, _MM_SHUFFLE(0,1,2,3));

        for( ; n*4 <= N; )
        {
            nx = n;
            n *= 4;
            dw0 /= 4;

            for( i = 0; i < n0; i += n )
            {
                Complexf *v0, *v1;

                v0 = dst + i;
                v1 = v0 + nx*2;

                x02 = _mm_loadl_pi(x02, (const __m64*)&v0[0]);
                x13 = _mm_loadl_pi(x13, (const __m64*)&v0[nx]);
                x02 = _mm_loadh_pi(x02, (const __m64*)&v1[0]);
                x13 = _mm_loadh_pi(x13, (const __m64*)&v1[nx]);

                y01 = _mm_add_ps(x02, x13);
                y23 = _mm_sub_ps(x02, x13);
                t1 = _mm_xor_ps(_mm_shuffle_ps(y01, y23, _MM_SHUFFLE(2,3,3,2)), neg3_mask);
                t0 = _mm_movelh_ps(y01, y23);
                y01 = _mm_add_ps(t0, t1);
                y23 = _mm_sub_ps(t0, t1);

                _mm_storel_pi((__m64*)&v0[0], y01);
                _mm_storeh_pi((__m64*)&v0[nx], y01);
                _mm_storel_pi((__m64*)&v1[0], y23);
                _mm_storeh_pi((__m64*)&v1[nx], y23);

                for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
                {
                    v0 = dst + i + j;
                    v1 = v0 + nx*2;

                    x13 = _mm_loadl_pi(x13, (const __m64*)&v0[nx]);
                    w23 = _mm_loadl_pi(w23, (const __m64*)&wave[dw*2]);
                    x13 = _mm_loadh_pi(x13, (const __m64*)&v1[nx]); // x1, x3 = r1 i1 r3 i3
                    w23 = _mm_loadh_pi(w23, (const __m64*)&wave[dw*3]); // w2, w3 = wr2 wi2 wr3 wi3

                    t0 = _mm_mul_ps(_mm_moveldup_ps(x13), w23);
                    t1 = _mm_mul_ps(_mm_movehdup_ps(x13), _mm_shuffle_ps(w23, w23, _MM_SHUFFLE(2,3,0,1)));
                    x13 = _mm_addsub_ps(t0, t1);
                    // re(x1*w2), im(x1*w2), re(x3*w3), im(x3*w3)
                    x02 = _mm_loadl_pi(x02, (const __m64*)&v1[0]); // x2 = r2 i2
                    w01 = _mm_loadl_pi(w01, (const __m64*)&wave[dw]); // w1 = wr1 wi1
                    x02 = _mm_shuffle_ps(x02, x02, _MM_SHUFFLE(0,0,1,1));
                    w01 = _mm_shuffle_ps(w01, w01, _MM_SHUFFLE(1,0,0,1));
                    x02 = _mm_mul_ps(x02, w01);
                    x02 = _mm_addsub_ps(x02, _mm_movelh_ps(x02, x02));
                    // re(x0) im(x0) re(x2*w1), im(x2*w1)
                    x02 = _mm_loadl_pi(x02, (const __m64*)&v0[0]);

                    y01 = _mm_add_ps(x02, x13);
                    y23 = _mm_sub_ps(x02, x13);
                    t1 = _mm_xor_ps(_mm_shuffle_ps(y01, y23, _MM_SHUFFLE(2,3,3,2)), neg3_mask);
                    t0 = _mm_movelh_ps(y01, y23);
                    y01 = _mm_add_ps(t0, t1);
                    y23 = _mm_sub_ps(t0, t1);

                    _mm_storel_pi((__m64*)&v0[0], y01);
                    _mm_storeh_pi((__m64*)&v0[nx], y01);
                    _mm_storel_pi((__m64*)&v1[0], y23);
                    _mm_storeh_pi((__m64*)&v1[nx], y23);
                }
            }
        }

        _dw0 = dw0;
        return n;
    }
};

不满足上述标准的应用程序在GPU上运行时可能会比CPU要慢。

数据格式:9879896行,30列,每列之间以空格符隔开,例如:

  其实也不复杂,总觉得里面的_mm_loadh_pi用的很好,这也得益于其数据源采用了Complex格式,这样数据的实部和虚部在内存中是连续的,用SSE加载数据是也就很方便了。

使用MATLAB进行GPU编程

0 2 0 2 0 0 0 0 0 0 0 40 0 0 35 0 0 53 0 44 0 0 0 0 0 0 0 0 0 0
0 0 1 148 0 0 0 0 0 0 0 0 0 0 1133 0 1 0 0 1820 0 0 0 2 0 0 0 1 0 0
0 0 0 33 1 0 0 0 0 0 0 0 0 0 231 0 0 0 0 402 0 0 0 0 0 0 0 0 0 0
0 0 6 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 6 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
... ...
... ...

  上面只是部分代码,为了配合该过程,还有很多的初始化,比如数据的shuffle等。详见Opencv的文件。

FFT,IFFT以及线性代数运算超过了100个内置的MATLAB函数,通过提供一个类型为GPUArray(由并行计算工具箱提供的特殊数组类型)的输入参数,这些函数就能够直接在GPU上运行。这些启用GPU的函数都是重载的,换句话说,这些函数根据传递的参数类型的不同而执行不同的操作。

空间复杂度:单纯计算下结果大概有多大吧,每组计算结果包含相关系数和P值,若都以float(占4字节)精度存储,需要占用内存:4.88亿亿×8B≈400TB,当然,我们不具备这么大内存,因此无论以何种方式计算,都需要一批批地重复将数据载入内存、计算、存入硬盘这个过程,直到运算完成。那么,存入硬盘的结果会有400TB吗?不然,P值小于或等于0.05的结果才会需要输出,因此实际上会远远小于这个值,具体会小多少,先运行一批数据后才能做出估算。

  对于2维的FFT变换,我没有去扣CV的代码,而是直接先每行进行一维的FFT1D,然后对结果在进行列方向的FFT1D,由于FFT1D算法需要处理的序列必须是连续的内存,因此,需要对中间的结果进行转置,处理完后在转置回来,为了节省时间,这个转置也应该用SSE优化。

例如,以下代码使用FFT算法查找CPU上伪随机数向量的离散傅里叶变换:

时间复杂度:计算的组数规模是(n×(n-1))÷2,那么就看程序能跑多快了。我想先看看MATLAB多线程、Python多线程、Spark分布式计算能跑多快,是否能在最快时间内解决问题。

  当2维的宽度和高度相同时,这个转置是不需要分配另外一份额外的内存的,这个叫In-Place转置,另外一个重要优点就是FFT1D算法也支持In-Place操作。

A = rand(2^16,1);

 

  由于是对Complex数据进行转置,我们可以换个角度考虑,因为一个Complex正好和double数据占用同样的内存,这样直接利用和double相关的SSE指令就可以方便的实现Complex相关数据的转置。比如当宽度和高度一样时,刻直接使用下述方式。

B = fft (A);

2. MATLAB多线程

inline void Inplace_TransposeSSE2X2D(double *Src, double *Dest, int Length)
{
    __m128d S0 = _mm_loadu_pd(Src);                        
    __m128d S1 = _mm_loadu_pd((Src + Length));            
    __m128d D0 = _mm_loadu_pd(Dest);                
    __m128d D1 = _mm_loadu_pd((Dest + Length));            
    _mm_storeu_pd(Dest, _mm_unpacklo_pd(S0, S1));
    _mm_storeu_pd(Dest + Length, _mm_unpackhi_pd(S0, S1));
    _mm_storeu_pd(Src, _mm_unpacklo_pd(D0, D1));
    _mm_storeu_pd(Src + Length, _mm_unpackhi_pd(D0, D1));
}

为在GPU上执行相同的操作,我们首先使用gpuArray命令将数据从MATLAB工作空间转移至GPU设备内存。然后我们能够运行重载函数fft:

  MATLAB写起来最简单,计算相关系数和P值都不用操心,一行自带的函数调用就完成。打开MATLAB左下角的并行池,MATLAB将会自动寻找到机子上有的物理核心,并分配与物理核心数相同的worker。比如我的电脑是4核8线程,它只能开4个worker,不识别虚拟核心。

  参考系统的Intrinsic的_MM_TRANSPOSE4_PS函数,可以发现double类型的转置方式不太一样,可以直接使用有关的unpack函数,而不使用shuffle,当然使用shuffle也是没有问题的。

A = gpuArray(rand(2^16,1));

  代码如下:

  那么最后我实现的FFT2D函数如下所示:

B = fft (A);

t1 = clock;
disp('>> loading ...');
A = importdata('D:/MASTER2016/5.CUDA/data-ID-top30-kv3.txt');
b = A'; %由于MATLAB只计算列与列之间的相关系数,因此需要转置操作
disp(etime(clock,t1));

num = size(b, 2);
disp('>> calculating ...');
fid = fopen('D:/MASTER2016/5.CUDA/result-matlab.txt', 'wt');

for i = 1 : num
    for j = i+1 : num
        [m, n] = corr(b(:, i), b(:, j), 'type', 'Spearman', 'tail', 'both');
        if isnan(n) || n>0.05
            continue;
        end

        fprintf(fid, 'X%d\tX%d\t%d\t%d\n', i, j, m, n);
    end
end
fclose(fid);
disp('>> OK!');
int FFT2D(Complex *Input, Complex *Output, int Width, int Height, bool Invert, int StartZeroLines = 0, int EndZeroLines = 0)
{
    if ((Input == NULL) || (Output == NULL))                                return IM_STATUS_NULLREFRENCE;
    if ((IsPowerOfTwo(Width) == false) || (IsPowerOfTwo(Height) == false))    return IM_STATUS_INVALIDPARAMETER;

    if (Width == Height)
    {
        void *Buffer = FFT_Init(Width);
        if (Buffer == NULL)    return IM_STATUS_OUTOFMEMORY;

        for (int Y = StartZeroLines; Y < Height - EndZeroLines; Y++)
            FFT1D(Input + Y * Width, Output + Y * Width, Buffer, Width, Invert);            //    先水平方向变换
        InPlace_TransposeD((double *)Output, Width);                                        //    再位转置

        for (int Y = 0; Y < Height; Y++)
            FFT1D(Output + Y * Width, Output + Y * Width, Buffer, Width, Invert);            //    再垂直方向变换
        InPlace_TransposeD((double *)Output, Width);
        FFT_Free(Buffer);
    }
    else
    {
        void *Buffer = FFT_Init(Width);
        if (Buffer == NULL)    return IM_STATUS_OUTOFMEMORY;

        for (int Y = StartZeroLines; Y < Height - EndZeroLines; Y++)
            FFT1D(Input + Y * Width, Output + Y * Width, Buffer, Width, Invert);
        Complex *Temp = (Complex *)malloc(Width * Height * sizeof(Complex));
        if (Temp == NULL)
        {
            FFT_Free(Buffer);
            return IM_STATUS_OUTOFMEMORY;
        }
        TransposeD((double *)Output, (double *)Temp, Width, Height);
        FFT_Free(Buffer);

        Buffer = FFT_Init(Height);                                                            //    必须要二次初始化
        if (Buffer == NULL)    return IM_STATUS_OUTOFMEMORY;

        for (int Y = 0; Y < Width; Y++)
            FFT1D(Temp + Y * Height, Temp + Y * Height, Buffer, Height, Invert);
        TransposeD((double *)Temp, (double *)Output, Height, Width);
        FFT_Free(Buffer);
        free(Temp);
    }
    return IM_STATUS_OK;
}

fft操作在GPU上而不是在CPU上执行,因为输入参数(GPUArray)位于GPU的内存中。

  这里我并没有考虑内存空间不够的问题,因为我只是想说明MATLAB的计算速度。开了多颗核心的情况下,MATLAB并没能完全压榨出所有的CPU性能,计算速度缓慢无比,更要命的是,它会越算越慢。据我估算,即使空间复杂度足够,MATLAB也要用超过20年的时间才能算完,这还是不考虑越算越慢的情况。

  注意在函数的最后我增加了StartZeroLines
和EndZeroLines
两个参数,他们主要是在进行行方向的FFT1D是忽略最前面的StartZeroLines和最后面的EndZeroLines
个全0的行,因为全0的变换后面的结果还是0,没有必要进行计算,减少计算时间。在opencv中也是有个nonzero_rows的,但是他只是针对前面的行,而实际中,比如很多情况是,先要扩展数据,然后把有效数据放置在左上角,这样只有右下角有非零元素,所以opencv这样做就无法起到加速作用了。

结果B存储在GPU当中。然而,B在MATLAB工作空间中依旧可见。通过运行class(B),我们看到B是一个GPUArray。

  好了,此方案仅是打酱油。

  对扣取的代码进行了实际的测试,1024*1024的数据,进行100次正反变换,耗时3000ms,使用matlab进行同样的操作,耗时约5500ms,并且观察任务管理器,在4核PC上CPU使用率100%,说明他内部使用了多线程,不过有一点就是matlab使用的是double类型的数据。听说matlab最新版使用的就是FFTW库,不过无论如何,这个速度还是可以接受和相当快的。

class(B)

 

  下面我们重点谈下基于FFT的图像卷积的实现,理论上如果图像a大小为N
* M,卷积核b大小为 X * Y,则卷积实现的过程如下:

ans =

3. Python多线程

  首先扩展数据,扩展后的大小为 (N + X –
1) * (M + Y –
1),将卷积核数据放置到扩展后的数据的左上角,其他元素填充0,得到bb,
对bb进行FFT2D正向变换得到B,然后也将图像数据放置到图像的左上角,其他元素填充为0,得到aa,对aa也进行FFT2D正向变换得到B,接着对A和B进行点乘得到C,最后对C进行逆向的FFT变换得到D,最后取D的中间部分有效数据就是卷积的结果。

parallel.gpu.GPUArray

  Python语言由于本身的体质问题,Cython下不能调用多核,只能用多线程。理论上是这样,但还是有很多扩展包能够充分压榨出多核CPU性能,例如multiprocessing是其中的佼佼者。multiprocessing用起来也非常简单,考虑到CPU的多核运算下,每颗核心的算力还是很可观的,所有不能把每个计算组都拆成并行线程,那样内存的读写开销反而会使CPU一直在等待状态,不能一直满负载工作。鉴于此,我设计9879895组线程,每组代表某个特定行与剩下的各个数据行形成的数据组。这样每组线程下的运算量还是比较大的,能使CPU尽可能全在满负载状态。

  举个例子,假设图像数据为:

我们能够使用启用GPU的函数继续对B进行操作。例如,为可视化操作结果,plot命令自动处理GPUArrays。

  代码如下:

       图片 1