1* fortran version of *dgpfa5* - 2* radix-5 section of self-sorting, in-place, 3* generalized pfa 4* 5*------------------------------------------------------------------- 6* 7 subroutine dgpfa5f(a,b,trigs,inc,jump,n,mm,lot,isign) 8 double precision a(*), b(*), trigs(*) 9 integer inc, jump, n, mm, lot, isign 10 double precision s, ax, bx, c1, c2, c3 11 double precision t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11 12 double precision u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11 13 double precision co1, co2, co3, co4, si1, si2, si3, si4 14 double precision aja, ajb, ajc, ajd, aje, bjb, bje, bjc 15 double precision bjd, bja, ajf, ajk, bjf, bjk, ajg, ajj 16 double precision ajh, aji, ajl, ajq, bjg, bjj, bjh, bji 17 double precision bjl, bjq, ajo, ajm, ajn, ajr, ajw, bjo 18 double precision bjm, bjn, bjr, bjw, ajt, ajs, ajx, ajp 19 double precision bjt, bjs, bjx, bjp, ajv, ajy, aju, bjv 20 double precision bjy, bju 21 data sin36/0.587785252292473d+0/, 22 * sin72/0.951056516295154d+0/, 23 * qrt5/0.559016994374947d+0/ 24 data lvr/128/ 25* 26* *************************************************************** 27* * * 28* * N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. * 29* * RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE * 30* * (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER. * 31* * * 32* *************************************************************** 33* 34 n5 = 5 ** mm 35 inq = n / n5 36 jstepx = (n5-n) * inc 37 ninc = n * inc 38 ink = inc * inq 39 mu = mod(inq,5) 40 if (isign.eq.-1) mu = 5 - mu 41* 42 m = mm 43 mh = (m+1)/2 44 s = dfloat(isign) 45 c1 = qrt5 46 c2 = sin72 47 c3 = sin36 48 if (mu.eq.2.or.mu.eq.3) then 49 c1 = -c1 50 c2 = sin36 51 c3 = sin72 52 endif 53 if (mu.eq.3.or.mu.eq.4) c2 = -c2 54 if (mu.eq.2.or.mu.eq.4) c3 = -c3 55* 56 nblox = 1 + (lot-1)/lvr 57 left = lot 58 s = dfloat(isign) 59 istart = 1 60* 61* loop on blocks of lvr transforms 62* -------------------------------- 63 do 500 nb = 1 , nblox 64* 65 if (left.le.lvr) then 66 nvex = left 67 else if (left.lt.(2*lvr)) then 68 nvex = left/2 69 nvex = nvex + mod(nvex,2) 70 else 71 nvex = lvr 72 endif 73 left = left - nvex 74* 75 la = 1 76* 77* loop on type I radix-5 passes 78* ----------------------------- 79 do 160 ipass = 1 , mh 80 jstep = (n*inc) / (5*la) 81 jstepl = jstep - ninc 82 kk = 0 83* 84* loop on k 85* --------- 86 do 150 k = 0 , jstep-ink , ink 87* 88 if (k.gt.0) then 89 co1 = trigs(kk+1) 90 si1 = s*trigs(kk+2) 91 co2 = trigs(2*kk+1) 92 si2 = s*trigs(2*kk+2) 93 co3 = trigs(3*kk+1) 94 si3 = s*trigs(3*kk+2) 95 co4 = trigs(4*kk+1) 96 si4 = s*trigs(4*kk+2) 97 endif 98* 99* loop along transform 100* -------------------- 101 do 140 jjj = k , (n-1)*inc , 5*jstep 102 ja = istart + jjj 103* 104* "transverse" loop 105* ----------------- 106 do 135 nu = 1 , inq 107 jb = ja + jstepl 108 if (jb.lt.istart) jb = jb + ninc 109 jc = jb + jstepl 110 if (jc.lt.istart) jc = jc + ninc 111 jd = jc + jstepl 112 if (jd.lt.istart) jd = jd + ninc 113 je = jd + jstepl 114 if (je.lt.istart) je = je + ninc 115 j = 0 116* 117* loop across transforms 118* ---------------------- 119 if (k.eq.0) then 120* 121cdir$ ivdep, shortloop 122 do 110 l = 1 , nvex 123 ajb = a(jb+j) 124 aje = a(je+j) 125 t1 = ajb + aje 126 ajc = a(jc+j) 127 ajd = a(jd+j) 128 t2 = ajc + ajd 129 t3 = ajb - aje 130 t4 = ajc - ajd 131 t5 = t1 + t2 132 t6 = c1 * ( t1 - t2 ) 133 aja = a(ja+j) 134 t7 = aja - 0.25 * t5 135 a(ja+j) = aja + t5 136 t8 = t7 + t6 137 t9 = t7 - t6 138 t10 = c3 * t3 - c2 * t4 139 t11 = c2 * t3 + c3 * t4 140 bjb = b(jb+j) 141 bje = b(je+j) 142 u1 = bjb + bje 143 bjc = b(jc+j) 144 bjd = b(jd+j) 145 u2 = bjc + bjd 146 u3 = bjb - bje 147 u4 = bjc - bjd 148 u5 = u1 + u2 149 u6 = c1 * ( u1 - u2 ) 150 bja = b(ja+j) 151 u7 = bja - 0.25 * u5 152 b(ja+j) = bja + u5 153 u8 = u7 + u6 154 u9 = u7 - u6 155 u10 = c3 * u3 - c2 * u4 156 u11 = c2 * u3 + c3 * u4 157 a(jb+j) = t8 - u11 158 b(jb+j) = u8 + t11 159 a(je+j) = t8 + u11 160 b(je+j) = u8 - t11 161 a(jc+j) = t9 - u10 162 b(jc+j) = u9 + t10 163 a(jd+j) = t9 + u10 164 b(jd+j) = u9 - t10 165 j = j + jump 166 110 continue 167* 168 else 169* 170cdir$ ivdep,shortloop 171 do 130 l = 1 , nvex 172 ajb = a(jb+j) 173 aje = a(je+j) 174 t1 = ajb + aje 175 ajc = a(jc+j) 176 ajd = a(jd+j) 177 t2 = ajc + ajd 178 t3 = ajb - aje 179 t4 = ajc - ajd 180 t5 = t1 + t2 181 t6 = c1 * ( t1 - t2 ) 182 aja = a(ja+j) 183 t7 = aja - 0.25 * t5 184 a(ja+j) = aja + t5 185 t8 = t7 + t6 186 t9 = t7 - t6 187 t10 = c3 * t3 - c2 * t4 188 t11 = c2 * t3 + c3 * t4 189 bjb = b(jb+j) 190 bje = b(je+j) 191 u1 = bjb + bje 192 bjc = b(jc+j) 193 bjd = b(jd+j) 194 u2 = bjc + bjd 195 u3 = bjb - bje 196 u4 = bjc - bjd 197 u5 = u1 + u2 198 u6 = c1 * ( u1 - u2 ) 199 bja = b(ja+j) 200 u7 = bja - 0.25 * u5 201 b(ja+j) = bja + u5 202 u8 = u7 + u6 203 u9 = u7 - u6 204 u10 = c3 * u3 - c2 * u4 205 u11 = c2 * u3 + c3 * u4 206 a(jb+j) = co1*(t8-u11) - si1*(u8+t11) 207 b(jb+j) = si1*(t8-u11) + co1*(u8+t11) 208 a(je+j) = co4*(t8+u11) - si4*(u8-t11) 209 b(je+j) = si4*(t8+u11) + co4*(u8-t11) 210 a(jc+j) = co2*(t9-u10) - si2*(u9+t10) 211 b(jc+j) = si2*(t9-u10) + co2*(u9+t10) 212 a(jd+j) = co3*(t9+u10) - si3*(u9-t10) 213 b(jd+j) = si3*(t9+u10) + co3*(u9-t10) 214 j = j + jump 215 130 continue 216* 217 endif 218* 219*-----( end of loop across transforms ) 220* 221 ja = ja + jstepx 222 if (ja.lt.istart) ja = ja + ninc 223 135 continue 224 140 continue 225*-----( end of loop along transforms ) 226 kk = kk + 2*la 227 150 continue 228*-----( end of loop on nonzero k ) 229 la = 5*la 230 160 continue 231*-----( end of loop on type I radix-5 passes) 232* 233 if (n.eq.5) go to 490 234* 235* loop on type II radix-5 passes 236* ------------------------------ 237 400 continue 238* 239 do 480 ipass = mh+1 , m 240 jstep = (n*inc) / (5*la) 241 jstepl = jstep - ninc 242 laincl = la * ink - ninc 243 kk = 0 244* 245* loop on k 246* --------- 247 do 470 k = 0 , jstep-ink , ink 248* 249 if (k.gt.0) then 250 co1 = trigs(kk+1) 251 si1 = s*trigs(kk+2) 252 co2 = trigs(2*kk+1) 253 si2 = s*trigs(2*kk+2) 254 co3 = trigs(3*kk+1) 255 si3 = s*trigs(3*kk+2) 256 co4 = trigs(4*kk+1) 257 si4 = s*trigs(4*kk+2) 258 endif 259* 260* double loop along first transform in block 261* ------------------------------------------ 262 do 460 ll = k , (la-1)*ink , 5*jstep 263* 264 do 450 jjj = ll , (n-1)*inc , 5*la*ink 265 ja = istart + jjj 266* 267* "transverse" loop 268* ----------------- 269 do 445 nu = 1 , inq 270 jb = ja + jstepl 271 if (jb.lt.istart) jb = jb + ninc 272 jc = jb + jstepl 273 if (jc.lt.istart) jc = jc + ninc 274 jd = jc + jstepl 275 if (jd.lt.istart) jd = jd + ninc 276 je = jd + jstepl 277 if (je.lt.istart) je = je + ninc 278 jf = ja + laincl 279 if (jf.lt.istart) jf = jf + ninc 280 jg = jf + jstepl 281 if (jg.lt.istart) jg = jg + ninc 282 jh = jg + jstepl 283 if (jh.lt.istart) jh = jh + ninc 284 ji = jh + jstepl 285 if (ji.lt.istart) ji = ji + ninc 286 jj = ji + jstepl 287 if (jj.lt.istart) jj = jj + ninc 288 jk = jf + laincl 289 if (jk.lt.istart) jk = jk + ninc 290 jl = jk + jstepl 291 if (jl.lt.istart) jl = jl + ninc 292 jm = jl + jstepl 293 if (jm.lt.istart) jm = jm + ninc 294 jn = jm + jstepl 295 if (jn.lt.istart) jn = jn + ninc 296 jo = jn + jstepl 297 if (jo.lt.istart) jo = jo + ninc 298 jp = jk + laincl 299 if (jp.lt.istart) jp = jp + ninc 300 jq = jp + jstepl 301 if (jq.lt.istart) jq = jq + ninc 302 jr = jq + jstepl 303 if (jr.lt.istart) jr = jr + ninc 304 js = jr + jstepl 305 if (js.lt.istart) js = js + ninc 306 jt = js + jstepl 307 if (jt.lt.istart) jt = jt + ninc 308 ju = jp + laincl 309 if (ju.lt.istart) ju = ju + ninc 310 jv = ju + jstepl 311 if (jv.lt.istart) jv = jv + ninc 312 jw = jv + jstepl 313 if (jw.lt.istart) jw = jw + ninc 314 jx = jw + jstepl 315 if (jx.lt.istart) jx = jx + ninc 316 jy = jx + jstepl 317 if (jy.lt.istart) jy = jy + ninc 318 j = 0 319* 320* loop across transforms 321* ---------------------- 322 if (k.eq.0) then 323* 324cdir$ ivdep, shortloop 325 do 410 l = 1 , nvex 326 ajb = a(jb+j) 327 aje = a(je+j) 328 t1 = ajb + aje 329 ajc = a(jc+j) 330 ajd = a(jd+j) 331 t2 = ajc + ajd 332 t3 = ajb - aje 333 t4 = ajc - ajd 334 ajf = a(jf+j) 335 ajb = ajf 336 t5 = t1 + t2 337 t6 = c1 * ( t1 - t2 ) 338 aja = a(ja+j) 339 t7 = aja - 0.25 * t5 340 a(ja+j) = aja + t5 341 t8 = t7 + t6 342 t9 = t7 - t6 343 ajk = a(jk+j) 344 ajc = ajk 345 t10 = c3 * t3 - c2 * t4 346 t11 = c2 * t3 + c3 * t4 347 bjb = b(jb+j) 348 bje = b(je+j) 349 u1 = bjb + bje 350 bjc = b(jc+j) 351 bjd = b(jd+j) 352 u2 = bjc + bjd 353 u3 = bjb - bje 354 u4 = bjc - bjd 355 bjf = b(jf+j) 356 bjb = bjf 357 u5 = u1 + u2 358 u6 = c1 * ( u1 - u2 ) 359 bja = b(ja+j) 360 u7 = bja - 0.25 * u5 361 b(ja+j) = bja + u5 362 u8 = u7 + u6 363 u9 = u7 - u6 364 bjk = b(jk+j) 365 bjc = bjk 366 u10 = c3 * u3 - c2 * u4 367 u11 = c2 * u3 + c3 * u4 368 a(jf+j) = t8 - u11 369 b(jf+j) = u8 + t11 370 aje = t8 + u11 371 bje = u8 - t11 372 a(jk+j) = t9 - u10 373 b(jk+j) = u9 + t10 374 ajd = t9 + u10 375 bjd = u9 - t10 376*---------------------- 377 ajg = a(jg+j) 378 ajj = a(jj+j) 379 t1 = ajg + ajj 380 ajh = a(jh+j) 381 aji = a(ji+j) 382 t2 = ajh + aji 383 t3 = ajg - ajj 384 t4 = ajh - aji 385 ajl = a(jl+j) 386 ajh = ajl 387 t5 = t1 + t2 388 t6 = c1 * ( t1 - t2 ) 389 t7 = ajb - 0.25 * t5 390 a(jb+j) = ajb + t5 391 t8 = t7 + t6 392 t9 = t7 - t6 393 ajq = a(jq+j) 394 aji = ajq 395 t10 = c3 * t3 - c2 * t4 396 t11 = c2 * t3 + c3 * t4 397 bjg = b(jg+j) 398 bjj = b(jj+j) 399 u1 = bjg + bjj 400 bjh = b(jh+j) 401 bji = b(ji+j) 402 u2 = bjh + bji 403 u3 = bjg - bjj 404 u4 = bjh - bji 405 bjl = b(jl+j) 406 bjh = bjl 407 u5 = u1 + u2 408 u6 = c1 * ( u1 - u2 ) 409 u7 = bjb - 0.25 * u5 410 b(jb+j) = bjb + u5 411 u8 = u7 + u6 412 u9 = u7 - u6 413 bjq = b(jq+j) 414 bji = bjq 415 u10 = c3 * u3 - c2 * u4 416 u11 = c2 * u3 + c3 * u4 417 a(jg+j) = t8 - u11 418 b(jg+j) = u8 + t11 419 ajj = t8 + u11 420 bjj = u8 - t11 421 a(jl+j) = t9 - u10 422 b(jl+j) = u9 + t10 423 a(jq+j) = t9 + u10 424 b(jq+j) = u9 - t10 425*---------------------- 426 ajo = a(jo+j) 427 t1 = ajh + ajo 428 ajm = a(jm+j) 429 ajn = a(jn+j) 430 t2 = ajm + ajn 431 t3 = ajh - ajo 432 t4 = ajm - ajn 433 ajr = a(jr+j) 434 ajn = ajr 435 t5 = t1 + t2 436 t6 = c1 * ( t1 - t2 ) 437 t7 = ajc - 0.25 * t5 438 a(jc+j) = ajc + t5 439 t8 = t7 + t6 440 t9 = t7 - t6 441 ajw = a(jw+j) 442 ajo = ajw 443 t10 = c3 * t3 - c2 * t4 444 t11 = c2 * t3 + c3 * t4 445 bjo = b(jo+j) 446 u1 = bjh + bjo 447 bjm = b(jm+j) 448 bjn = b(jn+j) 449 u2 = bjm + bjn 450 u3 = bjh - bjo 451 u4 = bjm - bjn 452 bjr = b(jr+j) 453 bjn = bjr 454 u5 = u1 + u2 455 u6 = c1 * ( u1 - u2 ) 456 u7 = bjc - 0.25 * u5 457 b(jc+j) = bjc + u5 458 u8 = u7 + u6 459 u9 = u7 - u6 460 bjw = b(jw+j) 461 bjo = bjw 462 u10 = c3 * u3 - c2 * u4 463 u11 = c2 * u3 + c3 * u4 464 a(jh+j) = t8 - u11 465 b(jh+j) = u8 + t11 466 a(jw+j) = t8 + u11 467 b(jw+j) = u8 - t11 468 a(jm+j) = t9 - u10 469 b(jm+j) = u9 + t10 470 a(jr+j) = t9 + u10 471 b(jr+j) = u9 - t10 472*---------------------- 473 ajt = a(jt+j) 474 t1 = aji + ajt 475 ajs = a(js+j) 476 t2 = ajn + ajs 477 t3 = aji - ajt 478 t4 = ajn - ajs 479 ajx = a(jx+j) 480 ajt = ajx 481 t5 = t1 + t2 482 t6 = c1 * ( t1 - t2 ) 483 ajp = a(jp+j) 484 t7 = ajp - 0.25 * t5 485 ax = ajp + t5 486 t8 = t7 + t6 487 t9 = t7 - t6 488 a(jp+j) = ajd 489 t10 = c3 * t3 - c2 * t4 490 t11 = c2 * t3 + c3 * t4 491 a(jd+j) = ax 492 bjt = b(jt+j) 493 u1 = bji + bjt 494 bjs = b(js+j) 495 u2 = bjn + bjs 496 u3 = bji - bjt 497 u4 = bjn - bjs 498 bjx = b(jx+j) 499 bjt = bjx 500 u5 = u1 + u2 501 u6 = c1 * ( u1 - u2 ) 502 bjp = b(jp+j) 503 u7 = bjp - 0.25 * u5 504 bx = bjp + u5 505 u8 = u7 + u6 506 u9 = u7 - u6 507 b(jp+j) = bjd 508 u10 = c3 * u3 - c2 * u4 509 u11 = c2 * u3 + c3 * u4 510 b(jd+j) = bx 511 a(ji+j) = t8 - u11 512 b(ji+j) = u8 + t11 513 a(jx+j) = t8 + u11 514 b(jx+j) = u8 - t11 515 a(jn+j) = t9 - u10 516 b(jn+j) = u9 + t10 517 a(js+j) = t9 + u10 518 b(js+j) = u9 - t10 519*---------------------- 520 ajv = a(jv+j) 521 ajy = a(jy+j) 522 t1 = ajv + ajy 523 t2 = ajo + ajt 524 t3 = ajv - ajy 525 t4 = ajo - ajt 526 a(jv+j) = ajj 527 t5 = t1 + t2 528 t6 = c1 * ( t1 - t2 ) 529 aju = a(ju+j) 530 t7 = aju - 0.25 * t5 531 ax = aju + t5 532 t8 = t7 + t6 533 t9 = t7 - t6 534 a(ju+j) = aje 535 t10 = c3 * t3 - c2 * t4 536 t11 = c2 * t3 + c3 * t4 537 a(je+j) = ax 538 bjv = b(jv+j) 539 bjy = b(jy+j) 540 u1 = bjv + bjy 541 u2 = bjo + bjt 542 u3 = bjv - bjy 543 u4 = bjo - bjt 544 b(jv+j) = bjj 545 u5 = u1 + u2 546 u6 = c1 * ( u1 - u2 ) 547 bju = b(ju+j) 548 u7 = bju - 0.25 * u5 549 bx = bju + u5 550 u8 = u7 + u6 551 u9 = u7 - u6 552 b(ju+j) = bje 553 u10 = c3 * u3 - c2 * u4 554 u11 = c2 * u3 + c3 * u4 555 b(je+j) = bx 556 a(jj+j) = t8 - u11 557 b(jj+j) = u8 + t11 558 a(jy+j) = t8 + u11 559 b(jy+j) = u8 - t11 560 a(jo+j) = t9 - u10 561 b(jo+j) = u9 + t10 562 a(jt+j) = t9 + u10 563 b(jt+j) = u9 - t10 564 j = j + jump 565 410 continue 566* 567 else 568* 569cdir$ ivdep, shortloop 570 do 440 l = 1 , nvex 571 ajb = a(jb+j) 572 aje = a(je+j) 573 t1 = ajb + aje 574 ajc = a(jc+j) 575 ajd = a(jd+j) 576 t2 = ajc + ajd 577 t3 = ajb - aje 578 t4 = ajc - ajd 579 ajf = a(jf+j) 580 ajb = ajf 581 t5 = t1 + t2 582 t6 = c1 * ( t1 - t2 ) 583 aja = a(ja+j) 584 t7 = aja - 0.25 * t5 585 a(ja+j) = aja + t5 586 t8 = t7 + t6 587 t9 = t7 - t6 588 ajk = a(jk+j) 589 ajc = ajk 590 t10 = c3 * t3 - c2 * t4 591 t11 = c2 * t3 + c3 * t4 592 bjb = b(jb+j) 593 bje = b(je+j) 594 u1 = bjb + bje 595 bjc = b(jc+j) 596 bjd = b(jd+j) 597 u2 = bjc + bjd 598 u3 = bjb - bje 599 u4 = bjc - bjd 600 bjf = b(jf+j) 601 bjb = bjf 602 u5 = u1 + u2 603 u6 = c1 * ( u1 - u2 ) 604 bja = b(ja+j) 605 u7 = bja - 0.25 * u5 606 b(ja+j) = bja + u5 607 u8 = u7 + u6 608 u9 = u7 - u6 609 bjk = b(jk+j) 610 bjc = bjk 611 u10 = c3 * u3 - c2 * u4 612 u11 = c2 * u3 + c3 * u4 613 a(jf+j) = co1*(t8-u11) - si1*(u8+t11) 614 b(jf+j) = si1*(t8-u11) + co1*(u8+t11) 615 aje = co4*(t8+u11) - si4*(u8-t11) 616 bje = si4*(t8+u11) + co4*(u8-t11) 617 a(jk+j) = co2*(t9-u10) - si2*(u9+t10) 618 b(jk+j) = si2*(t9-u10) + co2*(u9+t10) 619 ajd = co3*(t9+u10) - si3*(u9-t10) 620 bjd = si3*(t9+u10) + co3*(u9-t10) 621*---------------------- 622 ajg = a(jg+j) 623 ajj = a(jj+j) 624 t1 = ajg + ajj 625 ajh = a(jh+j) 626 aji = a(ji+j) 627 t2 = ajh + aji 628 t3 = ajg - ajj 629 t4 = ajh - aji 630 ajl = a(jl+j) 631 ajh = ajl 632 t5 = t1 + t2 633 t6 = c1 * ( t1 - t2 ) 634 t7 = ajb - 0.25 * t5 635 a(jb+j) = ajb + t5 636 t8 = t7 + t6 637 t9 = t7 - t6 638 ajq = a(jq+j) 639 aji = ajq 640 t10 = c3 * t3 - c2 * t4 641 t11 = c2 * t3 + c3 * t4 642 bjg = b(jg+j) 643 bjj = b(jj+j) 644 u1 = bjg + bjj 645 bjh = b(jh+j) 646 bji = b(ji+j) 647 u2 = bjh + bji 648 u3 = bjg - bjj 649 u4 = bjh - bji 650 bjl = b(jl+j) 651 bjh = bjl 652 u5 = u1 + u2 653 u6 = c1 * ( u1 - u2 ) 654 u7 = bjb - 0.25 * u5 655 b(jb+j) = bjb + u5 656 u8 = u7 + u6 657 u9 = u7 - u6 658 bjq = b(jq+j) 659 bji = bjq 660 u10 = c3 * u3 - c2 * u4 661 u11 = c2 * u3 + c3 * u4 662 a(jg+j) = co1*(t8-u11) - si1*(u8+t11) 663 b(jg+j) = si1*(t8-u11) + co1*(u8+t11) 664 ajj = co4*(t8+u11) - si4*(u8-t11) 665 bjj = si4*(t8+u11) + co4*(u8-t11) 666 a(jl+j) = co2*(t9-u10) - si2*(u9+t10) 667 b(jl+j) = si2*(t9-u10) + co2*(u9+t10) 668 a(jq+j) = co3*(t9+u10) - si3*(u9-t10) 669 b(jq+j) = si3*(t9+u10) + co3*(u9-t10) 670*---------------------- 671 ajo = a(jo+j) 672 t1 = ajh + ajo 673 ajm = a(jm+j) 674 ajn = a(jn+j) 675 t2 = ajm + ajn 676 t3 = ajh - ajo 677 t4 = ajm - ajn 678 ajr = a(jr+j) 679 ajn = ajr 680 t5 = t1 + t2 681 t6 = c1 * ( t1 - t2 ) 682 t7 = ajc - 0.25 * t5 683 a(jc+j) = ajc + t5 684 t8 = t7 + t6 685 t9 = t7 - t6 686 ajw = a(jw+j) 687 ajo = ajw 688 t10 = c3 * t3 - c2 * t4 689 t11 = c2 * t3 + c3 * t4 690 bjo = b(jo+j) 691 u1 = bjh + bjo 692 bjm = b(jm+j) 693 bjn = b(jn+j) 694 u2 = bjm + bjn 695 u3 = bjh - bjo 696 u4 = bjm - bjn 697 bjr = b(jr+j) 698 bjn = bjr 699 u5 = u1 + u2 700 u6 = c1 * ( u1 - u2 ) 701 u7 = bjc - 0.25 * u5 702 b(jc+j) = bjc + u5 703 u8 = u7 + u6 704 u9 = u7 - u6 705 bjw = b(jw+j) 706 bjo = bjw 707 u10 = c3 * u3 - c2 * u4 708 u11 = c2 * u3 + c3 * u4 709 a(jh+j) = co1*(t8-u11) - si1*(u8+t11) 710 b(jh+j) = si1*(t8-u11) + co1*(u8+t11) 711 a(jw+j) = co4*(t8+u11) - si4*(u8-t11) 712 b(jw+j) = si4*(t8+u11) + co4*(u8-t11) 713 a(jm+j) = co2*(t9-u10) - si2*(u9+t10) 714 b(jm+j) = si2*(t9-u10) + co2*(u9+t10) 715 a(jr+j) = co3*(t9+u10) - si3*(u9-t10) 716 b(jr+j) = si3*(t9+u10) + co3*(u9-t10) 717*---------------------- 718 ajt = a(jt+j) 719 t1 = aji + ajt 720 ajs = a(js+j) 721 t2 = ajn + ajs 722 t3 = aji - ajt 723 t4 = ajn - ajs 724 ajx = a(jx+j) 725 ajt = ajx 726 t5 = t1 + t2 727 t6 = c1 * ( t1 - t2 ) 728 ajp = a(jp+j) 729 t7 = ajp - 0.25 * t5 730 ax = ajp + t5 731 t8 = t7 + t6 732 t9 = t7 - t6 733 a(jp+j) = ajd 734 t10 = c3 * t3 - c2 * t4 735 t11 = c2 * t3 + c3 * t4 736 a(jd+j) = ax 737 bjt = b(jt+j) 738 u1 = bji + bjt 739 bjs = b(js+j) 740 u2 = bjn + bjs 741 u3 = bji - bjt 742 u4 = bjn - bjs 743 bjx = b(jx+j) 744 bjt = bjx 745 u5 = u1 + u2 746 u6 = c1 * ( u1 - u2 ) 747 bjp = b(jp+j) 748 u7 = bjp - 0.25 * u5 749 bx = bjp + u5 750 u8 = u7 + u6 751 u9 = u7 - u6 752 b(jp+j) = bjd 753 u10 = c3 * u3 - c2 * u4 754 u11 = c2 * u3 + c3 * u4 755 b(jd+j) = bx 756 a(ji+j) = co1*(t8-u11) - si1*(u8+t11) 757 b(ji+j) = si1*(t8-u11) + co1*(u8+t11) 758 a(jx+j) = co4*(t8+u11) - si4*(u8-t11) 759 b(jx+j) = si4*(t8+u11) + co4*(u8-t11) 760 a(jn+j) = co2*(t9-u10) - si2*(u9+t10) 761 b(jn+j) = si2*(t9-u10) + co2*(u9+t10) 762 a(js+j) = co3*(t9+u10) - si3*(u9-t10) 763 b(js+j) = si3*(t9+u10) + co3*(u9-t10) 764*---------------------- 765 ajv = a(jv+j) 766 ajy = a(jy+j) 767 t1 = ajv + ajy 768 t2 = ajo + ajt 769 t3 = ajv - ajy 770 t4 = ajo - ajt 771 a(jv+j) = ajj 772 t5 = t1 + t2 773 t6 = c1 * ( t1 - t2 ) 774 aju = a(ju+j) 775 t7 = aju - 0.25 * t5 776 ax = aju + t5 777 t8 = t7 + t6 778 t9 = t7 - t6 779 a(ju+j) = aje 780 t10 = c3 * t3 - c2 * t4 781 t11 = c2 * t3 + c3 * t4 782 a(je+j) = ax 783 bjv = b(jv+j) 784 bjy = b(jy+j) 785 u1 = bjv + bjy 786 u2 = bjo + bjt 787 u3 = bjv - bjy 788 u4 = bjo - bjt 789 b(jv+j) = bjj 790 u5 = u1 + u2 791 u6 = c1 * ( u1 - u2 ) 792 bju = b(ju+j) 793 u7 = bju - 0.25 * u5 794 bx = bju + u5 795 u8 = u7 + u6 796 u9 = u7 - u6 797 b(ju+j) = bje 798 u10 = c3 * u3 - c2 * u4 799 u11 = c2 * u3 + c3 * u4 800 b(je+j) = bx 801 a(jj+j) = co1*(t8-u11) - si1*(u8+t11) 802 b(jj+j) = si1*(t8-u11) + co1*(u8+t11) 803 a(jy+j) = co4*(t8+u11) - si4*(u8-t11) 804 b(jy+j) = si4*(t8+u11) + co4*(u8-t11) 805 a(jo+j) = co2*(t9-u10) - si2*(u9+t10) 806 b(jo+j) = si2*(t9-u10) + co2*(u9+t10) 807 a(jt+j) = co3*(t9+u10) - si3*(u9-t10) 808 b(jt+j) = si3*(t9+u10) + co3*(u9-t10) 809 j = j + jump 810 440 continue 811* 812 endif 813* 814*-----(end of loop across transforms) 815* 816 ja = ja + jstepx 817 if (ja.lt.istart) ja = ja + ninc 818 445 continue 819 450 continue 820 460 continue 821*-----( end of double loop for this k ) 822 kk = kk + 2*la 823 470 continue 824*-----( end of loop over values of k ) 825 la = 5*la 826 480 continue 827*-----( end of loop on type II radix-5 passes ) 828*-----( nvex transforms completed) 829 490 continue 830 istart = istart + nvex * jump 831 500 continue 832*-----( end of loop on blocks of transforms ) 833* 834 return 835 end 836