打印

C66x定点浮点混合DSP循环编程优化指南(2)

[复制链接]
1343|0
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
Jasmines|  楼主 | 2017-11-10 10:33 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
C66x定点浮点混合DSP循环编程优化指南


考虑C64x+平台的优化
        优化点包括:
        快速计算1/sqrt(x)的函数,需要inline的函数:可以考虑使用分段查找表来计算该函数,或者采用查表和其他的迭代算法结合的方法;
        考虑适当的scale缩放和Q值调整;
        考虑更高位宽的数据读写;
        考虑定点的复数乘法SIMD指令
  • _nassert((int) a % 8 == 0);
  • _nassert((int) ejalpha % 8 == 0);
  • _nassert((int) abs_a % 8 == 0);
  • #pragma MUST_ITERATE(4,100, 4);
  • #pragma UNROLL(2);
  • for ( i = 0; i < n; i++)
  • {
  • temp1 = _amem4(&a<span style="font-style: italic;"><span style="font-style: normal;">);
  • a_sqr = _dotp2(temp1, temp1);
  • /* 1/sqrt(a_sqr) */
  • normal = _norm(a_sqr);
  • normal = normal & 0xFFFFFFFE;
  • x_norm = _sshvl(a_sqr, normal);
  • normal = normal >> 1;
  • Index = _sshvr(_sadd(x_norm,0x800000),24);
  • oneOverAbs_a=_mpylir( xcbia[Index], x_norm );
  • oneOverAbs_a=_sadd((int)x3sa[Index]<<16,
  • _sshvr(oneOverAbs_a,ShiftValDifp1a[Index]));
  • normal =15 - ShiftVala[Index] + normal;
  • ejbeta_re =_sadd(_sshvl(_mpyhir(temp1, oneOverAbs_a), normal - 1), 0x8000);
  • ejbeta_im =_sadd(_sshvl(_mpylir(temp1, oneOverAbs_a), normal - 1), 0x8000);
  • _amem4(&ejalpha</span><span style="font-style: normal;">)= _packh2(ejbeta_re, ejbeta_im);
  • abs_a</span><span style="font-style: normal;"> = sshvr(_sadd(_mpyhir(oneOverAbs_a, a_sqr), 1<<(15 - normal)),16-normal) ;
  • }</span></span>

复制代码
       警告编译器的优化得到的循环迭代信息如下
  • ;* SOFTWARE PIPELINE INFORMATION
  • ;*
  • ;* Loop source line : 198
  • ;* Loop opening brace source line : 199
  • ;* Loop closing brace source line : 218
  • ;* Loop Unroll Multiple : 2x
  • ;* Known Minimum Trip Count : 2
  • ;* Known Maximum Trip Count : 50
  • ;* Known Max Trip Count Factor : 2
  • ;* Loop Carried Dependency Bound(^) : 0
  • ;* Unpartitioned Resource Bound : 9
  • ;* Partitioned Resource Bound(*) : 9
  • ;* Resource Partition:
  • ;* A-side B-side
  • ;* .L units 1 1
  • ;* .S units 6 6
  • ;* .D units 7 6
  • ;* .M units 9* 9*
  • ;* .X cross paths 5 1
  • ;* .T address paths 7 6
  • ;* Long read paths 0 0
  • ;* Long write paths 0 0
  • ;* Logical ops (.LS) 8 7 (.L or .S unit)
  • ;* Addition ops (.LSD) 5 7 (.L or .S or .D unit)
  • ;* Bound(.L .S .LS) 8 7
  • ;* Bound(.L .S .D .LS .LSD) 9* 9*
  • ;*
  • ;* Searching for software pipeline schedule at ...
  • ;* ii = 9 Register is live too long
  • ;* ii = 9 Did not find schedule
  • ;* ii = 10 Register is live too long
  • ;* ii = 10 Register is live too long
  • ;* ii = 10 Did not find schedule
  • ;* ii = 11 Schedule found with 5 iterations in parallel
  • ;* Done

