1! Fast Fourier/Cosine/Sine Transform 2! dimension :two 3! data length :power of 2 4! decimation :frequency 5! radix :4, 2, row-column 6! data :inplace 7! table :use 8! subroutines 9! cdft2d: Complex Discrete Fourier Transform 10! rdft2d: Real Discrete Fourier Transform 11! ddct2d: Discrete Cosine Transform 12! ddst2d: Discrete Sine Transform 13! 14! 15! -------- Complex DFT (Discrete Fourier Transform) -------- 16! [definition] 17! <case1> 18! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) * 19! exp(2*pi*i*j1*k1/n1) * 20! exp(2*pi*i*j2*k2/n2), 21! 0<=k1<n1, 0<=k2<n2 22! <case2> 23! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) * 24! exp(-2*pi*i*j1*k1/n1) * 25! exp(-2*pi*i*j2*k2/n2), 26! 0<=k1<n1, 0<=k2<n2 27! (notes: sum_j=0^n-1 is a summation from j=0 to n-1) 28! [usage] 29! <case1> 30! ip(0) = 0 ! first time only 31! call cdft2d(n1max, 2*n1, n2, 1, a, ip, w) 32! <case2> 33! ip(0) = 0 ! first time only 34! call cdft2d(n1max, 2*n1, n2, -1, a, ip, w) 35! [parameters] 36! n1max :row size of the 2D array (integer) 37! 2*n1 :data length (integer) 38! n1 >= 1, n1 = power of 2 39! n2 :data length (integer) 40! n2 >= 1, n2 = power of 2 41! a(0:2*n1-1,0:n2-1) 42! :input/output data (real*8) 43! input data 44! a(2*j1,j2) = Re(x(j1,j2)), 45! a(2*j1+1,j2) = Im(x(j1,j2)), 46! 0<=j1<n1, 0<=j2<n2 47! output data 48! a(2*k1,k2) = Re(X(k1,k2)), 49! a(2*k1+1,k2) = Im(X(k1,k2)), 50! 0<=k1<n1, 0<=k2<n2 51! ip(0:*):work area for bit reversal (integer) 52! length of ip >= 2+sqrt(n) 53! (n = max(n1, n2)) 54! ip(0),ip(1) are pointers of the cos/sin table. 55! w(0:*) :cos/sin table (real*8) 56! length of w >= max(n1/2, n2/2) 57! w(),ip() are initialized if ip(0) = 0. 58! [remark] 59! Inverse of 60! call cdft2d(n1max, 2*n1, n2, -1, a, ip, w) 61! is 62! call cdft2d(n1max, 2*n1, n2, 1, a, ip, w) 63! do j2 = 0, n2 - 1 64! do j1 = 0, 2 * n1 - 1 65! a(j1, j2) = a(j1, j2) * (1.0d0 / (n1 * n2)) 66! end do 67! end do 68! . 69! 70! 71! -------- Real DFT / Inverse of Real DFT -------- 72! [definition] 73! <case1> RDFT 74! R(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 75! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), 76! 0<=k1<n1, 0<=k2<n2 77! I(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 78! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), 79! 0<=k1<n1, 0<=k2<n2 80! <case2> IRDFT (excluding scale) 81! a(k1,k2) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 82! (R(j1,j2) * 83! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2) + 84! I(j1,j2) * 85! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2)), 86! 0<=k1<n1, 0<=k2<n2 87! (notes: R(n1-k1,n2-k2) = R(k1,k2), 88! I(n1-k1,n2-k2) = -I(k1,k2), 89! R(n1-k1,0) = R(k1,0), 90! I(n1-k1,0) = -I(k1,0), 91! R(0,n2-k2) = R(0,k2), 92! I(0,n2-k2) = -I(0,k2), 93! 0<k1<n1, 0<k2<n2) 94! [usage] 95! <case1> 96! ip(0) = 0 ! first time only 97! call rdft2d(n1max, n1, n2, 1, a, ip, w) 98! <case2> 99! ip(0) = 0 ! first time only 100! call rdft2d(n1max, n1, n2, -1, a, ip, w) 101! [parameters] 102! n1max :row size of the 2D array (integer) 103! n1 :data length (integer) 104! n1 >= 2, n1 = power of 2 105! n2 :data length (integer) 106! n2 >= 2, n2 = power of 2 107! a(0:n1-1,0:n2-1) 108! :input/output data (real*8) 109! <case1> 110! output data 111! a(2*k1,k2) = R(k1,k2) = R(n1-k1,n2-k2), 112! a(2*k1+1,k2) = I(k1,k2) = -I(n1-k1,n2-k2), 113! 0<k1<n1/2, 0<k2<n2, 114! a(2*k1,0) = R(k1,0) = R(n1-k1,0), 115! a(2*k1+1,0) = I(k1,0) = -I(n1-k1,0), 116! 0<k1<n1/2, 117! a(0,k2) = R(0,k2) = R(0,n2-k2), 118! a(1,k2) = I(0,k2) = -I(0,n2-k2), 119! a(1,n2-k2) = R(n1/2,k2) = R(n1/2,n2-k2), 120! a(0,n2-k2) = -I(n1/2,k2) = I(n1/2,n2-k2), 121! 0<k2<n2/2, 122! a(0,0) = R(0,0), 123! a(1,0) = R(n1/2,0), 124! a(0,n2/2) = R(0,n2/2), 125! a(1,n2/2) = R(n1/2,n2/2) 126! <case2> 127! input data 128! a(2*j1,j2) = R(j1,j2) = R(n1-j1,n2-j2), 129! a(2*j1+1,j2) = I(j1,j2) = -I(n1-j1,n2-j2), 130! 0<j1<n1/2, 0<j2<n2, 131! a(2*j1,0) = R(j1,0) = R(n1-j1,0), 132! a(2*j1+1,0) = I(j1,0) = -I(n1-j1,0), 133! 0<j1<n1/2, 134! a(0,j2) = R(0,j2) = R(0,n2-j2), 135! a(1,j2) = I(0,j2) = -I(0,n2-j2), 136! a(1,n2-j2) = R(n1/2,j2) = R(n1/2,n2-j2), 137! a(0,n2-j2) = -I(n1/2,j2) = I(n1/2,n2-j2), 138! 0<j2<n2/2, 139! a(0,0) = R(0,0), 140! a(1,0) = R(n1/2,0), 141! a(0,n2/2) = R(0,n2/2), 142! a(1,n2/2) = R(n1/2,n2/2) 143! ip(0:*):work area for bit reversal (integer) 144! length of ip >= 2+sqrt(n) 145! (n = max(n1/2, n2)) 146! ip(0),ip(1) are pointers of the cos/sin table. 147! w(0:*) :cos/sin table (real*8) 148! length of w >= max(n1/4, n2/2) + n1/4 149! w(),ip() are initialized if ip(0) = 0. 150! [remark] 151! Inverse of 152! call rdft2d(n1max, n1, n2, 1, a, ip, w) 153! is 154! call rdft2d(n1max, n1, n2, -1, a, ip, w) 155! do j2 = 0, n2 - 1 156! do j1 = 0, n1 - 1 157! a(j1, j2) = a(j1, j2) * (2.0d0 / (n1 * n2)) 158! end do 159! end do 160! . 161! 162! 163! -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- 164! [definition] 165! <case1> IDCT (excluding scale) 166! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 167! cos(pi*j1*(k1+1/2)/n1) * 168! cos(pi*j2*(k2+1/2)/n2), 169! 0<=k1<n1, 0<=k2<n2 170! <case2> DCT 171! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 172! cos(pi*(j1+1/2)*k1/n1) * 173! cos(pi*(j2+1/2)*k2/n2), 174! 0<=k1<n1, 0<=k2<n2 175! [usage] 176! <case1> 177! ip(0) = 0 ! first time only 178! call ddct2d(n1max, n1, n2, 1, a, t, ip, w) 179! <case2> 180! ip(0) = 0 ! first time only 181! call ddct2d(n1max, n1, n2, -1, a, t, ip, w) 182! [parameters] 183! n1max :row size of the 2D array (integer) 184! n1 :data length (integer) 185! n1 >= 2, n1 = power of 2 186! n2 :data length (integer) 187! n2 >= 2, n2 = power of 2 188! a(0:n1-1,0:n2-1) 189! :input/output data (real*8) 190! output data 191! a(k1,k2) = C(k1,k2), 0<=k1<n1, 0<=k2<n2 192! t(0:n1-1,0:n2-1) 193! :work area (real*8) 194! ip(0:*):work area for bit reversal (integer) 195! length of ip >= 2+sqrt(n) 196! (n = max(n1/2, n2)) 197! ip(0),ip(1) are pointers of the cos/sin table. 198! w(0:*) :cos/sin table (real*8) 199! length of w >= max(n1/4, n2/2) + max(n1, n2) 200! w(),ip() are initialized if ip(0) = 0. 201! [remark] 202! Inverse of 203! call ddct2d(n1max, n1, n2, -1, a, t, ip, w) 204! is 205! do j1 = 0, n1 - 1 206! a(j1, 0) = a(j1, 0) * 0.5d0 207! end do 208! do j2 = 0, n2 - 1 209! a(0, j2) = a(0, j2) * 0.5d0 210! end do 211! call ddct2d(n1max, n1, n2, 1, a, t, ip, w) 212! do j2 = 0, n2 - 1 213! do j1 = 0, n1 - 1 214! a(j1, j2) = a(j1, j2) * (4.0d0 / (n1 * n2)) 215! end do 216! end do 217! . 218! 219! 220! -------- DST (Discrete Sine Transform) / Inverse of DST -------- 221! [definition] 222! <case1> IDST (excluding scale) 223! S(k1,k2) = sum_j1=1^n1 sum_j2=1^n2 A(j1,j2) * 224! sin(pi*j1*(k1+1/2)/n1) * 225! sin(pi*j2*(k2+1/2)/n2), 226! 0<=k1<n1, 0<=k2<n2 227! <case2> DST 228! S(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 229! sin(pi*(j1+1/2)*k1/n1) * 230! sin(pi*(j2+1/2)*k2/n2), 231! 0<k1<=n1, 0<k2<=n2 232! [usage] 233! <case1> 234! ip(0) = 0 ! first time only 235! call ddst2d(n1max, n1, n2, 1, a, t, ip, w) 236! <case2> 237! ip(0) = 0 ! first time only 238! call ddst2d(n1max, n1, n2, -1, a, t, ip, w) 239! [parameters] 240! n1max :row size of the 2D array (integer) 241! n1 :data length (integer) 242! n1 >= 2, n1 = power of 2 243! n2 :data length (integer) 244! n2 >= 2, n2 = power of 2 245! a(0:n1-1,0:n2-1) 246! :input/output data (real*8) 247! <case1> 248! input data 249! a(j1,j2) = A(j1,j2), 0<j1<n1, 0<j2<n2, 250! a(j1,0) = A(j1,n2), 0<j1<n1, 251! a(0,j2) = A(n1,j2), 0<j2<n2, 252! a(0,0) = A(n1,n2) 253! (i.e. A(j1,j2) = a(mod(j1,n1),mod(j2,n2))) 254! output data 255! a(k1,k2) = S(k1,k2), 0<=k1<n1, 0<=k2<n2 256! <case2> 257! output data 258! a(k1,k2) = S(k1,k2), 0<k1<n1, 0<k2<n2, 259! a(k1,0) = S(k1,n2), 0<k1<n1, 260! a(0,k2) = S(n1,k2), 0<k2<n2, 261! a(0,0) = S(n1,n2) 262! (i.e. S(k1,k2) = a(mod(k1,n1),mod(k2,n2))) 263! t(0:n1-1,0:n2-1) 264! :work area (real*8) 265! ip(0:*):work area for bit reversal (integer) 266! length of ip >= 2+sqrt(n) 267! (n = max(n1/2, n2)) 268! ip(0),ip(1) are pointers of the cos/sin table. 269! w(0:*) :cos/sin table (real*8) 270! length of w >= max(n1/4, n2/2) + max(n1, n2) 271! w(),ip() are initialized if ip(0) = 0. 272! [remark] 273! Inverse of 274! call ddst2d(n1max, n1, n2, -1, a, t, ip, w) 275! is 276! do j1 = 0, n1 - 1 277! a(j1, 0) = a(j1, 0) * 0.5d0 278! end do 279! do j2 = 0, n2 - 1 280! a(0, j2) = a(0, j2) * 0.5d0 281! end do 282! call ddst2d(n1max, n1, n2, 1, a, t, ip, w) 283! do j2 = 0, n2 - 1 284! do j1 = 0, n1 - 1 285! a(j1, j2) = a(j1, j2) * (4.0d0 / (n1 * n2)) 286! end do 287! end do 288! . 289! 290! 291 subroutine cdft2d(n1max, n1, n2, isgn, a, ip, w) 292 integer n1max, n1, n2, isgn, ip(0 : *), n 293 real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *) 294 n = max(n1, 2 * n2) 295 if (n .gt. 4 * ip(0)) then 296 call makewt(n / 4, ip, w) 297 end if 298 if (n1 .gt. 4) then 299 call bitrv2row(n1max, n1, n2, ip(2), a) 300 end if 301 if (n2 .gt. 2) then 302 call bitrv2col(n1max, n1, n2, ip(2), a) 303 end if 304 if (isgn .lt. 0) then 305 call cftfrow(n1max, n1, n2, a, w) 306 call cftfcol(n1max, n1, n2, a, w) 307 else 308 call cftbrow(n1max, n1, n2, a, w) 309 call cftbcol(n1max, n1, n2, a, w) 310 end if 311 end 312! 313 subroutine rdft2d(n1max, n1, n2, isgn, a, ip, w) 314 integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, n2h, i, j 315 real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *), xi 316 n = max(n1, 2 * n2) 317 nw = ip(0) 318 if (n .gt. 4 * nw) then 319 nw = n / 4 320 call makewt(nw, ip, w) 321 end if 322 nc = ip(1) 323 if (n1 .gt. 4 * nc) then 324 nc = n1 / 4 325 call makect(nc, ip, w(nw)) 326 end if 327 n2h = n2 / 2 328 if (isgn .lt. 0) then 329 do i = 1, n2h - 1 330 j = n2 - i 331 xi = a(0, i) - a(0, j) 332 a(0, i) = a(0, i) + a(0, j) 333 a(0, j) = xi 334 xi = a(1, j) - a(1, i) 335 a(1, i) = a(1, i) + a(1, j) 336 a(1, j) = xi 337 end do 338 if (n2 .gt. 2) then 339 call bitrv2col(n1max, n1, n2, ip(2), a) 340 end if 341 call cftfcol(n1max, n1, n2, a, w) 342 do i = 0, n2 - 1 343 a(1, i) = 0.5d0 * (a(0, i) - a(1, i)) 344 a(0, i) = a(0, i) - a(1, i) 345 end do 346 if (n1 .gt. 4) then 347 call rftfrow(n1max, n1, n2, a, nc, w(nw)) 348 call bitrv2row(n1max, n1, n2, ip(2), a) 349 end if 350 call cftfrow(n1max, n1, n2, a, w) 351 else 352 if (n1 .gt. 4) then 353 call bitrv2row(n1max, n1, n2, ip(2), a) 354 end if 355 call cftbrow(n1max, n1, n2, a, w) 356 if (n1 .gt. 4) then 357 call rftbrow(n1max, n1, n2, a, nc, w(nw)) 358 end if 359 do i = 0, n2 - 1 360 xi = a(0, i) - a(1, i) 361 a(0, i) = a(0, i) + a(1, i) 362 a(1, i) = xi 363 end do 364 if (n2 .gt. 2) then 365 call bitrv2col(n1max, n1, n2, ip(2), a) 366 end if 367 call cftbcol(n1max, n1, n2, a, w) 368 do i = 1, n2h - 1 369 j = n2 - i 370 a(0, j) = 0.5d0 * (a(0, i) - a(0, j)) 371 a(0, i) = a(0, i) - a(0, j) 372 a(1, j) = 0.5d0 * (a(1, i) + a(1, j)) 373 a(1, i) = a(1, i) - a(1, j) 374 end do 375 end if 376 end 377! 378 subroutine ddct2d(n1max, n1, n2, isgn, a, t, ip, w) 379 integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, n1h, n2h, 380 & i, ix, ic, j, jx, jc 381 real*8 a(0 : n1max - 1, 0 : n2 - 1), 382 & t(0 : n1max - 1, 0 : n2 - 1), w(0 : *), xi 383 n = max(n1, 2 * n2) 384 nw = ip(0) 385 if (n .gt. 4 * nw) then 386 nw = n / 4 387 call makewt(nw, ip, w) 388 end if 389 nc = ip(1) 390 if (n1 .gt. nc .or. n2 .gt. nc) then 391 nc = max(n1, n2) 392 call makect(nc, ip, w(nw)) 393 end if 394 n1h = n1 / 2 395 n2h = n2 / 2 396 if (isgn .ge. 0) then 397 do i = 0, n2 - 1 398 do j = 1, n1h - 1 399 jx = 2 * j 400 t(jx, i) = a(j, i) 401 t(jx + 1, i) = a(n1 - j, i) 402 end do 403 end do 404 t(0, 0) = a(0, 0) 405 t(1, 0) = a(n1h, 0) 406 t(0, n2h) = a(0, n2h) 407 t(1, n2h) = a(n1h, n2h) 408 do i = 1, n2h - 1 409 ic = n2 - i 410 t(0, i) = a(0, i) 411 t(1, ic) = a(n1h, i) 412 t(1, i) = a(0, ic) 413 t(0, ic) = a(n1h, ic) 414 end do 415 call dctfsub(n1max, n1, n2, t, nc, w(nw)) 416 if (n2 .gt. 2) then 417 call bitrv2col(n1max, n1, n2, ip(2), t) 418 end if 419 call cftfcol(n1max, n1, n2, t, w) 420 do i = 0, n2 - 1 421 t(1, i) = 0.5d0 * (t(0, i) - t(1, i)) 422 t(0, i) = t(0, i) - t(1, i) 423 end do 424 if (n1 .gt. 4) then 425 call rftfrow(n1max, n1, n2, t, nc, w(nw)) 426 call bitrv2row(n1max, n1, n2, ip(2), t) 427 end if 428 call cftfrow(n1max, n1, n2, t, w) 429 do i = 0, n2h - 1 430 ix = 2 * i 431 ic = n2 - 1 - i 432 do j = 0, n1h - 1 433 jx = 2 * j 434 jc = n1 - 1 - j 435 a(jx, ix) = t(j, i) 436 a(jx + 1, ix) = t(jc, i) 437 a(jx, ix + 1) = t(j, ic) 438 a(jx + 1, ix + 1) = t(jc, ic) 439 end do 440 end do 441 else 442 do i = 0, n2h - 1 443 ix = 2 * i 444 ic = n2 - 1 - i 445 do j = 0, n1h - 1 446 jx = 2 * j 447 jc = n1 - 1 - j 448 t(j, i) = a(jx, ix) 449 t(jc, i) = a(jx + 1, ix) 450 t(j, ic) = a(jx, ix + 1) 451 t(jc, ic) = a(jx + 1, ix + 1) 452 end do 453 end do 454 if (n1 .gt. 4) then 455 call bitrv2row(n1max, n1, n2, ip(2), t) 456 end if 457 call cftbrow(n1max, n1, n2, t, w) 458 if (n1 .gt. 4) then 459 call rftbrow(n1max, n1, n2, t, nc, w(nw)) 460 end if 461 do i = 0, n2 - 1 462 xi = t(0, i) - t(1, i) 463 t(0, i) = t(0, i) + t(1, i) 464 t(1, i) = xi 465 end do 466 if (n2 .gt. 2) then 467 call bitrv2col(n1max, n1, n2, ip(2), t) 468 end if 469 call cftbcol(n1max, n1, n2, t, w) 470 call dctbsub(n1max, n1, n2, t, nc, w(nw)) 471 do i = 0, n2 - 1 472 do j = 1, n1h - 1 473 jx = 2 * j 474 a(j, i) = t(jx, i) 475 a(n1 - j, i) = t(jx + 1, i) 476 end do 477 end do 478 a(0, 0) = t(0, 0) 479 a(n1h, 0) = t(1, 0) 480 a(0, n2h) = t(0, n2h) 481 a(n1h, n2h) = t(1, n2h) 482 do i = 1, n2h - 1 483 ic = n2 - i 484 a(0, i) = t(0, i) 485 a(n1h, i) = t(1, ic) 486 a(0, ic) = t(1, i) 487 a(n1h, ic) = t(0, ic) 488 end do 489 end if 490 end 491! 492 subroutine ddst2d(n1max, n1, n2, isgn, a, t, ip, w) 493 integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, n1h, n2h, 494 & i, ix, ic, j, jx, jc 495 real*8 a(0 : n1max - 1, 0 : n2 - 1), 496 & t(0 : n1max - 1, 0 : n2 - 1), w(0 : *), xi 497 n = max(n1, 2 * n2) 498 nw = ip(0) 499 if (n .gt. 4 * nw) then 500 nw = n / 4 501 call makewt(nw, ip, w) 502 end if 503 nc = ip(1) 504 if (n1 .gt. nc .or. n2 .gt. nc) then 505 nc = max(n1, n2) 506 call makect(nc, ip, w(nw)) 507 end if 508 n1h = n1 / 2 509 n2h = n2 / 2 510 if (isgn .ge. 0) then 511 do i = 0, n2 - 1 512 do j = 1, n1h - 1 513 jx = 2 * j 514 t(jx, i) = a(j, i) 515 t(jx + 1, i) = a(n1 - j, i) 516 end do 517 end do 518 t(0, 0) = a(0, 0) 519 t(1, 0) = a(n1h, 0) 520 t(0, n2h) = a(0, n2h) 521 t(1, n2h) = a(n1h, n2h) 522 do i = 1, n2h - 1 523 ic = n2 - i 524 t(0, i) = a(0, i) 525 t(1, ic) = a(n1h, i) 526 t(1, i) = a(0, ic) 527 t(0, ic) = a(n1h, ic) 528 end do 529 call dstfsub(n1max, n1, n2, t, nc, w(nw)) 530 if (n2 .gt. 2) then 531 call bitrv2col(n1max, n1, n2, ip(2), t) 532 end if 533 call cftfcol(n1max, n1, n2, t, w) 534 do i = 0, n2 - 1 535 t(1, i) = 0.5d0 * (t(0, i) - t(1, i)) 536 t(0, i) = t(0, i) - t(1, i) 537 end do 538 if (n1 .gt. 4) then 539 call rftfrow(n1max, n1, n2, t, nc, w(nw)) 540 call bitrv2row(n1max, n1, n2, ip(2), t) 541 end if 542 call cftfrow(n1max, n1, n2, t, w) 543 do i = 0, n2h - 1 544 ix = 2 * i 545 ic = n2 - 1 - i 546 do j = 0, n1h - 1 547 jx = 2 * j 548 jc = n1 - 1 - j 549 a(jx, ix) = t(j, i) 550 a(jx + 1, ix) = -t(jc, i) 551 a(jx, ix + 1) = -t(j, ic) 552 a(jx + 1, ix + 1) = t(jc, ic) 553 end do 554 end do 555 else 556 do i = 0, n2h - 1 557 ix = 2 * i 558 ic = n2 - 1 - i 559 do j = 0, n1h - 1 560 jx = 2 * j 561 jc = n1 - 1 - j 562 t(j, i) = a(jx, ix) 563 t(jc, i) = -a(jx + 1, ix) 564 t(j, ic) = -a(jx, ix + 1) 565 t(jc, ic) = a(jx + 1, ix + 1) 566 end do 567 end do 568 if (n1 .gt. 4) then 569 call bitrv2row(n1max, n1, n2, ip(2), t) 570 end if 571 call cftbrow(n1max, n1, n2, t, w) 572 if (n1 .gt. 4) then 573 call rftbrow(n1max, n1, n2, t, nc, w(nw)) 574 end if 575 do i = 0, n2 - 1 576 xi = t(0, i) - t(1, i) 577 t(0, i) = t(0, i) + t(1, i) 578 t(1, i) = xi 579 end do 580 if (n2 .gt. 2) then 581 call bitrv2col(n1max, n1, n2, ip(2), t) 582 end if 583 call cftbcol(n1max, n1, n2, t, w) 584 call dstbsub(n1max, n1, n2, t, nc, w(nw)) 585 do i = 0, n2 - 1 586 do j = 1, n1h - 1 587 jx = 2 * j 588 a(j, i) = t(jx, i) 589 a(n1 - j, i) = t(jx + 1, i) 590 end do 591 end do 592 a(0, 0) = t(0, 0) 593 a(n1h, 0) = t(1, 0) 594 a(0, n2h) = t(0, n2h) 595 a(n1h, n2h) = t(1, n2h) 596 do i = 1, n2h - 1 597 ic = n2 - i 598 a(0, i) = t(0, i) 599 a(n1h, i) = t(1, ic) 600 a(0, ic) = t(1, i) 601 a(n1h, ic) = t(0, ic) 602 end do 603 end if 604 end 605! 606! -------- initializing routines -------- 607! 608 subroutine makewt(nw, ip, w) 609 integer nw, ip(0 : *), nwh, j 610 real*8 w(0 : nw - 1), delta, x, y 611 ip(0) = nw 612 ip(1) = 1 613 if (nw .gt. 2) then 614 nwh = nw / 2 615 delta = atan(1.0d0) / nwh 616 w(0) = 1 617 w(1) = 0 618 w(nwh) = cos(delta * nwh) 619 w(nwh + 1) = w(nwh) 620 do j = 2, nwh - 2, 2 621 x = cos(delta * j) 622 y = sin(delta * j) 623 w(j) = x 624 w(j + 1) = y 625 w(nw - j) = y 626 w(nw - j + 1) = x 627 end do 628 call bitrv2(nw, ip(2), w) 629 end if 630 end 631! 632 subroutine makect(nc, ip, c) 633 integer nc, ip(0 : *), nch, j 634 real*8 c(0 : nc - 1), delta 635 ip(1) = nc 636 if (nc .gt. 1) then 637 nch = nc / 2 638 delta = atan(1.0d0) / nch 639 c(0) = 0.5d0 640 c(nch) = 0.5d0 * cos(delta * nch) 641 do j = 1, nch - 1 642 c(j) = 0.5d0 * cos(delta * j) 643 c(nc - j) = 0.5d0 * sin(delta * j) 644 end do 645 end if 646 end 647! 648! -------- child routines -------- 649! 650 subroutine bitrv2(n, ip, a) 651 integer n, ip(0 : *), j, j1, k, k1, l, m, m2 652 real*8 a(0 : n - 1), xr, xi 653 ip(0) = 0 654 l = n 655 m = 1 656 do while (4 * m .lt. l) 657 l = l / 2 658 do j = 0, m - 1 659 ip(m + j) = ip(j) + l 660 end do 661 m = m * 2 662 end do 663 if (4 * m .gt. l) then 664 do k = 1, m - 1 665 do j = 0, k - 1 666 j1 = 2 * j + ip(k) 667 k1 = 2 * k + ip(j) 668 xr = a(j1) 669 xi = a(j1 + 1) 670 a(j1) = a(k1) 671 a(j1 + 1) = a(k1 + 1) 672 a(k1) = xr 673 a(k1 + 1) = xi 674 end do 675 end do 676 else 677 m2 = 2 * m 678 do k = 1, m - 1 679 do j = 0, k - 1 680 j1 = 2 * j + ip(k) 681 k1 = 2 * k + ip(j) 682 xr = a(j1) 683 xi = a(j1 + 1) 684 a(j1) = a(k1) 685 a(j1 + 1) = a(k1 + 1) 686 a(k1) = xr 687 a(k1 + 1) = xi 688 j1 = j1 + m2 689 k1 = k1 + m2 690 xr = a(j1) 691 xi = a(j1 + 1) 692 a(j1) = a(k1) 693 a(j1 + 1) = a(k1 + 1) 694 a(k1) = xr 695 a(k1 + 1) = xi 696 end do 697 end do 698 end if 699 end 700! 701 subroutine bitrv2row(n1max, n, n2, ip, a) 702 integer n1max, n, n2, ip(0 : *), i, j, j1, k, k1, l, m, m2 703 real*8 a(0 : n1max - 1, 0 : n2 - 1), xr, xi 704 ip(0) = 0 705 l = n 706 m = 1 707 do while (4 * m .lt. l) 708 l = l / 2 709 do j = 0, m - 1 710 ip(m + j) = ip(j) + l 711 end do 712 m = m * 2 713 end do 714 if (4 * m .gt. l) then 715 do i = 0, n2 - 1 716 do k = 1, m - 1 717 do j = 0, k - 1 718 j1 = 2 * j + ip(k) 719 k1 = 2 * k + ip(j) 720 xr = a(j1, i) 721 xi = a(j1 + 1, i) 722 a(j1, i) = a(k1, i) 723 a(j1 + 1, i) = a(k1 + 1, i) 724 a(k1, i) = xr 725 a(k1 + 1, i) = xi 726 end do 727 end do 728 end do 729 else 730 m2 = 2 * m 731 do i = 0, n2 - 1 732 do k = 1, m - 1 733 do j = 0, k - 1 734 j1 = 2 * j + ip(k) 735 k1 = 2 * k + ip(j) 736 xr = a(j1, i) 737 xi = a(j1 + 1, i) 738 a(j1, i) = a(k1, i) 739 a(j1 + 1, i) = a(k1 + 1, i) 740 a(k1, i) = xr 741 a(k1 + 1, i) = xi 742 j1 = j1 + m2 743 k1 = k1 + m2 744 xr = a(j1, i) 745 xi = a(j1 + 1, i) 746 a(j1, i) = a(k1, i) 747 a(j1 + 1, i) = a(k1 + 1, i) 748 a(k1, i) = xr 749 a(k1 + 1, i) = xi 750 end do 751 end do 752 end do 753 end if 754 end 755! 756 subroutine bitrv2col(n1max, n1, n, ip, a) 757 integer n1max, n1, n, ip(0 : *), i, j, j1, k, k1, l, m 758 real*8 a(0 : n1max - 1, 0 : n - 1), xr, xi 759 ip(0) = 0 760 l = n 761 m = 1 762 do while (2 * m .lt. l) 763 l = l / 2 764 do j = 0, m - 1 765 ip(m + j) = ip(j) + l 766 end do 767 m = m * 2 768 end do 769 if (2 * m .gt. l) then 770 do k = 1, m - 1 771 do j = 0, k - 1 772 j1 = j + ip(k) 773 k1 = k + ip(j) 774 do i = 0, n1 - 2, 2 775 xr = a(i, j1) 776 xi = a(i + 1, j1) 777 a(i, j1) = a(i, k1) 778 a(i + 1, j1) = a(i + 1, k1) 779 a(i, k1) = xr 780 a(i + 1, k1) = xi 781 end do 782 end do 783 end do 784 else 785 do k = 1, m - 1 786 do j = 0, k - 1 787 j1 = j + ip(k) 788 k1 = k + ip(j) 789 do i = 0, n1 - 2, 2 790 xr = a(i, j1) 791 xi = a(i + 1, j1) 792 a(i, j1) = a(i, k1) 793 a(i + 1, j1) = a(i + 1, k1) 794 a(i, k1) = xr 795 a(i + 1, k1) = xi 796 end do 797 j1 = j1 + m 798 k1 = k1 + m 799 do i = 0, n1 - 2, 2 800 xr = a(i, j1) 801 xi = a(i + 1, j1) 802 a(i, j1) = a(i, k1) 803 a(i + 1, j1) = a(i + 1, k1) 804 a(i, k1) = xr 805 a(i + 1, k1) = xi 806 end do 807 end do 808 end do 809 end if 810 end 811! 812 subroutine cftbrow(n1max, n, n2, a, w) 813 integer n1max, n, n2, i, j, j1, j2, j3, k, k1, ks, l, m 814 real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *) 815 real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i 816 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 817 do i = 0, n2 - 1 818 l = 2 819 do while (2 * l .lt. n) 820 m = 4 * l 821 do j = 0, l - 2, 2 822 j1 = j + l 823 j2 = j1 + l 824 j3 = j2 + l 825 x0r = a(j, i) + a(j1, i) 826 x0i = a(j + 1, i) + a(j1 + 1, i) 827 x1r = a(j, i) - a(j1, i) 828 x1i = a(j + 1, i) - a(j1 + 1, i) 829 x2r = a(j2, i) + a(j3, i) 830 x2i = a(j2 + 1, i) + a(j3 + 1, i) 831 x3r = a(j2, i) - a(j3, i) 832 x3i = a(j2 + 1, i) - a(j3 + 1, i) 833 a(j, i) = x0r + x2r 834 a(j + 1, i) = x0i + x2i 835 a(j2, i) = x0r - x2r 836 a(j2 + 1, i) = x0i - x2i 837 a(j1, i) = x1r - x3i 838 a(j1 + 1, i) = x1i + x3r 839 a(j3, i) = x1r + x3i 840 a(j3 + 1, i) = x1i - x3r 841 end do 842 if (m .lt. n) then 843 wk1r = w(2) 844 do j = m, l + m - 2, 2 845 j1 = j + l 846 j2 = j1 + l 847 j3 = j2 + l 848 x0r = a(j, i) + a(j1, i) 849 x0i = a(j + 1, i) + a(j1 + 1, i) 850 x1r = a(j, i) - a(j1, i) 851 x1i = a(j + 1, i) - a(j1 + 1, i) 852 x2r = a(j2, i) + a(j3, i) 853 x2i = a(j2 + 1, i) + a(j3 + 1, i) 854 x3r = a(j2, i) - a(j3, i) 855 x3i = a(j2 + 1, i) - a(j3 + 1, i) 856 a(j, i) = x0r + x2r 857 a(j + 1, i) = x0i + x2i 858 a(j2, i) = x2i - x0i 859 a(j2 + 1, i) = x0r - x2r 860 x0r = x1r - x3i 861 x0i = x1i + x3r 862 a(j1, i) = wk1r * (x0r - x0i) 863 a(j1 + 1, i) = wk1r * (x0r + x0i) 864 x0r = x3i + x1r 865 x0i = x3r - x1i 866 a(j3, i) = wk1r * (x0i - x0r) 867 a(j3 + 1, i) = wk1r * (x0i + x0r) 868 end do 869 k1 = 1 870 ks = -1 871 do k = 2 * m, n - m, m 872 k1 = k1 + 1 873 ks = -ks 874 wk1r = w(2 * k1) 875 wk1i = w(2 * k1 + 1) 876 wk2r = ks * w(k1) 877 wk2i = w(k1 + ks) 878 wk3r = wk1r - 2 * wk2i * wk1i 879 wk3i = 2 * wk2i * wk1r - wk1i 880 do j = k, l + k - 2, 2 881 j1 = j + l 882 j2 = j1 + l 883 j3 = j2 + l 884 x0r = a(j, i) + a(j1, i) 885 x0i = a(j + 1, i) + a(j1 + 1, i) 886 x1r = a(j, i) - a(j1, i) 887 x1i = a(j + 1, i) - a(j1 + 1, i) 888 x2r = a(j2, i) + a(j3, i) 889 x2i = a(j2 + 1, i) + a(j3 + 1, i) 890 x3r = a(j2, i) - a(j3, i) 891 x3i = a(j2 + 1, i) - a(j3 + 1, i) 892 a(j, i) = x0r + x2r 893 a(j + 1, i) = x0i + x2i 894 x0r = x0r - x2r 895 x0i = x0i - x2i 896 a(j2, i) = wk2r * x0r - wk2i * x0i 897 a(j2 + 1, i) = wk2r * x0i + wk2i * x0r 898 x0r = x1r - x3i 899 x0i = x1i + x3r 900 a(j1, i) = wk1r * x0r - wk1i * x0i 901 a(j1 + 1, i) = wk1r * x0i + wk1i * x0r 902 x0r = x1r + x3i 903 x0i = x1i - x3r 904 a(j3, i) = wk3r * x0r - wk3i * x0i 905 a(j3 + 1, i) = wk3r * x0i + wk3i * x0r 906 end do 907 end do 908 end if 909 l = m 910 end do 911 if (l .lt. n) then 912 do j = 0, l - 2, 2 913 j1 = j + l 914 x0r = a(j, i) - a(j1, i) 915 x0i = a(j + 1, i) - a(j1 + 1, i) 916 a(j, i) = a(j, i) + a(j1, i) 917 a(j + 1, i) = a(j + 1, i) + a(j1 + 1, i) 918 a(j1, i) = x0r 919 a(j1 + 1, i) = x0i 920 end do 921 end if 922 end do 923 end 924! 925 subroutine cftbcol(n1max, n1, n, a, w) 926 integer n1max, n1, n, i, j, j1, j2, j3, k, k1, ks, l, m 927 real*8 a(0 : n1max - 1, 0 : n - 1), w(0 : *) 928 real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i 929 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 930 l = 1 931 do while (2 * l .lt. n) 932 m = 4 * l 933 do j = 0, l - 1 934 j1 = j + l 935 j2 = j1 + l 936 j3 = j2 + l 937 do i = 0, n1 - 2, 2 938 x0r = a(i, j) + a(i, j1) 939 x0i = a(i + 1, j) + a(i + 1, j1) 940 x1r = a(i, j) - a(i, j1) 941 x1i = a(i + 1, j) - a(i + 1, j1) 942 x2r = a(i, j2) + a(i, j3) 943 x2i = a(i + 1, j2) + a(i + 1, j3) 944 x3r = a(i, j2) - a(i, j3) 945 x3i = a(i + 1, j2) - a(i + 1, j3) 946 a(i, j) = x0r + x2r 947 a(i + 1, j) = x0i + x2i 948 a(i, j2) = x0r - x2r 949 a(i + 1, j2) = x0i - x2i 950 a(i, j1) = x1r - x3i 951 a(i + 1, j1) = x1i + x3r 952 a(i, j3) = x1r + x3i 953 a(i + 1, j3) = x1i - x3r 954 end do 955 end do 956 if (m .lt. n) then 957 wk1r = w(2) 958 do j = m, l + m - 1 959 j1 = j + l 960 j2 = j1 + l 961 j3 = j2 + l 962 do i = 0, n1 - 2, 2 963 x0r = a(i, j) + a(i, j1) 964 x0i = a(i + 1, j) + a(i + 1, j1) 965 x1r = a(i, j) - a(i, j1) 966 x1i = a(i + 1, j) - a(i + 1, j1) 967 x2r = a(i, j2) + a(i, j3) 968 x2i = a(i + 1, j2) + a(i + 1, j3) 969 x3r = a(i, j2) - a(i, j3) 970 x3i = a(i + 1, j2) - a(i + 1, j3) 971 a(i, j) = x0r + x2r 972 a(i + 1, j) = x0i + x2i 973 a(i, j2) = x2i - x0i 974 a(i + 1, j2) = x0r - x2r 975 x0r = x1r - x3i 976 x0i = x1i + x3r 977 a(i, j1) = wk1r * (x0r - x0i) 978 a(i + 1, j1) = wk1r * (x0r + x0i) 979 x0r = x3i + x1r 980 x0i = x3r - x1i 981 a(i, j3) = wk1r * (x0i - x0r) 982 a(i + 1, j3) = wk1r * (x0i + x0r) 983 end do 984 end do 985 k1 = 1 986 ks = -1 987 do k = 2 * m, n - m, m 988 k1 = k1 + 1 989 ks = -ks 990 wk1r = w(2 * k1) 991 wk1i = w(2 * k1 + 1) 992 wk2r = ks * w(k1) 993 wk2i = w(k1 + ks) 994 wk3r = wk1r - 2 * wk2i * wk1i 995 wk3i = 2 * wk2i * wk1r - wk1i 996 do j = k, l + k - 1 997 j1 = j + l 998 j2 = j1 + l 999 j3 = j2 + l 1000 do i = 0, n1 - 2, 2 1001 x0r = a(i, j) + a(i, j1) 1002 x0i = a(i + 1, j) + a(i + 1, j1) 1003 x1r = a(i, j) - a(i, j1) 1004 x1i = a(i + 1, j) - a(i + 1, j1) 1005 x2r = a(i, j2) + a(i, j3) 1006 x2i = a(i + 1, j2) + a(i + 1, j3) 1007 x3r = a(i, j2) - a(i, j3) 1008 x3i = a(i + 1, j2) - a(i + 1, j3) 1009 a(i, j) = x0r + x2r 1010 a(i + 1, j) = x0i + x2i 1011 x0r = x0r - x2r 1012 x0i = x0i - x2i 1013 a(i, j2) = wk2r * x0r - wk2i * x0i 1014 a(i + 1, j2) = wk2r * x0i + wk2i * x0r 1015 x0r = x1r - x3i 1016 x0i = x1i + x3r 1017 a(i, j1) = wk1r * x0r - wk1i * x0i 1018 a(i + 1, j1) = wk1r * x0i + wk1i * x0r 1019 x0r = x1r + x3i 1020 x0i = x1i - x3r 1021 a(i, j3) = wk3r * x0r - wk3i * x0i 1022 a(i + 1, j3) = wk3r * x0i + wk3i * x0r 1023 end do 1024 end do 1025 end do 1026 end if 1027 l = m 1028 end do 1029 if (l .lt. n) then 1030 do j = 0, l - 1 1031 j1 = j + l 1032 do i = 0, n1 - 2, 2 1033 x0r = a(i, j) - a(i, j1) 1034 x0i = a(i + 1, j) - a(i + 1, j1) 1035 a(i, j) = a(i, j) + a(i, j1) 1036 a(i + 1, j) = a(i + 1, j) + a(i + 1, j1) 1037 a(i, j1) = x0r 1038 a(i + 1, j1) = x0i 1039 end do 1040 end do 1041 end if 1042 end 1043! 1044 subroutine cftfrow(n1max, n, n2, a, w) 1045 integer n1max, n, n2, i, j, j1, j2, j3, k, k1, ks, l, m 1046 real*8 a(0 : n1max - 1, 0 : n2 - 1), w(0 : *) 1047 real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i 1048 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 1049 do i = 0, n2 - 1 1050 l = 2 1051 do while (2 * l .lt. n) 1052 m = 4 * l 1053 do j = 0, l - 2, 2 1054 j1 = j + l 1055 j2 = j1 + l 1056 j3 = j2 + l 1057 x0r = a(j, i) + a(j1, i) 1058 x0i = a(j + 1, i) + a(j1 + 1, i) 1059 x1r = a(j, i) - a(j1, i) 1060 x1i = a(j + 1, i) - a(j1 + 1, i) 1061 x2r = a(j2, i) + a(j3, i) 1062 x2i = a(j2 + 1, i) + a(j3 + 1, i) 1063 x3r = a(j2, i) - a(j3, i) 1064 x3i = a(j2 + 1, i) - a(j3 + 1, i) 1065 a(j, i) = x0r + x2r 1066 a(j + 1, i) = x0i + x2i 1067 a(j2, i) = x0r - x2r 1068 a(j2 + 1, i) = x0i - x2i 1069 a(j1, i) = x1r + x3i 1070 a(j1 + 1, i) = x1i - x3r 1071 a(j3, i) = x1r - x3i 1072 a(j3 + 1, i) = x1i + x3r 1073 end do 1074 if (m .lt. n) then 1075 wk1r = w(2) 1076 do j = m, l + m - 2, 2 1077 j1 = j + l 1078 j2 = j1 + l 1079 j3 = j2 + l 1080 x0r = a(j, i) + a(j1, i) 1081 x0i = a(j + 1, i) + a(j1 + 1, i) 1082 x1r = a(j, i) - a(j1, i) 1083 x1i = a(j + 1, i) - a(j1 + 1, i) 1084 x2r = a(j2, i) + a(j3, i) 1085 x2i = a(j2 + 1, i) + a(j3 + 1, i) 1086 x3r = a(j2, i) - a(j3, i) 1087 x3i = a(j2 + 1, i) - a(j3 + 1, i) 1088 a(j, i) = x0r + x2r 1089 a(j + 1, i) = x0i + x2i 1090 a(j2, i) = x0i - x2i 1091 a(j2 + 1, i) = x2r - x0r 1092 x0r = x1r + x3i 1093 x0i = x1i - x3r 1094 a(j1, i) = wk1r * (x0i + x0r) 1095 a(j1 + 1, i) = wk1r * (x0i - x0r) 1096 x0r = x3i - x1r 1097 x0i = x3r + x1i 1098 a(j3, i) = wk1r * (x0r + x0i) 1099 a(j3 + 1, i) = wk1r * (x0r - x0i) 1100 end do 1101 k1 = 1 1102 ks = -1 1103 do k = 2 * m, n - m, m 1104 k1 = k1 + 1 1105 ks = -ks 1106 wk1r = w(2 * k1) 1107 wk1i = w(2 * k1 + 1) 1108 wk2r = ks * w(k1) 1109 wk2i = w(k1 + ks) 1110 wk3r = wk1r - 2 * wk2i * wk1i 1111 wk3i = 2 * wk2i * wk1r - wk1i 1112 do j = k, l + k - 2, 2 1113 j1 = j + l 1114 j2 = j1 + l 1115 j3 = j2 + l 1116 x0r = a(j, i) + a(j1, i) 1117 x0i = a(j + 1, i) + a(j1 + 1, i) 1118 x1r = a(j, i) - a(j1, i) 1119 x1i = a(j + 1, i) - a(j1 + 1, i) 1120 x2r = a(j2, i) + a(j3, i) 1121 x2i = a(j2 + 1, i) + a(j3 + 1, i) 1122 x3r = a(j2, i) - a(j3, i) 1123 x3i = a(j2 + 1, i) - a(j3 + 1, i) 1124 a(j, i) = x0r + x2r 1125 a(j + 1, i) = x0i + x2i 1126 x0r = x0r - x2r 1127 x0i = x0i - x2i 1128 a(j2, i) = wk2r * x0r + wk2i * x0i 1129 a(j2 + 1, i) = wk2r * x0i - wk2i * x0r 1130 x0r = x1r + x3i 1131 x0i = x1i - x3r 1132 a(j1, i) = wk1r * x0r + wk1i * x0i 1133 a(j1 + 1, i) = wk1r * x0i - wk1i * x0r 1134 x0r = x1r - x3i 1135 x0i = x1i + x3r 1136 a(j3, i) = wk3r * x0r + wk3i * x0i 1137 a(j3 + 1, i) = wk3r * x0i - wk3i * x0r 1138 end do 1139 end do 1140 end if 1141 l = m 1142 end do 1143 if (l .lt. n) then 1144 do j = 0, l - 2, 2 1145 j1 = j + l 1146 x0r = a(j, i) - a(j1, i) 1147 x0i = a(j + 1, i) - a(j1 + 1, i) 1148 a(j, i) = a(j, i) + a(j1, i) 1149 a(j + 1, i) = a(j + 1, i) + a(j1 + 1, i) 1150 a(j1, i) = x0r 1151 a(j1 + 1, i) = x0i 1152 end do 1153 end if 1154 end do 1155 end 1156! 1157 subroutine cftfcol(n1max, n1, n, a, w) 1158 integer n1max, n1, n, i, j, j1, j2, j3, k, k1, ks, l, m 1159 real*8 a(0 : n1max - 1, 0 : n - 1), w(0 : *) 1160 real*8 wk1r, wk1i, wk2r, wk2i, wk3r, wk3i 1161 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 1162 l = 1 1163 do while (2 * l .lt. n) 1164 m = 4 * l 1165 do j = 0, l - 1 1166 j1 = j + l 1167 j2 = j1 + l 1168 j3 = j2 + l 1169 do i = 0, n1 - 2, 2 1170 x0r = a(i, j) + a(i, j1) 1171 x0i = a(i + 1, j) + a(i + 1, j1) 1172 x1r = a(i, j) - a(i, j1) 1173 x1i = a(i + 1, j) - a(i + 1, j1) 1174 x2r = a(i, j2) + a(i, j3) 1175 x2i = a(i + 1, j2) + a(i + 1, j3) 1176 x3r = a(i, j2) - a(i, j3) 1177 x3i = a(i + 1, j2) - a(i + 1, j3) 1178 a(i, j) = x0r + x2r 1179 a(i + 1, j) = x0i + x2i 1180 a(i, j2) = x0r - x2r 1181 a(i + 1, j2) = x0i - x2i 1182 a(i, j1) = x1r + x3i 1183 a(i + 1, j1) = x1i - x3r 1184 a(i, j3) = x1r - x3i 1185 a(i + 1, j3) = x1i + x3r 1186 end do 1187 end do 1188 if (m .lt. n) then 1189 wk1r = w(2) 1190 do j = m, l + m - 1 1191 j1 = j + l 1192 j2 = j1 + l 1193 j3 = j2 + l 1194 do i = 0, n1 - 2, 2 1195 x0r = a(i, j) + a(i, j1) 1196 x0i = a(i + 1, j) + a(i + 1, j1) 1197 x1r = a(i, j) - a(i, j1) 1198 x1i = a(i + 1, j) - a(i + 1, j1) 1199 x2r = a(i, j2) + a(i, j3) 1200 x2i = a(i + 1, j2) + a(i + 1, j3) 1201 x3r = a(i, j2) - a(i, j3) 1202 x3i = a(i + 1, j2) - a(i + 1, j3) 1203 a(i, j) = x0r + x2r 1204 a(i + 1, j) = x0i + x2i 1205 a(i, j2) = x0i - x2i 1206 a(i + 1, j2) = x2r - x0r 1207 x0r = x1r + x3i 1208 x0i = x1i - x3r 1209 a(i, j1) = wk1r * (x0i + x0r) 1210 a(i + 1, j1) = wk1r * (x0i - x0r) 1211 x0r = x3i - x1r 1212 x0i = x3r + x1i 1213 a(i, j3) = wk1r * (x0r + x0i) 1214 a(i + 1, j3) = wk1r * (x0r - x0i) 1215 end do 1216 end do 1217 k1 = 1 1218 ks = -1 1219 do k = 2 * m, n - m, m 1220 k1 = k1 + 1 1221 ks = -ks 1222 wk1r = w(2 * k1) 1223 wk1i = w(2 * k1 + 1) 1224 wk2r = ks * w(k1) 1225 wk2i = w(k1 + ks) 1226 wk3r = wk1r - 2 * wk2i * wk1i 1227 wk3i = 2 * wk2i * wk1r - wk1i 1228 do j = k, l + k - 1 1229 j1 = j + l 1230 j2 = j1 + l 1231 j3 = j2 + l 1232 do i = 0, n1 - 2, 2 1233 x0r = a(i, j) + a(i, j1) 1234 x0i = a(i + 1, j) + a(i + 1, j1) 1235 x1r = a(i, j) - a(i, j1) 1236 x1i = a(i + 1, j) - a(i + 1, j1) 1237 x2r = a(i, j2) + a(i, j3) 1238 x2i = a(i + 1, j2) + a(i + 1, j3) 1239 x3r = a(i, j2) - a(i, j3) 1240 x3i = a(i + 1, j2) - a(i + 1, j3) 1241 a(i, j) = x0r + x2r 1242 a(i + 1, j) = x0i + x2i 1243 x0r = x0r - x2r 1244 x0i = x0i - x2i 1245 a(i, j2) = wk2r * x0r + wk2i * x0i 1246 a(i + 1, j2) = wk2r * x0i - wk2i * x0r 1247 x0r = x1r + x3i 1248 x0i = x1i - x3r 1249 a(i, j1) = wk1r * x0r + wk1i * x0i 1250 a(i + 1, j1) = wk1r * x0i - wk1i * x0r 1251 x0r = x1r - x3i 1252 x0i = x1i + x3r 1253 a(i, j3) = wk3r * x0r + wk3i * x0i 1254 a(i + 1, j3) = wk3r * x0i - wk3i * x0r 1255 end do 1256 end do 1257 end do 1258 end if 1259 l = m 1260 end do 1261 if (l .lt. n) then 1262 do j = 0, l - 1 1263 j1 = j + l 1264 do i = 0, n1 - 2, 2 1265 x0r = a(i, j) - a(i, j1) 1266 x0i = a(i + 1, j) - a(i + 1, j1) 1267 a(i, j) = a(i, j) + a(i, j1) 1268 a(i + 1, j) = a(i + 1, j) + a(i + 1, j1) 1269 a(i, j1) = x0r 1270 a(i + 1, j1) = x0i 1271 end do 1272 end do 1273 end if 1274 end 1275! 1276 subroutine rftbrow(n1max, n, n2, a, nc, c) 1277 integer n1max, n, n2, nc, i, j, k, kk, ks 1278 real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), 1279 & wkr, wki, xr, xi, yr, yi 1280 ks = 4 * nc / n 1281 do i = 0, n2 - 1 1282 kk = 0 1283 do k = n / 2 - 2, 2, -2 1284 j = n - k 1285 kk = kk + ks 1286 wkr = 0.5d0 - c(kk) 1287 wki = c(nc - kk) 1288 xr = a(k, i) - a(j, i) 1289 xi = a(k + 1, i) + a(j + 1, i) 1290 yr = wkr * xr - wki * xi 1291 yi = wkr * xi + wki * xr 1292 a(k, i) = a(k, i) - yr 1293 a(k + 1, i) = a(k + 1, i) - yi 1294 a(j, i) = a(j, i) + yr 1295 a(j + 1, i) = a(j + 1, i) - yi 1296 end do 1297 end do 1298 end 1299! 1300 subroutine rftfrow(n1max, n, n2, a, nc, c) 1301 integer n1max, n, n2, nc, i, j, k, kk, ks 1302 real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), 1303 & wkr, wki, xr, xi, yr, yi 1304 ks = 4 * nc / n 1305 do i = 0, n2 - 1 1306 kk = 0 1307 do k = n / 2 - 2, 2, -2 1308 j = n - k 1309 kk = kk + ks 1310 wkr = 0.5d0 - c(kk) 1311 wki = c(nc - kk) 1312 xr = a(k, i) - a(j, i) 1313 xi = a(k + 1, i) + a(j + 1, i) 1314 yr = wkr * xr + wki * xi 1315 yi = wkr * xi - wki * xr 1316 a(k, i) = a(k, i) - yr 1317 a(k + 1, i) = a(k + 1, i) - yi 1318 a(j, i) = a(j, i) + yr 1319 a(j + 1, i) = a(j + 1, i) - yi 1320 end do 1321 end do 1322 end 1323! 1324 subroutine dctbsub(n1max, n1, n2, a, nc, c) 1325 integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, 1326 & k1, k2 1327 real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), 1328 & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i 1329 ks1 = nc / n1 1330 ks2 = nc / n2 1331 n2h = n2 / 2 1332 kk2 = ks2 1333 do k2 = 1, n2h - 1 1334 j2 = n2 - k2 1335 w2r = 2 * c(kk2) 1336 w2i = 2 * c(nc - kk2) 1337 kk2 = kk2 + ks2 1338 kk1 = ks1 1339 do k1 = 2, n1 - 2, 2 1340 x0r = w2r * c(kk1) 1341 x0i = w2i * c(kk1) 1342 x1r = w2r * c(nc - kk1) 1343 x1i = w2i * c(nc - kk1) 1344 wkr = x0r - x1i 1345 wki = x0i + x1r 1346 wji = x0r + x1i 1347 wjr = x0i - x1r 1348 kk1 = kk1 + ks1 1349 x0r = wkr * a(k1, k2) - wki * a(k1 + 1, k2) 1350 x0i = wkr * a(k1 + 1, k2) + wki * a(k1, k2) 1351 x1r = wjr * a(k1, j2) - wji * a(k1 + 1, j2) 1352 x1i = wjr * a(k1 + 1, j2) + wji * a(k1, j2) 1353 a(k1, k2) = x0r + x1i 1354 a(k1 + 1, k2) = x0i - x1r 1355 a(k1, j2) = x1r + x0i 1356 a(k1 + 1, j2) = x1i - x0r 1357 end do 1358 wkr = w2r * 0.5d0 1359 wki = w2i * 0.5d0 1360 wjr = w2r * c(kk1) 1361 wji = w2i * c(kk1) 1362 x0r = a(0, k2) + a(0, j2) 1363 x0i = a(1, k2) - a(1, j2) 1364 x1r = a(0, k2) - a(0, j2) 1365 x1i = a(1, k2) + a(1, j2) 1366 a(0, k2) = wkr * x0r - wki * x0i 1367 a(1, k2) = wkr * x0i + wki * x0r 1368 a(0, j2) = -wjr * x1r + wji * x1i 1369 a(1, j2) = wjr * x1i + wji * x1r 1370 end do 1371 w2r = 2 * c(kk2) 1372 kk1 = ks1 1373 do k1 = 2, n1 - 2, 2 1374 wkr = 2 * c(kk1) 1375 wki = 2 * c(nc - kk1) 1376 wjr = w2r * wkr 1377 wji = w2r * wki 1378 kk1 = kk1 + ks1 1379 x0i = wkr * a(k1 + 1, 0) + wki * a(k1, 0) 1380 a(k1, 0) = wkr * a(k1, 0) - wki * a(k1 + 1, 0) 1381 a(k1 + 1, 0) = x0i 1382 x0i = wjr * a(k1 + 1, n2h) + wji * a(k1, n2h) 1383 a(k1, n2h) = wjr * a(k1, n2h) - wji * a(k1 + 1, n2h) 1384 a(k1 + 1, n2h) = x0i 1385 end do 1386 a(1, 0) = a(1, 0) * w2r 1387 a(0, n2h) = a(0, n2h) * w2r 1388 a(1, n2h) = a(1, n2h) * 0.5d0 1389 end 1390! 1391 subroutine dctfsub(n1max, n1, n2, a, nc, c) 1392 integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, 1393 & k1, k2 1394 real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), 1395 & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i 1396 ks1 = nc / n1 1397 ks2 = nc / n2 1398 n2h = n2 / 2 1399 kk2 = ks2 1400 do k2 = 1, n2h - 1 1401 j2 = n2 - k2 1402 w2r = 2 * c(kk2) 1403 w2i = 2 * c(nc - kk2) 1404 kk2 = kk2 + ks2 1405 kk1 = ks1 1406 do k1 = 2, n1 - 2, 2 1407 x0r = w2r * c(kk1) 1408 x0i = w2i * c(kk1) 1409 x1r = w2r * c(nc - kk1) 1410 x1i = w2i * c(nc - kk1) 1411 wkr = x0r - x1i 1412 wki = x0i + x1r 1413 wji = x0r + x1i 1414 wjr = x0i - x1r 1415 kk1 = kk1 + ks1 1416 x0r = a(k1, k2) - a(k1 + 1, j2) 1417 x0i = a(k1, j2) + a(k1 + 1, k2) 1418 x1r = a(k1, j2) - a(k1 + 1, k2) 1419 x1i = a(k1, k2) + a(k1 + 1, j2) 1420 a(k1, k2) = wkr * x0r + wki * x0i 1421 a(k1 + 1, k2) = wkr * x0i - wki * x0r 1422 a(k1, j2) = wjr * x1r + wji * x1i 1423 a(k1 + 1, j2) = wjr * x1i - wji * x1r 1424 end do 1425 x0r = 2 * c(kk1) 1426 wjr = x0r * w2r 1427 wji = x0r * w2i 1428 x0r = w2r * a(0, k2) + w2i * a(1, k2) 1429 x0i = w2r * a(1, k2) - w2i * a(0, k2) 1430 x1r = -wjr * a(0, j2) + wji * a(1, j2) 1431 x1i = wjr * a(1, j2) + wji * a(0, j2) 1432 a(0, k2) = x0r + x1r 1433 a(1, k2) = x1i + x0i 1434 a(0, j2) = x0r - x1r 1435 a(1, j2) = x1i - x0i 1436 end do 1437 w2r = 2 * c(kk2) 1438 kk1 = ks1 1439 do k1 = 2, n1 - 2, 2 1440 wkr = 2 * c(kk1) 1441 wki = 2 * c(nc - kk1) 1442 wjr = w2r * wkr 1443 wji = w2r * wki 1444 kk1 = kk1 + ks1 1445 x0i = wkr * a(k1 + 1, 0) - wki * a(k1, 0) 1446 a(k1, 0) = wkr * a(k1, 0) + wki * a(k1 + 1, 0) 1447 a(k1 + 1, 0) = x0i 1448 x0i = wjr * a(k1 + 1, n2h) - wji * a(k1, n2h) 1449 a(k1, n2h) = wjr * a(k1, n2h) + wji * a(k1 + 1, n2h) 1450 a(k1 + 1, n2h) = x0i 1451 end do 1452 w2r = w2r * 2 1453 a(0, 0) = a(0, 0) * 2 1454 a(1, 0) = a(1, 0) * w2r 1455 a(0, n2h) = a(0, n2h) * w2r 1456 end 1457! 1458 subroutine dstbsub(n1max, n1, n2, a, nc, c) 1459 integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, 1460 & k1, k2 1461 real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), 1462 & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i 1463 ks1 = nc / n1 1464 ks2 = nc / n2 1465 n2h = n2 / 2 1466 kk2 = ks2 1467 do k2 = 1, n2h - 1 1468 j2 = n2 - k2 1469 w2r = 2 * c(kk2) 1470 w2i = 2 * c(nc - kk2) 1471 kk2 = kk2 + ks2 1472 kk1 = ks1 1473 do k1 = 2, n1 - 2, 2 1474 x0r = w2r * c(kk1) 1475 x0i = w2i * c(kk1) 1476 x1r = w2r * c(nc - kk1) 1477 x1i = w2i * c(nc - kk1) 1478 wkr = x0r - x1i 1479 wki = x0i + x1r 1480 wji = x0r + x1i 1481 wjr = x0i - x1r 1482 kk1 = kk1 + ks1 1483 x0r = wkr * a(k1, k2) - wki * a(k1 + 1, k2) 1484 x0i = wkr * a(k1 + 1, k2) + wki * a(k1, k2) 1485 x1r = wjr * a(k1, j2) - wji * a(k1 + 1, j2) 1486 x1i = wjr * a(k1 + 1, j2) + wji * a(k1, j2) 1487 a(k1, k2) = x1i - x0r 1488 a(k1 + 1, k2) = x1r + x0i 1489 a(k1, j2) = x0i - x1r 1490 a(k1 + 1, j2) = x0r + x1i 1491 end do 1492 wkr = w2r * 0.5d0 1493 wki = w2i * 0.5d0 1494 wjr = w2r * c(kk1) 1495 wji = w2i * c(kk1) 1496 x0r = a(0, k2) + a(0, j2) 1497 x0i = a(1, k2) - a(1, j2) 1498 x1r = a(0, k2) - a(0, j2) 1499 x1i = a(1, k2) + a(1, j2) 1500 a(1, k2) = wkr * x0r - wki * x0i 1501 a(0, k2) = wkr * x0i + wki * x0r 1502 a(1, j2) = -wjr * x1r + wji * x1i 1503 a(0, j2) = wjr * x1i + wji * x1r 1504 end do 1505 w2r = 2 * c(kk2) 1506 kk1 = ks1 1507 do k1 = 2, n1 - 2, 2 1508 wkr = 2 * c(kk1) 1509 wki = 2 * c(nc - kk1) 1510 wjr = w2r * wkr 1511 wji = w2r * wki 1512 kk1 = kk1 + ks1 1513 x0i = wkr * a(k1 + 1, 0) + wki * a(k1, 0) 1514 a(k1 + 1, 0) = wkr * a(k1, 0) - wki * a(k1 + 1, 0) 1515 a(k1, 0) = x0i 1516 x0i = wjr * a(k1 + 1, n2h) + wji * a(k1, n2h) 1517 a(k1 + 1, n2h) = wjr * a(k1, n2h) - wji * a(k1 + 1, n2h) 1518 a(k1, n2h) = x0i 1519 end do 1520 a(1, 0) = a(1, 0) * w2r 1521 a(0, n2h) = a(0, n2h) * w2r 1522 a(1, n2h) = a(1, n2h) * 0.5d0 1523 end 1524! 1525 subroutine dstfsub(n1max, n1, n2, a, nc, c) 1526 integer n1max, n1, n2, nc, kk1, kk2, ks1, ks2, n2h, j2, 1527 & k1, k2 1528 real*8 a(0 : n1max - 1, 0 : n2 - 1), c(0 : nc - 1), 1529 & w2r, w2i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i 1530 ks1 = nc / n1 1531 ks2 = nc / n2 1532 n2h = n2 / 2 1533 kk2 = ks2 1534 do k2 = 1, n2h - 1 1535 j2 = n2 - k2 1536 w2r = 2 * c(kk2) 1537 w2i = 2 * c(nc - kk2) 1538 kk2 = kk2 + ks2 1539 kk1 = ks1 1540 do k1 = 2, n1 - 2, 2 1541 x0r = w2r * c(kk1) 1542 x0i = w2i * c(kk1) 1543 x1r = w2r * c(nc - kk1) 1544 x1i = w2i * c(nc - kk1) 1545 wkr = x0r - x1i 1546 wki = x0i + x1r 1547 wji = x0r + x1i 1548 wjr = x0i - x1r 1549 kk1 = kk1 + ks1 1550 x0r = a(k1 + 1, j2) - a(k1, k2) 1551 x0i = a(k1 + 1, k2) + a(k1, j2) 1552 x1r = a(k1 + 1, k2) - a(k1, j2) 1553 x1i = a(k1 + 1, j2) + a(k1, k2) 1554 a(k1, k2) = wkr * x0r + wki * x0i 1555 a(k1 + 1, k2) = wkr * x0i - wki * x0r 1556 a(k1, j2) = wjr * x1r + wji * x1i 1557 a(k1 + 1, j2) = wjr * x1i - wji * x1r 1558 end do 1559 x0r = 2 * c(kk1) 1560 wjr = x0r * w2r 1561 wji = x0r * w2i 1562 x0r = w2r * a(1, k2) + w2i * a(0, k2) 1563 x0i = w2r * a(0, k2) - w2i * a(1, k2) 1564 x1r = -wjr * a(1, j2) + wji * a(0, j2) 1565 x1i = wjr * a(0, j2) + wji * a(1, j2) 1566 a(0, k2) = x0r + x1r 1567 a(1, k2) = x1i + x0i 1568 a(0, j2) = x0r - x1r 1569 a(1, j2) = x1i - x0i 1570 end do 1571 w2r = 2 * c(kk2) 1572 kk1 = ks1 1573 do k1 = 2, n1 - 2, 2 1574 wkr = 2 * c(kk1) 1575 wki = 2 * c(nc - kk1) 1576 wjr = w2r * wkr 1577 wji = w2r * wki 1578 kk1 = kk1 + ks1 1579 x0i = wkr * a(k1, 0) - wki * a(k1 + 1, 0) 1580 a(k1, 0) = wkr * a(k1 + 1, 0) + wki * a(k1, 0) 1581 a(k1 + 1, 0) = x0i 1582 x0i = wjr * a(k1, n2h) - wji * a(k1 + 1, n2h) 1583 a(k1, n2h) = wjr * a(k1 + 1, n2h) + wji * a(k1, n2h) 1584 a(k1 + 1, n2h) = x0i 1585 end do 1586 w2r = w2r * 2 1587 a(0, 0) = a(0, 0) * 2 1588 a(1, 0) = a(1, 0) * w2r 1589 a(0, n2h) = a(0, n2h) * w2r 1590 end 1591! 1592