1! Fast Fourier/Cosine/Sine Transform 2! dimension :three 3! data length :power of 2 4! decimation :frequency 5! radix :split-radix, row-column 6! data :inplace 7! table :use 8! subroutines 9! cdft3d: Complex Discrete Fourier Transform 10! rdft3d: Real Discrete Fourier Transform 11! ddct3d: Discrete Cosine Transform 12! ddst3d: Discrete Sine Transform 13! necessary package 14! fftsg.f : 1D-FFT package 15! 16! 17! -------- Complex DFT (Discrete Fourier Transform) -------- 18! [definition] 19! <case1> 20! X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 21! x(j1,j2,j3) * 22! exp(2*pi*i*j1*k1/n1) * 23! exp(2*pi*i*j2*k2/n2) * 24! exp(2*pi*i*j3*k3/n3), 25! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 26! <case2> 27! X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 28! x(j1,j2,j3) * 29! exp(-2*pi*i*j1*k1/n1) * 30! exp(-2*pi*i*j2*k2/n2) * 31! exp(-2*pi*i*j3*k3/n3), 32! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 33! (notes: sum_j=0^n-1 is a summation from j=0 to n-1) 34! [usage] 35! <case1> 36! ip(0) = 0 ! first time only 37! call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w) 38! <case2> 39! ip(0) = 0 ! first time only 40! call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w) 41! [parameters] 42! n1max :row1 size of the 3D array (integer) 43! n2max :row2 size of the 3D array (integer) 44! 2*n1 :data length (integer) 45! n1 >= 1, n1 = power of 2 46! n2 :data length (integer) 47! n2 >= 1, n2 = power of 2 48! n3 :data length (integer) 49! n3 >= 1, n3 = power of 2 50! a(0:2*n1-1,0:n2-1,0:n3-1) 51! :input/output data (real*8) 52! input data 53! a(2*j1,j2,j3) = Re(x(j1,j2,j3)), 54! a(2*j1+1,j2,j3) = Im(x(j1,j2,j3)), 55! 0<=j1<n1, 0<=j2<n2, 0<=j3<n3 56! output data 57! a(2*k1,k2,k3) = Re(X(k1,k2,k3)), 58! a(2*k1+1,k2,k3) = Im(X(k1,k2,k3)), 59! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 60! t(0:*) :work area (real*8) 61! length of t >= max(8*n2, 8*n3) 62! ip(0:*):work area for bit reversal (integer) 63! length of ip >= 2+sqrt(n) 64! (n = max(n1, n2, n3)) 65! ip(0),ip(1) are pointers of the cos/sin table. 66! w(0:*) :cos/sin table (real*8) 67! length of w >= max(n1/2, n2/2, n3/2) 68! w(),ip() are initialized if ip(0) = 0. 69! [remark] 70! Inverse of 71! call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w) 72! is 73! call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w) 74! do j3 = 0, n3 - 1 75! do j2 = 0, n2 - 1 76! do j1 = 0, 2 * n1 - 1 77! a(j1,j2,j3) = a(j1,j2,j3) * (1.0d0/n1/n2/n3) 78! end do 79! end do 80! end do 81! . 82! 83! 84! -------- Real DFT / Inverse of Real DFT -------- 85! [definition] 86! <case1> RDFT 87! R(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 88! a(j1,j2,j3) * 89! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 90! 2*pi*j3*k3/n3), 91! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 92! I(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 93! a(j1,j2,j3) * 94! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 95! 2*pi*j3*k3/n3), 96! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 97! <case2> IRDFT (excluding scale) 98! a(k1,k2,k3) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 99! (R(j1,j2,j3) * 100! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 101! 2*pi*j3*k3/n3) + 102! I(j1,j2,j3) * 103! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 104! 2*pi*j3*k3/n3)), 105! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 106! (notes: R(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = R(k1,k2,k3), 107! I(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = -I(k1,k2,k3), 108! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3) 109! [usage] 110! <case1> 111! ip(0) = 0 ! first time only 112! call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) 113! <case2> 114! ip(0) = 0 ! first time only 115! call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) 116! [parameters] 117! n1max :row1 size of the 3D array (integer) 118! n2max :row2 size of the 3D array (integer) 119! n1 :data length (integer) 120! n1 >= 2, n1 = power of 2 121! n2 :data length (integer) 122! n2 >= 2, n2 = power of 2 123! n3 :data length (integer) 124! n3 >= 2, n3 = power of 2 125! a(0:n1-1,0:n2-1,0:n3-1) 126! :input/output data (real*8) 127! <case1> 128! output data 129! a(2*k1,k2,k3) = R(k1,k2,k3) 130! = R(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), 131! a(2*k1+1,k2,k3) = I(k1,k2,k3) 132! = -I(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), 133! 0<k1<n1/2, 0<=k2<n2, 0<=k3<n3, 134! a(0,k2,k3) = R(0,k2,k3) 135! = R(0,n2-k2,mod(n3-k3,n3)), 136! a(1,k2,k3) = I(0,k2,k3) 137! = -I(0,n2-k2,mod(n3-k3,n3)), 138! a(1,n2-k2,k3) = R(n1/2,k2,mod(n3-k3,n3)) 139! = R(n1/2,n2-k2,k3), 140! a(0,n2-k2,k3) = -I(n1/2,k2,mod(n3-k3,n3)) 141! = I(n1/2,n2-k2,k3), 142! 0<k2<n2/2, 0<=k3<n3, 143! a(0,0,k3) = R(0,0,k3) 144! = R(0,0,n3-k3), 145! a(1,0,k3) = I(0,0,k3) 146! = -I(0,0,n3-k3), 147! a(0,n2/2,k3) = R(0,n2/2,k3) 148! = R(0,n2/2,n3-k3), 149! a(1,n2/2,k3) = I(0,n2/2,k3) 150! = -I(0,n2/2,n3-k3), 151! a(1,0,n3-k3) = R(n1/2,0,k3) 152! = R(n1/2,0,n3-k3), 153! a(0,0,n3-k3) = -I(n1/2,0,k3) 154! = I(n1/2,0,n3-k3), 155! a(1,n2/2,n3-k3) = R(n1/2,n2/2,k3) 156! = R(n1/2,n2/2,n3-k3), 157! a(0,n2/2,n3-k3) = -I(n1/2,n2/2,k3) 158! = I(n1/2,n2/2,n3-k3), 159! 0<k3<n3/2, 160! a(0,0,0) = R(0,0,0), 161! a(1,0,0) = R(n1/2,0,0), 162! a(0,0,n3/2) = R(0,0,n3/2), 163! a(1,0,n3/2) = R(n1/2,0,n3/2), 164! a(0,n2/2,0) = R(0,n2/2,0), 165! a(1,n2/2,0) = R(n1/2,n2/2,0), 166! a(0,n2/2,n3/2) = R(0,n2/2,n3/2), 167! a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2) 168! <case2> 169! input data 170! a(2*j1,j2,j3) = R(j1,j2,j3) 171! = R(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), 172! a(2*j1+1,j2,j3) = I(j1,j2,j3) 173! = -I(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), 174! 0<j1<n1/2, 0<=j2<n2, 0<=j3<n3, 175! a(0,j2,j3) = R(0,j2,j3) 176! = R(0,n2-j2,mod(n3-j3,n3)), 177! a(1,j2,j3) = I(0,j2,j3) 178! = -I(0,n2-j2,mod(n3-j3,n3)), 179! a(1,n2-j2,j3) = R(n1/2,j2,mod(n3-j3,n3)) 180! = R(n1/2,n2-j2,j3), 181! a(0,n2-j2,j3) = -I(n1/2,j2,mod(n3-j3,n3)) 182! = I(n1/2,n2-j2,j3), 183! 0<j2<n2/2, 0<=j3<n3, 184! a(0,0,j3) = R(0,0,j3) 185! = R(0,0,n3-j3), 186! a(1,0,j3) = I(0,0,j3) 187! = -I(0,0,n3-j3), 188! a(0,n2/2,j3) = R(0,n2/2,j3) 189! = R(0,n2/2,n3-j3), 190! a(1,n2/2,j3) = I(0,n2/2,j3) 191! = -I(0,n2/2,n3-j3), 192! a(1,0,n3-j3) = R(n1/2,0,j3) 193! = R(n1/2,0,n3-j3), 194! a(0,0,n3-j3) = -I(n1/2,0,j3) 195! = I(n1/2,0,n3-j3), 196! a(1,n2/2,n3-j3) = R(n1/2,n2/2,j3) 197! = R(n1/2,n2/2,n3-j3), 198! a(0,n2/2,n3-j3) = -I(n1/2,n2/2,j3) 199! = I(n1/2,n2/2,n3-j3), 200! 0<j3<n3/2, 201! a(0,0,0) = R(0,0,0), 202! a(1,0,0) = R(n1/2,0,0), 203! a(0,0,n3/2) = R(0,0,n3/2), 204! a(1,0,n3/2) = R(n1/2,0,n3/2), 205! a(0,n2/2,0) = R(0,n2/2,0), 206! a(1,n2/2,0) = R(n1/2,n2/2,0), 207! a(0,n2/2,n3/2) = R(0,n2/2,n3/2), 208! a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2) 209! ---- output ordering ---- 210! call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) 211! call rdft3dsort(n1max, n2max, n1, n2, n3, 1, a) 212! ! stored data is a(0:n1-1,0:n2-1,0:n3+1): 213! ! a(2*k1,k2,k3) = R(k1,k2,k3), 214! ! a(2*k1+1,k2,k3) = I(k1,k2,k3), 215! ! 0<=k1<=n1/2, 0<=k2<n2, 0<=k3<n3. 216! ! the stored data is larger than the input data! 217! ---- input ordering ---- 218! call rdft3dsort(n1max, n2max, n1, n2, n3, -1, a) 219! call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) 220! t(0:*) :work area (real*8) 221! length of t >= max(8*n2, 8*n3) 222! ip(0:*):work area for bit reversal (integer) 223! length of ip >= 2+sqrt(n) 224! (n = max(n1/2, n2, n3)) 225! ip(0),ip(1) are pointers of the cos/sin table. 226! w(0:*) :cos/sin table (real*8) 227! length of w >= max(n1/4, n2/2, n3/2) + n1/4 228! w(),ip() are initialized if ip(0) = 0. 229! [remark] 230! Inverse of 231! call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) 232! is 233! call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) 234! do j3 = 0, n3 - 1 235! do j2 = 0, n2 - 1 236! do j1 = 0, n1 - 1 237! a(j1,j2,j3) = a(j1,j2,j3) * (2.0d0/n1/n2/n3) 238! end do 239! end do 240! end do 241! . 242! 243! 244! -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- 245! [definition] 246! <case1> IDCT (excluding scale) 247! C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 248! a(j1,j2,j3) * 249! cos(pi*j1*(k1+1/2)/n1) * 250! cos(pi*j2*(k2+1/2)/n2) * 251! cos(pi*j3*(k3+1/2)/n3), 252! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 253! <case2> DCT 254! C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 255! a(j1,j2,j3) * 256! cos(pi*(j1+1/2)*k1/n1) * 257! cos(pi*(j2+1/2)*k2/n2) * 258! cos(pi*(j3+1/2)*k3/n3), 259! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 260! [usage] 261! <case1> 262! ip(0) = 0 ! first time only 263! call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) 264! <case2> 265! ip(0) = 0 ! first time only 266! call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) 267! [parameters] 268! n1max :row1 size of the 3D array (integer) 269! n2max :row2 size of the 3D array (integer) 270! n1 :data length (integer) 271! n1 >= 2, n1 = power of 2 272! n2 :data length (integer) 273! n2 >= 2, n2 = power of 2 274! n3 :data length (integer) 275! n3 >= 2, n3 = power of 2 276! a(0:n1-1,0:n2-1,0:n3-1) 277! :input/output data (real*8) 278! output data 279! a(k1,k2,k3) = C(k1,k2,k3), 280! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 281! t(0:*) :work area (real*8) 282! length of t >= max(4*n2, 4*n3) 283! ip(0:*):work area for bit reversal (integer) 284! length of ip >= 2+sqrt(n) 285! (n = max(n1/2, n2/2, n3/2)) 286! ip(0),ip(1) are pointers of the cos/sin table. 287! w(0:*) :cos/sin table (real*8) 288! length of w >= max(n1*3/2, n2*3/2, n3*3/2) 289! w(),ip() are initialized if ip(0) = 0. 290! [remark] 291! Inverse of 292! call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) 293! is 294! do j3 = 0, n3 - 1 295! do j2 = 0, n2 - 1 296! a(0, j2, j3) = a(0, j2, j3) * 0.5d0 297! end do 298! do j1 = 0, n1 - 1 299! a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0 300! end do 301! end do 302! do j2 = 0, n2 - 1 303! do j1 = 0, n1 - 1 304! a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0 305! end do 306! end do 307! call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) 308! do j3 = 0, n3 - 1 309! do j2 = 0, n2 - 1 310! do j1 = 0, n1 - 1 311! a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3) 312! end do 313! end do 314! end do 315! . 316! 317! 318! -------- DST (Discrete Sine Transform) / Inverse of DST -------- 319! [definition] 320! <case1> IDST (excluding scale) 321! S(k1,k2,k3) = sum_j1=1^n1 sum_j2=1^n2 sum_j3=1^n3 322! A(j1,j2,j3) * 323! sin(pi*j1*(k1+1/2)/n1) * 324! sin(pi*j2*(k2+1/2)/n2) * 325! sin(pi*j3*(k3+1/2)/n3), 326! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 327! <case2> DST 328! S(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 329! a(j1,j2,j3) * 330! sin(pi*(j1+1/2)*k1/n1) * 331! sin(pi*(j2+1/2)*k2/n2) * 332! sin(pi*(j3+1/2)*k3/n3), 333! 0<k1<=n1, 0<k2<=n2, 0<k3<=n3 334! [usage] 335! <case1> 336! ip(0) = 0 ! first time only 337! call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) 338! <case2> 339! ip(0) = 0 ! first time only 340! call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) 341! [parameters] 342! n1max :row1 size of the 3D array (integer) 343! n2max :row2 size of the 3D array (integer) 344! n1 :data length (integer) 345! n1 >= 2, n1 = power of 2 346! n2 :data length (integer) 347! n2 >= 2, n2 = power of 2 348! n3 :data length (integer) 349! n3 >= 2, n3 = power of 2 350! a(0:n1-1,0:n2-1,0:n3-1) 351! :input/output data (real*8) 352! <case1> 353! input data 354! a(mod(j1,n1),mod(j2,n2),mod(j3,n3)) = A(j1,j2,j3), 355! 0<j1<=n1, 0<j2<=n2, 0<j3<=n3 356! output data 357! a(k1,k2,k3) = S(k1,k2,k3), 358! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 359! <case2> 360! output data 361! a(mod(k1,n1),mod(k2,n2),mod(k3,n3)) = S(k1,k2,k3), 362! 0<k1<=n1, 0<k2<=n2, 0<k3<=n3 363! t(0:*) :work area (real*8) 364! length of t >= max(4*n2, 4*n3) 365! ip(0:*):work area for bit reversal (integer) 366! length of ip >= 2+sqrt(n) 367! (n = max(n1/2, n2/2, n3/2)) 368! ip(0),ip(1) are pointers of the cos/sin table. 369! w(0:*) :cos/sin table (real*8) 370! length of w >= max(n1*3/2, n2*3/2, n3*3/2) 371! w(),ip() are initialized if ip(0) = 0. 372! [remark] 373! Inverse of 374! call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) 375! is 376! do j3 = 0, n3 - 1 377! do j2 = 0, n2 - 1 378! a(0, j2, j3) = a(0, j2, j3) * 0.5d0 379! end do 380! do j1 = 0, n1 - 1 381! a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0 382! end do 383! end do 384! do j2 = 0, n2 - 1 385! do j1 = 0, n1 - 1 386! a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0 387! end do 388! end do 389! call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) 390! do j3 = 0, n3 - 1 391! do j2 = 0, n2 - 1 392! do j1 = 0, n1 - 1 393! a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3) 394! end do 395! end do 396! end do 397! . 398! 399! 400 subroutine cdft3d(n1max, n2max, n1, n2, n3, isgn, a, 401 & t, ip, w) 402 integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), n 403 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 404 & t(0 : *), w(0 : *) 405 n = 2 * max(n2, n3) 406 n = max(n, n1) 407 if (n .gt. 4 * ip(0)) then 408 call makewt(n / 4, ip, w) 409 end if 410 call xdft3da_sub(n1max, n2max, n1, n2, n3, 0, 411 & isgn, a, t, ip, w) 412 call cdft3db_sub(n1max, n2max, n1, n2, n3, 413 & isgn, a, t, ip, w) 414 end 415! 416 subroutine rdft3d(n1max, n2max, n1, n2, n3, isgn, a, 417 & t, ip, w) 418 integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 419 & n, nw, nc 420 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 421 & t(0 : *), w(0 : *) 422 n = 2 * max(n2, n3) 423 n = max(n, n1) 424 nw = ip(0) 425 if (n .gt. 4 * nw) then 426 nw = n / 4 427 call makewt(nw, ip, w) 428 end if 429 nc = ip(1) 430 if (n1 .gt. 4 * nc) then 431 nc = n1 / 4 432 call makect(nc, ip, w(nw)) 433 end if 434 if (isgn .lt. 0) then 435 call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a) 436 call cdft3db_sub(n1max, n2max, n1, n2, n3, 437 & isgn, a, t, ip, w) 438 call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, 439 & isgn, a, t, ip, w) 440 else 441 call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, 442 & isgn, a, t, ip, w) 443 call cdft3db_sub(n1max, n2max, n1, n2, n3, 444 & isgn, a, t, ip, w) 445 call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a) 446 end if 447 end 448! 449 subroutine rdft3dsort(n1max, n2max, n1, n2, n3, isgn, a) 450 integer n1max, n2max, n1, n2, n3, isgn, n2h, n3h, j, k 451 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), x, y 452 n2h = n2 / 2 453 n3h = n3 / 2 454 if (isgn .lt. 0) then 455 do k = 0, n3 - 1 456 do j = n2h + 1, n2 - 1 457 a(0, j, k) = a(n1 + 1, j, k) 458 a(1, j, k) = a(n1, j, k) 459 end do 460 end do 461 do k = n3h + 1, n3 - 1 462 a(0, 0, k) = a(n1 + 1, 0, k) 463 a(1, 0, k) = a(n1, 0, k) 464 a(0, n2h, k) = a(n1 + 1, n2h, k) 465 a(1, n2h, k) = a(n1, n2h, k) 466 end do 467 a(1, 0, 0) = a(n1, 0, 0) 468 a(1, n2h, 0) = a(n1, n2h, 0) 469 a(1, 0, n3h) = a(n1, 0, n3h) 470 a(1, n2h, n3h) = a(n1, n2h, n3h) 471 else 472 do j = n2h + 1, n2 - 1 473 y = a(0, j, 0) 474 x = a(1, j, 0) 475 a(n1, j, 0) = x 476 a(n1 + 1, j, 0) = y 477 a(n1, n2 - j, 0) = x 478 a(n1 + 1, n2 - j, 0) = -y 479 a(0, j, 0) = a(0, n2 - j, 0) 480 a(1, j, 0) = -a(1, n2 - j, 0) 481 end do 482 do k = 1, n3 - 1 483 do j = n2h + 1, n2 - 1 484 y = a(0, j, k) 485 x = a(1, j, k) 486 a(n1, j, k) = x 487 a(n1 + 1, j, k) = y 488 a(n1, n2 - j, n3 - k) = x 489 a(n1 + 1, n2 - j, n3 - k) = -y 490 a(0, j, k) = a(0, n2 - j, n3 - k) 491 a(1, j, k) = -a(1, n2 - j, n3 - k) 492 end do 493 end do 494 do k = n3h + 1, n3 - 1 495 y = a(0, 0, k) 496 x = a(1, 0, k) 497 a(n1, 0, k) = x 498 a(n1 + 1, 0, k) = y 499 a(n1, 0, n3 - k) = x 500 a(n1 + 1, 0, n3 - k) = -y 501 a(0, 0, k) = a(0, 0, n3 - k) 502 a(1, 0, k) = -a(1, 0, n3 - k) 503 y = a(0, n2h, k) 504 x = a(1, n2h, k) 505 a(n1, n2h, k) = x 506 a(n1 + 1, n2h, k) = y 507 a(n1, n2h, n3 - k) = x 508 a(n1 + 1, n2h, n3 - k) = -y 509 a(0, n2h, k) = a(0, n2h, n3 - k) 510 a(1, n2h, k) = -a(1, n2h, n3 - k) 511 end do 512 a(n1, 0, 0) = a(1, 0, 0) 513 a(n1 + 1, 0, 0) = 0 514 a(1, 0, 0) = 0 515 a(n1, n2h, 0) = a(1, n2h, 0) 516 a(n1 + 1, n2h, 0) = 0 517 a(1, n2h, 0) = 0 518 a(n1, 0, n3h) = a(1, 0, n3h) 519 a(n1 + 1, 0, n3h) = 0 520 a(1, 0, n3h) = 0 521 a(n1, n2h, n3h) = a(1, n2h, n3h) 522 a(n1 + 1, n2h, n3h) = 0 523 a(1, n2h, n3h) = 0 524 end if 525 end 526! 527 subroutine ddct3d(n1max, n2max, n1, n2, n3, isgn, a, 528 & t, ip, w) 529 integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 530 & n, nw, nc 531 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 532 & t(0 : *), w(0 : *) 533 n = max(n2, n3) 534 n = max(n, n1) 535 nw = ip(0) 536 if (n .gt. 4 * nw) then 537 nw = n / 4 538 call makewt(nw, ip, w) 539 end if 540 nc = ip(1) 541 if (n .gt. nc) then 542 nc = n 543 call makect(nc, ip, w(nw)) 544 end if 545 call ddxt3da_sub(n1max, n2max, n1, n2, n3, 0, 546 & isgn, a, t, ip, w) 547 call ddxt3db_sub(n1max, n2max, n1, n2, n3, 0, 548 & isgn, a, t, ip, w) 549 end 550! 551 subroutine ddst3d(n1max, n2max, n1, n2, n3, isgn, a, 552 & t, ip, w) 553 integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 554 & n, nw, nc 555 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 556 & t(0 : *), w(0 : *) 557 n = max(n2, n3) 558 n = max(n, n1) 559 nw = ip(0) 560 if (n .gt. 4 * nw) then 561 nw = n / 4 562 call makewt(nw, ip, w) 563 end if 564 nc = ip(1) 565 if (n .gt. nc) then 566 nc = n 567 call makect(nc, ip, w(nw)) 568 end if 569 call ddxt3da_sub(n1max, n2max, n1, n2, n3, 1, 570 & isgn, a, t, ip, w) 571 call ddxt3db_sub(n1max, n2max, n1, n2, n3, 1, 572 & isgn, a, t, ip, w) 573 end 574! 575! -------- child routines -------- 576! 577 subroutine xdft3da_sub(n1max, n2max, n1, n2, n3, icr, 578 & isgn, a, t, ip, w) 579 integer n1max, n2max, n1, n2, n3, icr, isgn, 580 & ip(0 : *), i, j, k 581 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 582 & t(0 : *), w(0 : *) 583 do k = 0, n3 - 1 584 if (icr .eq. 0) then 585 do j = 0, n2 - 1 586 call cdft(n1, isgn, a(0, j, k), ip, w) 587 end do 588 else if (isgn .ge. 0) then 589 do j = 0, n2 - 1 590 call rdft(n1, isgn, a(0, j, k), ip, w) 591 end do 592 end if 593 if (n1 .gt. 4) then 594 do i = 0, n1 - 8, 8 595 do j = 0, n2 - 1 596 t(2 * j) = a(i, j, k) 597 t(2 * j + 1) = a(i + 1, j, k) 598 t(2 * n2 + 2 * j) = a(i + 2, j, k) 599 t(2 * n2 + 2 * j + 1) = a(i + 3, j, k) 600 t(4 * n2 + 2 * j) = a(i + 4, j, k) 601 t(4 * n2 + 2 * j + 1) = a(i + 5, j, k) 602 t(6 * n2 + 2 * j) = a(i + 6, j, k) 603 t(6 * n2 + 2 * j + 1) = a(i + 7, j, k) 604 end do 605 call cdft(2 * n2, isgn, t, ip, w) 606 call cdft(2 * n2, isgn, t(2 * n2), ip, w) 607 call cdft(2 * n2, isgn, t(4 * n2), ip, w) 608 call cdft(2 * n2, isgn, t(6 * n2), ip, w) 609 do j = 0, n2 - 1 610 a(i, j, k) = t(2 * j) 611 a(i + 1, j, k) = t(2 * j + 1) 612 a(i + 2, j, k) = t(2 * n2 + 2 * j) 613 a(i + 3, j, k) = t(2 * n2 + 2 * j + 1) 614 a(i + 4, j, k) = t(4 * n2 + 2 * j) 615 a(i + 5, j, k) = t(4 * n2 + 2 * j + 1) 616 a(i + 6, j, k) = t(6 * n2 + 2 * j) 617 a(i + 7, j, k) = t(6 * n2 + 2 * j + 1) 618 end do 619 end do 620 else if (n1 .eq. 4) then 621 do j = 0, n2 - 1 622 t(2 * j) = a(0, j, k) 623 t(2 * j + 1) = a(1, j, k) 624 t(2 * n2 + 2 * j) = a(2, j, k) 625 t(2 * n2 + 2 * j + 1) = a(3, j, k) 626 end do 627 call cdft(2 * n2, isgn, t, ip, w) 628 call cdft(2 * n2, isgn, t(2 * n2), ip, w) 629 do j = 0, n2 - 1 630 a(0, j, k) = t(2 * j) 631 a(1, j, k) = t(2 * j + 1) 632 a(2, j, k) = t(2 * n2 + 2 * j) 633 a(3, j, k) = t(2 * n2 + 2 * j + 1) 634 end do 635 else if (n1 .eq. 2) then 636 do j = 0, n2 - 1 637 t(2 * j) = a(0, j, k) 638 t(2 * j + 1) = a(1, j, k) 639 end do 640 call cdft(2 * n2, isgn, t, ip, w) 641 do j = 0, n2 - 1 642 a(0, j, k) = t(2 * j) 643 a(1, j, k) = t(2 * j + 1) 644 end do 645 end if 646 if (icr .ne. 0 .and. isgn .lt. 0) then 647 do j = 0, n2 - 1 648 call rdft(n1, isgn, a(0, j, k), ip, w) 649 end do 650 end if 651 end do 652 end 653! 654 subroutine cdft3db_sub(n1max, n2max, n1, n2, n3, 655 & isgn, a, t, ip, w) 656 integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 657 & i, j, k 658 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 659 & t(0 : *), w(0 : *) 660 if (n1 .gt. 4) then 661 do j = 0, n2 - 1 662 do i = 0, n1 - 8, 8 663 do k = 0, n3 - 1 664 t(2 * k) = a(i, j, k) 665 t(2 * k + 1) = a(i + 1, j, k) 666 t(2 * n3 + 2 * k) = a(i + 2, j, k) 667 t(2 * n3 + 2 * k + 1) = a(i + 3, j, k) 668 t(4 * n3 + 2 * k) = a(i + 4, j, k) 669 t(4 * n3 + 2 * k + 1) = a(i + 5, j, k) 670 t(6 * n3 + 2 * k) = a(i + 6, j, k) 671 t(6 * n3 + 2 * k + 1) = a(i + 7, j, k) 672 end do 673 call cdft(2 * n3, isgn, t, ip, w) 674 call cdft(2 * n3, isgn, t(2 * n3), ip, w) 675 call cdft(2 * n3, isgn, t(4 * n3), ip, w) 676 call cdft(2 * n3, isgn, t(6 * n3), ip, w) 677 do k = 0, n3 - 1 678 a(i, j, k) = t(2 * k) 679 a(i + 1, j, k) = t(2 * k + 1) 680 a(i + 2, j, k) = t(2 * n3 + 2 * k) 681 a(i + 3, j, k) = t(2 * n3 + 2 * k + 1) 682 a(i + 4, j, k) = t(4 * n3 + 2 * k) 683 a(i + 5, j, k) = t(4 * n3 + 2 * k + 1) 684 a(i + 6, j, k) = t(6 * n3 + 2 * k) 685 a(i + 7, j, k) = t(6 * n3 + 2 * k + 1) 686 end do 687 end do 688 end do 689 else if (n1 .eq. 4) then 690 do j = 0, n2 - 1 691 do k = 0, n3 - 1 692 t(2 * k) = a(0, j, k) 693 t(2 * k + 1) = a(1, j, k) 694 t(2 * n3 + 2 * k) = a(2, j, k) 695 t(2 * n3 + 2 * k + 1) = a(3, j, k) 696 end do 697 call cdft(2 * n3, isgn, t, ip, w) 698 call cdft(2 * n3, isgn, t(2 * n3), ip, w) 699 do k = 0, n3 - 1 700 a(0, j, k) = t(2 * k) 701 a(1, j, k) = t(2 * k + 1) 702 a(2, j, k) = t(2 * n3 + 2 * k) 703 a(3, j, k) = t(2 * n3 + 2 * k + 1) 704 end do 705 end do 706 else if (n1 .eq. 2) then 707 do j = 0, n2 - 1 708 do k = 0, n3 - 1 709 t(2 * k) = a(0, j, k) 710 t(2 * k + 1) = a(1, j, k) 711 end do 712 call cdft(2 * n3, isgn, t, ip, w) 713 do k = 0, n3 - 1 714 a(0, j, k) = t(2 * k) 715 a(1, j, k) = t(2 * k + 1) 716 end do 717 end do 718 end if 719 end 720! 721 subroutine rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a) 722 integer n1max, n2max, n1, n2, n3, isgn, 723 & n2h, n3h, i, j, k, l 724 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), xi 725 n2h = n2 / 2 726 n3h = n3 / 2 727 if (isgn .lt. 0) then 728 do k = 1, n3h - 1 729 l = n3 - k 730 xi = a(0, 0, k) - a(0, 0, l) 731 a(0, 0, k) = a(0, 0, k) + a(0, 0, l) 732 a(0, 0, l) = xi 733 xi = a(1, 0, l) - a(1, 0, k) 734 a(1, 0, k) = a(1, 0, k) + a(1, 0, l) 735 a(1, 0, l) = xi 736 xi = a(0, n2h, k) - a(0, n2h, l) 737 a(0, n2h, k) = a(0, n2h, k) + a(0, n2h, l) 738 a(0, n2h, l) = xi 739 xi = a(1, n2h, l) - a(1, n2h, k) 740 a(1, n2h, k) = a(1, n2h, k) + a(1, n2h, l) 741 a(1, n2h, l) = xi 742 do i = 1, n2h - 1 743 j = n2 - i 744 xi = a(0, i, k) - a(0, j, l) 745 a(0, i, k) = a(0, i, k) + a(0, j, l) 746 a(0, j, l) = xi 747 xi = a(1, j, l) - a(1, i, k) 748 a(1, i, k) = a(1, i, k) + a(1, j, l) 749 a(1, j, l) = xi 750 xi = a(0, i, l) - a(0, j, k) 751 a(0, i, l) = a(0, i, l) + a(0, j, k) 752 a(0, j, k) = xi 753 xi = a(1, j, k) - a(1, i, l) 754 a(1, i, l) = a(1, i, l) + a(1, j, k) 755 a(1, j, k) = xi 756 end do 757 end do 758 do i = 1, n2h - 1 759 j = n2 - i 760 xi = a(0, i, 0) - a(0, j, 0) 761 a(0, i, 0) = a(0, i, 0) + a(0, j, 0) 762 a(0, j, 0) = xi 763 xi = a(1, j, 0) - a(1, i, 0) 764 a(1, i, 0) = a(1, i, 0) + a(1, j, 0) 765 a(1, j, 0) = xi 766 xi = a(0, i, n3h) - a(0, j, n3h) 767 a(0, i, n3h) = a(0, i, n3h) + a(0, j, n3h) 768 a(0, j, n3h) = xi 769 xi = a(1, j, n3h) - a(1, i, n3h) 770 a(1, i, n3h) = a(1, i, n3h) + a(1, j, n3h) 771 a(1, j, n3h) = xi 772 end do 773 else 774 do k = 1, n3h - 1 775 l = n3 - k 776 a(0, 0, l) = 0.5d0 * (a(0, 0, k) - a(0, 0, l)) 777 a(0, 0, k) = a(0, 0, k) - a(0, 0, l) 778 a(1, 0, l) = 0.5d0 * (a(1, 0, k) + a(1, 0, l)) 779 a(1, 0, k) = a(1, 0, k) - a(1, 0, l) 780 a(0, n2h, l) = 0.5d0 * (a(0, n2h, k) - a(0, n2h, l)) 781 a(0, n2h, k) = a(0, n2h, k) - a(0, n2h, l) 782 a(1, n2h, l) = 0.5d0 * (a(1, n2h, k) + a(1, n2h, l)) 783 a(1, n2h, k) = a(1, n2h, k) - a(1, n2h, l) 784 do i = 1, n2h - 1 785 j = n2 - i 786 a(0, j, l) = 0.5d0 * (a(0, i, k) - a(0, j, l)) 787 a(0, i, k) = a(0, i, k) - a(0, j, l) 788 a(1, j, l) = 0.5d0 * (a(1, i, k) + a(1, j, l)) 789 a(1, i, k) = a(1, i, k) - a(1, j, l) 790 a(0, j, k) = 0.5d0 * (a(0, i, l) - a(0, j, k)) 791 a(0, i, l) = a(0, i, l) - a(0, j, k) 792 a(1, j, k) = 0.5d0 * (a(1, i, l) + a(1, j, k)) 793 a(1, i, l) = a(1, i, l) - a(1, j, k) 794 end do 795 end do 796 do i = 1, n2h - 1 797 j = n2 - i 798 a(0, j, 0) = 0.5d0 * (a(0, i, 0) - a(0, j, 0)) 799 a(0, i, 0) = a(0, i, 0) - a(0, j, 0) 800 a(1, j, 0) = 0.5d0 * (a(1, i, 0) + a(1, j, 0)) 801 a(1, i, 0) = a(1, i, 0) - a(1, j, 0) 802 a(0, j, n3h) = 0.5d0 * (a(0, i, n3h) - a(0, j, n3h)) 803 a(0, i, n3h) = a(0, i, n3h) - a(0, j, n3h) 804 a(1, j, n3h) = 0.5d0 * (a(1, i, n3h) + a(1, j, n3h)) 805 a(1, i, n3h) = a(1, i, n3h) - a(1, j, n3h) 806 end do 807 end if 808 end 809! 810 subroutine ddxt3da_sub(n1max, n2max, n1, n2, n3, ics, 811 & isgn, a, t, ip, w) 812 integer n1max, n2max, n1, n2, n3, ics, isgn, 813 & ip(0 : *), i, j, k 814 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 815 & t(0 : *), w(0 : *) 816 do k = 0, n3 - 1 817 if (ics .eq. 0) then 818 do j = 0, n2 - 1 819 call ddct(n1, isgn, a(0, j, k), ip, w) 820 end do 821 else 822 do j = 0, n2 - 1 823 call ddst(n1, isgn, a(0, j, k), ip, w) 824 end do 825 end if 826 if (n1 .gt. 2) then 827 do i = 0, n1 - 4, 4 828 do j = 0, n2 - 1 829 t(j) = a(i, j, k) 830 t(n2 + j) = a(i + 1, j, k) 831 t(2 * n2 + j) = a(i + 2, j, k) 832 t(3 * n2 + j) = a(i + 3, j, k) 833 end do 834 if (ics .eq. 0) then 835 call ddct(n2, isgn, t, ip, w) 836 call ddct(n2, isgn, t(n2), ip, w) 837 call ddct(n2, isgn, t(2 * n2), ip, w) 838 call ddct(n2, isgn, t(3 * n2), ip, w) 839 else 840 call ddst(n2, isgn, t, ip, w) 841 call ddst(n2, isgn, t(n2), ip, w) 842 call ddst(n2, isgn, t(2 * n2), ip, w) 843 call ddst(n2, isgn, t(3 * n2), ip, w) 844 end if 845 do j = 0, n2 - 1 846 a(i, j, k) = t(j) 847 a(i + 1, j, k) = t(n2 + j) 848 a(i + 2, j, k) = t(2 * n2 + j) 849 a(i + 3, j, k) = t(3 * n2 + j) 850 end do 851 end do 852 else if (n1 .eq. 2) then 853 do j = 0, n2 - 1 854 t(j) = a(0, j, k) 855 t(n2 + j) = a(1, j, k) 856 end do 857 if (ics .eq. 0) then 858 call ddct(n2, isgn, t, ip, w) 859 call ddct(n2, isgn, t(n2), ip, w) 860 else 861 call ddst(n2, isgn, t, ip, w) 862 call ddst(n2, isgn, t(n2), ip, w) 863 end if 864 do j = 0, n2 - 1 865 a(0, j, k) = t(j) 866 a(1, j, k) = t(n2 + j) 867 end do 868 end if 869 end do 870 end 871! 872 subroutine ddxt3db_sub(n1max, n2max, n1, n2, n3, ics, 873 & isgn, a, t, ip, w) 874 integer n1max, n2max, n1, n2, n3, ics, isgn, 875 & ip(0 : *), i, j, k 876 real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 877 & t(0 : *), w(0 : *) 878 if (n1 .gt. 2) then 879 do j = 0, n2 - 1 880 do i = 0, n1 - 4, 4 881 do k = 0, n3 - 1 882 t(k) = a(i, j, k) 883 t(n3 + k) = a(i + 1, j, k) 884 t(2 * n3 + k) = a(i + 2, j, k) 885 t(3 * n3 + k) = a(i + 3, j, k) 886 end do 887 if (ics .eq. 0) then 888 call ddct(n3, isgn, t, ip, w) 889 call ddct(n3, isgn, t(n3), ip, w) 890 call ddct(n3, isgn, t(2 * n3), ip, w) 891 call ddct(n3, isgn, t(3 * n3), ip, w) 892 else 893 call ddst(n3, isgn, t, ip, w) 894 call ddst(n3, isgn, t(n3), ip, w) 895 call ddst(n3, isgn, t(2 * n3), ip, w) 896 call ddst(n3, isgn, t(3 * n3), ip, w) 897 end if 898 do k = 0, n3 - 1 899 a(i, j, k) = t(k) 900 a(i + 1, j, k) = t(n3 + k) 901 a(i + 2, j, k) = t(2 * n3 + k) 902 a(i + 3, j, k) = t(3 * n3 + k) 903 end do 904 end do 905 end do 906 else if (n1 .eq. 2) then 907 do j = 0, n2 - 1 908 do k = 0, n3 - 1 909 t(k) = a(0, j, k) 910 t(n3 + k) = a(1, j, k) 911 end do 912 if (ics .eq. 0) then 913 call ddct(n3, isgn, t, ip, w) 914 call ddct(n3, isgn, t(n3), ip, w) 915 else 916 call ddst(n3, isgn, t, ip, w) 917 call ddst(n3, isgn, t(n3), ip, w) 918 end if 919 do k = 0, n3 - 1 920 a(0, j, k) = t(k) 921 a(1, j, k) = t(n3 + k) 922 end do 923 end do 924 end if 925 end 926! 927