复制代码
      
         以上处理中1/sqrt(x)被inline到循环内,_nassert(), restrict关键字以及#pragmas来告诉编译器尽可能多的关于buffer重叠、数据对其、循环次数和unroll信息。采用_amem4()来进行数据读写,使用 MPYLIR和MPYHIR左16-bit乘以32-bit,使用DOTP2计算16-bit复数的能量。其他的还有round和scale的操作来做Q值调整。现在100个数据的求解幅值和相位仅需要746个周期。
       考虑C66x平台的优化
        针对C66x平台,有单独的指令RSQRSP来计算1/sqrt(x),复数的乘法和浮点的运算以及定点和浮点的转换还可以考虑并行的SIMD指令。C66x提供了单周期的指令RSQRSP 来计算 1/sqrt(x),RCPSP来计算1/x。另外就是采用编译器友好的关键字_nassert()和restrict以及#pragmas来告诉编译器尽可能多的关于buffer重叠、数据对其、循环次数和unroll信息,另外就是RSQRSP和RCPSP以及双精度格式的RSQRDP、RCPDP能得到正确的指数exponent值和只有8比特精度的尾数mantissa值,因而要更精确的值,需要一些迭代算法,如Newton-Raphson迭代公式x[n+1]=x[n]*(2 - v *x[n]) 来计算倒数,用x[n+1]= x[n]*(1.5 - (v/2)*x[n]*x[n])迭代来计算1/sqrt(v)。多一次迭代多8比特精度,两次迭代24比特尾数精度,3次迭代32比特精度。同时这种实现也不需要查找表了。

  • _nassert((int) a % 8 == 0);
  • _nassert((int) ejalpha % 8 == 0);
  • _nassert((int) abs_a % 8 == 0);
  • #pragma MUST_ITERATE(4,100, 4);
  • #pragma UNROLL(2);
  • for ( i = 0; i < n; i++)
  • {
  • <span style="font-style: italic;"><span style="font-style: normal;">a_sqr = a.real * a.real + a</span><span style="font-style: italic;"><span style="font-style: normal;">.imag * a</span><span style="font-style: italic;"><span style="font-style: normal;">.imag;
  • oneOverAbs_a = _rsqrsp(a_sqr); /* 1/sqrt() instruction 8-bit mantissa
  • precision*/
  • /* One interpolation*/
  • oneOverAbs_a = oneOverAbs_a * (1.5f - (a_sqr/2.f)* oneOverAbs_a
  • *oneOverAbs_a);
  • abs_a</span><span style="font-style: italic;"><span style="font-style: normal;">= a_sqr * oneOverAbs_a;
  • ejalpha</span><span style="font-style: italic;"><span style="font-style: normal;">.real =a</span><span style="font-style: italic;"><span style="font-style: normal;">.real * oneOverAbs_a;
  • ejalpha</span><span style="font-style: italic;"><span style="font-style: normal;">.imag =a</span><span style="font-style: italic;"><span style="font-style: normal;">.imag * oneOverAbs_a;
  • }</span></span></span></span></span></span></span></span></span>

复制代码
       该循环只需要496个周期,而且不用额外的存储查找表。

        继续C66x平台的优化
        C66x有更强的SIMD处理能力,进一步的优化可以考虑浮点乘法、数据加载等。如数据加载采用_amem8(addr)来加载8字节对齐的整型数据,_amemd8(addr) 来加载8字节对齐的浮点数据。因而定义如下的复数数据结构,就能用_amemd8(addr)一次加载一个复数的实部和虚部了,注意一下定义同时考虑了大小端,保证加载的高4字节是虚步,低4字节是实部。
  • #ifdef _LITTLE_ENDIAN
  • typedef struct _CPLXF
  • {
  • float imag;
  • float real;
  • } cplxf_t;
  • #else
  • typedef struct _CPLXF
  • {
  • float real;
  • float imag;
  • } cplxf_t;
  • #endif

