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

只看该作者 · 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);
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)= _packh2(ejbeta_re, ejbeta_im);
abs_a = sshvr(_sadd(_mpyhir(oneOverAbs_a, a_sqr), 1<<(15 - normal)),16-normal) ;
}

复制代码
警告编译器的优化得到的循环迭代信息如下

;* 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++)
{
a_sqr = a.real * a.real + a.imag * a.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= a_sqr * oneOverAbs_a;
ejalpha.real =a.real * oneOverAbs_a;
ejalpha.imag =a.imag * oneOverAbs_a;
}

复制代码
   该循环只需要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]
C0 = A[0] * B[1]

复制代码
为得到复数的实部和虚部，定义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)
QMPYSP: C = A * B for i=0 to 3.

复制代码
使用SIMD修改的程序如下

_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++)
{
dtemp =_amemd8(&a);
/* 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 =a_sqr * oneOverAbs_a;
dtemp1 =_ftod(oneOverAbs_a, oneOverAbs_a);
/* using SIMD DMPYSP for the following operations */
/* ejalpha.real = a.real * oneOverAbs_a;*/
/* ejalpha.imag = a.imag * oneOverAbs_a;*/
_amemd8(&ejalpha) = _dmpysp(dtemp, dtemp1);
}

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

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

相关下载

相关帖子