| 1560 | #define CH2(a,b) ch[(a)+idl1*(b)] |
| 1561 | |
| 1562 | NOINLINE static void radbg(size_t ido, size_t ip, size_t l1, |
| 1563 | double * restrict cc, double * restrict ch, const double * restrict wa, |
| 1564 | const double * restrict csarr) |
| 1565 | { |
| 1566 | const size_t cdim=ip; |
| 1567 | size_t ipph=(ip+1)/ 2; |
| 1568 | size_t idl1 = ido*l1; |
| 1569 | |
| 1570 | for (size_t k=0; k<l1; ++k) // 102 |
| 1571 | for (size_t i=0; i<ido; ++i) // 101 |
| 1572 | CH(i,k,0) = CC(i,0,k); |
| 1573 | for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc) // 108 |
| 1574 | { |
| 1575 | size_t j2=2*j-1; |
| 1576 | for (size_t k=0; k<l1; ++k) |
| 1577 | { |
| 1578 | CH(0,k,j ) = 2*CC(ido-1,j2,k); |
| 1579 | CH(0,k,jc) = 2*CC(0,j2+1,k); |
| 1580 | } |
| 1581 | } |
| 1582 | |
| 1583 | if (ido!=1) |
| 1584 | { |
| 1585 | for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) // 111 |
| 1586 | { |
| 1587 | size_t j2=2*j-1; |
| 1588 | for (size_t k=0; k<l1; ++k) |
| 1589 | for (size_t i=1, ic=ido-i-2; i<=ido-2; i+=2, ic-=2) // 109 |
| 1590 | { |
| 1591 | CH(i ,k,j ) = CC(i ,j2+1,k)+CC(ic ,j2,k); |
| 1592 | CH(i ,k,jc) = CC(i ,j2+1,k)-CC(ic ,j2,k); |
| 1593 | CH(i+1,k,j ) = CC(i+1,j2+1,k)-CC(ic+1,j2,k); |
| 1594 | CH(i+1,k,jc) = CC(i+1,j2+1,k)+CC(ic+1,j2,k); |
| 1595 | } |
| 1596 | } |
| 1597 | } |
| 1598 | for (size_t l=1,lc=ip-1; l<ipph; ++l,--lc) |
| 1599 | { |
| 1600 | for (size_t ik=0; ik<idl1; ++ik) |
| 1601 | { |
| 1602 | C2(ik,l ) = CH2(ik,0)+csarr[2*l]*CH2(ik,1)+csarr[4*l]*CH2(ik,2); |
| 1603 | C2(ik,lc) = csarr[2*l+1]*CH2(ik,ip-1)+csarr[4*l+1]*CH2(ik,ip-2); |
| 1604 | } |
| 1605 | size_t iang=2*l; |
| 1606 | size_t j=3,jc=ip-3; |
| 1607 | for(; j<ipph-3; j+=4,jc-=4) |
| 1608 | { |
| 1609 | iang+=l; if(iang>ip) iang-=ip; |
| 1610 | double ar1=csarr[2*iang], ai1=csarr[2*iang+1]; |
| 1611 | iang+=l; if(iang>ip) iang-=ip; |
| 1612 | double ar2=csarr[2*iang], ai2=csarr[2*iang+1]; |
| 1613 | iang+=l; if(iang>ip) iang-=ip; |
| 1614 | double ar3=csarr[2*iang], ai3=csarr[2*iang+1]; |
| 1615 | iang+=l; if(iang>ip) iang-=ip; |
| 1616 | double ar4=csarr[2*iang], ai4=csarr[2*iang+1]; |
| 1617 | for (size_t ik=0; ik<idl1; ++ik) |
| 1618 | { |
| 1619 | C2(ik,l ) += ar1*CH2(ik,j )+ar2*CH2(ik,j +1) |