复制代码
       对于加载的寄存器对的数据,可以采用_hif(src)来得到实部,用_lof(src)来得到虚部,类似于C64x+中的_hill(src) 和 _loll(src)。组成一个寄存器对,可以采用类似_itoll(srcq, src2) 的_fod(src1, scr2).
        C66x处理器除了提供C674x+ DSP已经包含的MPYSP (SPxSP->SP), MPYDP (DPxDP->DP), MPYSPDP(SPxDP->DP), 以及MPYSP2DP (SPxSP->DP)外,还有新的DMPYSP以及CMPYSP和QMPYSP指令。

  •        DMPYSP:浮点的C = A * B for i=0 to 1
  •         CMPYSP: 用于浮点数据的复数乘法,结果为128-bit格式
  •         C3 = A[1] * B[1]
  •         C2 = A[1] * B[0]
  •         C1 = -A[0] * B[0]
  • <span style="font-style: italic;"><span style="font-style: normal;">        </span><span style="font-style: italic;"><span style="font-style: normal;">C0 = A[0] * B[1]</span></span></span>

复制代码
       为得到复数的实部和虚部,定义128-bit的数据C:
  • __x128_t C_128;
  • C_128 = _cmpysp(A, B);
  • C = _daddsp(_hid128(C_128), _lod128(C_128)),直接得到C3+C1 和 C2+C0。

复制代码
        或者使用
  • intrinsics _complex_mpysp():
  • C=_complex_mpysp(A, B)

复制代码
       得到A和B的共轭的乘积,考虑如下步骤:
  • __x128_t C_128;
  • C_128 = _cmpysp(B, A);
  • C3 = B[1] * A[1]
  • C2 = B[1] * A[0]
  • C1 = -B[0] * A[0]
  • C0 = B[0] * A[1]
  • C = _dsubsp(_hid128(C_128), _lod128(C_128))

复制代码

        直接得到C3-C1 和 C2-C0;
        或者采用简单的
  • intrinsics _complex_conjugate_mpysp():
  • C=_complex_conjugate_mpysp(A, B)
  • <span style="font-style: italic;"><span style="font-style: normal;">QMPYSP: C = A * B</span><span style="font-style: italic;"><span style="font-style: normal;"> for i=0 to 3.</span></span></span>

复制代码
       使用SIMD修改的程序如下
  • <span style="font-style: normal;">_nassert(n % 4 == 0);
  • _nassert((int) a % 8 == 0);
  • _nassert((int) ejalpha % 8 == 0);
  • _nassert((int) abs_a % 8 == 0);
  • #pragma MUST_ITERATE(4,100, 4);
  • for ( i = 0; i < n; i++)
  • {
  • </span>dtemp =_amemd8(&a<span style="font-style: italic;"><span style="font-style: normal;">);
  • /* using SIMD CMPYSP with conjugation for power calculation */
  • a_sqr =_hif(_complex_conjugate_mpysp(dtemp, dtemp));
  • /* or use the following */
  • /* dtemp2 = _dmpysp(dtemp, dtemp); */
  • /* a_sqr = _hif(dtemp2) + _lof(dtemp2); */
  • oneOverAbs_a = _rsqrsp(a_sqr); /* 1/sqrt() instruction 8-bit mantissa precision*/
  • /* 1st interpolation*/
  • oneOverAbs_a = oneOverAbs_a * (1.5f - (a_sqr/2.f)* oneOverAbs_a*oneOverAbs_a);
  • abs_a</span><span style="font-style: italic;"><span style="font-style: normal;"> =a_sqr * oneOverAbs_a;
  • dtemp1 =_ftod(oneOverAbs_a, oneOverAbs_a);
  • /* using SIMD DMPYSP for the following operations */
  • /* ejalpha</span><span style="font-style: italic;"><span style="font-style: normal;">.real = a</span><span style="font-style: italic;"><span style="font-style: normal;">.real * oneOverAbs_a;*/
  • /* ejalpha</span><span style="font-style: italic;"><span style="font-style: normal;">.imag = a</span><span style="font-style: italic;"><span style="font-style: normal;">.imag * oneOverAbs_a;*/
  • _amemd8(&ejalpha</span><span style="font-style: italic;"><span style="font-style: normal;">) = _dmpysp(dtemp, dtemp1);
  • }</span></span></span></span></span></span></span></span>

复制代码
       _complex_conjugate_mpysp(a, a)计算能量,只保存32 MSB到寄存器,还可以通过DMPYSP(a, a) 和_hif() + _lof()来计算能量。
        得到的结果是需要419个周期。


相关帖子

发新帖 我要提问
您需要登录后才可以回帖 登录 | 注册

本版积分规则

745

主题

1077

帖子

10

粉丝