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个周期。
|
|