| 1454 | } |
| 1455 | |
| 1456 | NOINLINE static void radb4(size_t ido, size_t l1, const double * restrict cc, |
| 1457 | double * restrict ch, const double * restrict wa) |
| 1458 | { |
| 1459 | const size_t cdim=4; |
| 1460 | static const double sqrt2=1.41421356237309504880; |
| 1461 | |
| 1462 | for (size_t k=0; k<l1; k++) |
| 1463 | { |
| 1464 | double tr1, tr2; |
| 1465 | PM (tr2,tr1,CC(0,0,k),CC(ido-1,3,k)) |
| 1466 | double tr3=2.*CC(ido-1,1,k); |
| 1467 | double tr4=2.*CC(0,2,k); |
| 1468 | PM (CH(0,k,0),CH(0,k,2),tr2,tr3) |
| 1469 | PM (CH(0,k,3),CH(0,k,1),tr1,tr4) |
| 1470 | } |
| 1471 | if ((ido&1)==0) |
| 1472 | for (size_t k=0; k<l1; k++) |
| 1473 | { |
| 1474 | double tr1,tr2,ti1,ti2; |
| 1475 | PM (ti1,ti2,CC(0 ,3,k),CC(0 ,1,k)) |
| 1476 | PM (tr2,tr1,CC(ido-1,0,k),CC(ido-1,2,k)) |
| 1477 | CH(ido-1,k,0)=tr2+tr2; |
| 1478 | CH(ido-1,k,1)=sqrt2*(tr1-ti1); |
| 1479 | CH(ido-1,k,2)=ti2+ti2; |
| 1480 | CH(ido-1,k,3)=-sqrt2*(tr1+ti1); |
| 1481 | } |
| 1482 | if (ido<=2) return; |
| 1483 | for (size_t k=0; k<l1;++k) |
| 1484 | for (size_t i=2; i<ido; i+=2) |
| 1485 | { |
| 1486 | double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; |
| 1487 | size_t ic=ido-i; |
| 1488 | PM (tr2,tr1,CC(i-1,0,k),CC(ic-1,3,k)) |
| 1489 | PM (ti1,ti2,CC(i ,0,k),CC(ic ,3,k)) |
| 1490 | PM (tr4,ti3,CC(i ,2,k),CC(ic ,1,k)) |
| 1491 | PM (tr3,ti4,CC(i-1,2,k),CC(ic-1,1,k)) |
| 1492 | PM (CH(i-1,k,0),cr3,tr2,tr3) |
| 1493 | PM (CH(i ,k,0),ci3,ti2,ti3) |
| 1494 | PM (cr4,cr2,tr1,tr4) |
| 1495 | PM (ci2,ci4,ti1,ti4) |
| 1496 | MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ci2,cr2) |
| 1497 | MULPM (CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),ci3,cr3) |
| 1498 | MULPM (CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),ci4,cr4) |
| 1499 | } |
| 1500 | } |
| 1501 | |
| 1502 | NOINLINE static void radb5(size_t ido, size_t l1, const double * restrict cc, |
| 1503 | double * restrict ch, const double * restrict wa) |