| ! Fast Fourier/Cosine/Sine Transform | |
| ! dimension :three | |
| ! data length :power of 2 | |
| ! decimation :frequency | |
| ! radix :split-radix, row-column | |
| ! data :inplace | |
| ! table :use | |
| ! subroutines | |
| ! cdft3d: Complex Discrete Fourier Transform | |
| ! rdft3d: Real Discrete Fourier Transform | |
| ! ddct3d: Discrete Cosine Transform | |
| ! ddst3d: Discrete Sine Transform | |
| ! necessary package | |
| ! fftsg.f : 1D-FFT package | |
| ! | |
| ! | |
| ! -------- Complex DFT (Discrete Fourier Transform) -------- | |
| ! [definition] | |
| ! <case1> | |
| ! X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! x(j1,j2,j3) * | |
| ! exp(2*pi*i*j1*k1/n1) * | |
| ! exp(2*pi*i*j2*k2/n2) * | |
| ! exp(2*pi*i*j3*k3/n3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! <case2> | |
| ! X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! x(j1,j2,j3) * | |
| ! exp(-2*pi*i*j1*k1/n1) * | |
| ! exp(-2*pi*i*j2*k2/n2) * | |
| ! exp(-2*pi*i*j3*k3/n3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! (notes: sum_j=0^n-1 is a summation from j=0 to n-1) | |
| ! [usage] | |
| ! <case1> | |
| ! ip(0) = 0 ! first time only | |
| ! call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w) | |
| ! <case2> | |
| ! ip(0) = 0 ! first time only | |
| ! call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w) | |
| ! [parameters] | |
| ! n1max :row1 size of the 3D array (integer) | |
| ! n2max :row2 size of the 3D array (integer) | |
| ! 2*n1 :data length (integer) | |
| ! n1 >= 1, n1 = power of 2 | |
| ! n2 :data length (integer) | |
| ! n2 >= 1, n2 = power of 2 | |
| ! n3 :data length (integer) | |
| ! n3 >= 1, n3 = power of 2 | |
| ! a(0:2*n1-1,0:n2-1,0:n3-1) | |
| ! :input/output data (real*8) | |
| ! input data | |
| ! a(2*j1,j2,j3) = Re(x(j1,j2,j3)), | |
| ! a(2*j1+1,j2,j3) = Im(x(j1,j2,j3)), | |
| ! 0<=j1<n1, 0<=j2<n2, 0<=j3<n3 | |
| ! output data | |
| ! a(2*k1,k2,k3) = Re(X(k1,k2,k3)), | |
| ! a(2*k1+1,k2,k3) = Im(X(k1,k2,k3)), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! t(0:*) :work area (real*8) | |
| ! length of t >= max(8*n2, 8*n3) | |
| ! ip(0:*):work area for bit reversal (integer) | |
| ! length of ip >= 2+sqrt(n) | |
| ! (n = max(n1, n2, n3)) | |
| ! ip(0),ip(1) are pointers of the cos/sin table. | |
| ! w(0:*) :cos/sin table (real*8) | |
| ! length of w >= max(n1/2, n2/2, n3/2) | |
| ! w(),ip() are initialized if ip(0) = 0. | |
| ! [remark] | |
| ! Inverse of | |
| ! call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w) | |
| ! is | |
| ! call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w) | |
| ! do j3 = 0, n3 - 1 | |
| ! do j2 = 0, n2 - 1 | |
| ! do j1 = 0, 2 * n1 - 1 | |
| ! a(j1,j2,j3) = a(j1,j2,j3) * (1.0d0/n1/n2/n3) | |
| ! end do | |
| ! end do | |
| ! end do | |
| ! . | |
| ! | |
| ! | |
| ! -------- Real DFT / Inverse of Real DFT -------- | |
| ! [definition] | |
| ! <case1> RDFT | |
| ! R(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! a(j1,j2,j3) * | |
| ! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + | |
| ! 2*pi*j3*k3/n3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! I(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! a(j1,j2,j3) * | |
| ! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + | |
| ! 2*pi*j3*k3/n3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! <case2> IRDFT (excluding scale) | |
| ! a(k1,k2,k3) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! (R(j1,j2,j3) * | |
| ! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + | |
| ! 2*pi*j3*k3/n3) + | |
| ! I(j1,j2,j3) * | |
| ! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + | |
| ! 2*pi*j3*k3/n3)), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! (notes: R(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = R(k1,k2,k3), | |
| ! I(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = -I(k1,k2,k3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3) | |
| ! [usage] | |
| ! <case1> | |
| ! ip(0) = 0 ! first time only | |
| ! call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) | |
| ! <case2> | |
| ! ip(0) = 0 ! first time only | |
| ! call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) | |
| ! [parameters] | |
| ! n1max :row1 size of the 3D array (integer) | |
| ! n2max :row2 size of the 3D array (integer) | |
| ! n1 :data length (integer) | |
| ! n1 >= 2, n1 = power of 2 | |
| ! n2 :data length (integer) | |
| ! n2 >= 2, n2 = power of 2 | |
| ! n3 :data length (integer) | |
| ! n3 >= 2, n3 = power of 2 | |
| ! a(0:n1-1,0:n2-1,0:n3-1) | |
| ! :input/output data (real*8) | |
| ! <case1> | |
| ! output data | |
| ! a(2*k1,k2,k3) = R(k1,k2,k3) | |
| ! = R(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), | |
| ! a(2*k1+1,k2,k3) = I(k1,k2,k3) | |
| ! = -I(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), | |
| ! 0<k1<n1/2, 0<=k2<n2, 0<=k3<n3, | |
| ! a(0,k2,k3) = R(0,k2,k3) | |
| ! = R(0,n2-k2,mod(n3-k3,n3)), | |
| ! a(1,k2,k3) = I(0,k2,k3) | |
| ! = -I(0,n2-k2,mod(n3-k3,n3)), | |
| ! a(1,n2-k2,k3) = R(n1/2,k2,mod(n3-k3,n3)) | |
| ! = R(n1/2,n2-k2,k3), | |
| ! a(0,n2-k2,k3) = -I(n1/2,k2,mod(n3-k3,n3)) | |
| ! = I(n1/2,n2-k2,k3), | |
| ! 0<k2<n2/2, 0<=k3<n3, | |
| ! a(0,0,k3) = R(0,0,k3) | |
| ! = R(0,0,n3-k3), | |
| ! a(1,0,k3) = I(0,0,k3) | |
| ! = -I(0,0,n3-k3), | |
| ! a(0,n2/2,k3) = R(0,n2/2,k3) | |
| ! = R(0,n2/2,n3-k3), | |
| ! a(1,n2/2,k3) = I(0,n2/2,k3) | |
| ! = -I(0,n2/2,n3-k3), | |
| ! a(1,0,n3-k3) = R(n1/2,0,k3) | |
| ! = R(n1/2,0,n3-k3), | |
| ! a(0,0,n3-k3) = -I(n1/2,0,k3) | |
| ! = I(n1/2,0,n3-k3), | |
| ! a(1,n2/2,n3-k3) = R(n1/2,n2/2,k3) | |
| ! = R(n1/2,n2/2,n3-k3), | |
| ! a(0,n2/2,n3-k3) = -I(n1/2,n2/2,k3) | |
| ! = I(n1/2,n2/2,n3-k3), | |
| ! 0<k3<n3/2, | |
| ! a(0,0,0) = R(0,0,0), | |
| ! a(1,0,0) = R(n1/2,0,0), | |
| ! a(0,0,n3/2) = R(0,0,n3/2), | |
| ! a(1,0,n3/2) = R(n1/2,0,n3/2), | |
| ! a(0,n2/2,0) = R(0,n2/2,0), | |
| ! a(1,n2/2,0) = R(n1/2,n2/2,0), | |
| ! a(0,n2/2,n3/2) = R(0,n2/2,n3/2), | |
| ! a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2) | |
| ! <case2> | |
| ! input data | |
| ! a(2*j1,j2,j3) = R(j1,j2,j3) | |
| ! = R(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), | |
| ! a(2*j1+1,j2,j3) = I(j1,j2,j3) | |
| ! = -I(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), | |
| ! 0<j1<n1/2, 0<=j2<n2, 0<=j3<n3, | |
| ! a(0,j2,j3) = R(0,j2,j3) | |
| ! = R(0,n2-j2,mod(n3-j3,n3)), | |
| ! a(1,j2,j3) = I(0,j2,j3) | |
| ! = -I(0,n2-j2,mod(n3-j3,n3)), | |
| ! a(1,n2-j2,j3) = R(n1/2,j2,mod(n3-j3,n3)) | |
| ! = R(n1/2,n2-j2,j3), | |
| ! a(0,n2-j2,j3) = -I(n1/2,j2,mod(n3-j3,n3)) | |
| ! = I(n1/2,n2-j2,j3), | |
| ! 0<j2<n2/2, 0<=j3<n3, | |
| ! a(0,0,j3) = R(0,0,j3) | |
| ! = R(0,0,n3-j3), | |
| ! a(1,0,j3) = I(0,0,j3) | |
| ! = -I(0,0,n3-j3), | |
| ! a(0,n2/2,j3) = R(0,n2/2,j3) | |
| ! = R(0,n2/2,n3-j3), | |
| ! a(1,n2/2,j3) = I(0,n2/2,j3) | |
| ! = -I(0,n2/2,n3-j3), | |
| ! a(1,0,n3-j3) = R(n1/2,0,j3) | |
| ! = R(n1/2,0,n3-j3), | |
| ! a(0,0,n3-j3) = -I(n1/2,0,j3) | |
| ! = I(n1/2,0,n3-j3), | |
| ! a(1,n2/2,n3-j3) = R(n1/2,n2/2,j3) | |
| ! = R(n1/2,n2/2,n3-j3), | |
| ! a(0,n2/2,n3-j3) = -I(n1/2,n2/2,j3) | |
| ! = I(n1/2,n2/2,n3-j3), | |
| ! 0<j3<n3/2, | |
| ! a(0,0,0) = R(0,0,0), | |
| ! a(1,0,0) = R(n1/2,0,0), | |
| ! a(0,0,n3/2) = R(0,0,n3/2), | |
| ! a(1,0,n3/2) = R(n1/2,0,n3/2), | |
| ! a(0,n2/2,0) = R(0,n2/2,0), | |
| ! a(1,n2/2,0) = R(n1/2,n2/2,0), | |
| ! a(0,n2/2,n3/2) = R(0,n2/2,n3/2), | |
| ! a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2) | |
| ! ---- output ordering ---- | |
| ! call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) | |
| ! call rdft3dsort(n1max, n2max, n1, n2, n3, 1, a) | |
| ! ! stored data is a(0:n1-1,0:n2-1,0:n3+1): | |
| ! ! a(2*k1,k2,k3) = R(k1,k2,k3), | |
| ! ! a(2*k1+1,k2,k3) = I(k1,k2,k3), | |
| ! ! 0<=k1<=n1/2, 0<=k2<n2, 0<=k3<n3. | |
| ! ! the stored data is larger than the input data! | |
| ! ---- input ordering ---- | |
| ! call rdft3dsort(n1max, n2max, n1, n2, n3, -1, a) | |
| ! call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) | |
| ! t(0:*) :work area (real*8) | |
| ! length of t >= max(8*n2, 8*n3) | |
| ! ip(0:*):work area for bit reversal (integer) | |
| ! length of ip >= 2+sqrt(n) | |
| ! (n = max(n1/2, n2, n3)) | |
| ! ip(0),ip(1) are pointers of the cos/sin table. | |
| ! w(0:*) :cos/sin table (real*8) | |
| ! length of w >= max(n1/4, n2/2, n3/2) + n1/4 | |
| ! w(),ip() are initialized if ip(0) = 0. | |
| ! [remark] | |
| ! Inverse of | |
| ! call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) | |
| ! is | |
| ! call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) | |
| ! do j3 = 0, n3 - 1 | |
| ! do j2 = 0, n2 - 1 | |
| ! do j1 = 0, n1 - 1 | |
| ! a(j1,j2,j3) = a(j1,j2,j3) * (2.0d0/n1/n2/n3) | |
| ! end do | |
| ! end do | |
| ! end do | |
| ! . | |
| ! | |
| ! | |
| ! -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- | |
| ! [definition] | |
| ! <case1> IDCT (excluding scale) | |
| ! C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! a(j1,j2,j3) * | |
| ! cos(pi*j1*(k1+1/2)/n1) * | |
| ! cos(pi*j2*(k2+1/2)/n2) * | |
| ! cos(pi*j3*(k3+1/2)/n3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! <case2> DCT | |
| ! C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! a(j1,j2,j3) * | |
| ! cos(pi*(j1+1/2)*k1/n1) * | |
| ! cos(pi*(j2+1/2)*k2/n2) * | |
| ! cos(pi*(j3+1/2)*k3/n3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! [usage] | |
| ! <case1> | |
| ! ip(0) = 0 ! first time only | |
| ! call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) | |
| ! <case2> | |
| ! ip(0) = 0 ! first time only | |
| ! call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) | |
| ! [parameters] | |
| ! n1max :row1 size of the 3D array (integer) | |
| ! n2max :row2 size of the 3D array (integer) | |
| ! n1 :data length (integer) | |
| ! n1 >= 2, n1 = power of 2 | |
| ! n2 :data length (integer) | |
| ! n2 >= 2, n2 = power of 2 | |
| ! n3 :data length (integer) | |
| ! n3 >= 2, n3 = power of 2 | |
| ! a(0:n1-1,0:n2-1,0:n3-1) | |
| ! :input/output data (real*8) | |
| ! output data | |
| ! a(k1,k2,k3) = C(k1,k2,k3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! t(0:*) :work area (real*8) | |
| ! length of t >= max(4*n2, 4*n3) | |
| ! ip(0:*):work area for bit reversal (integer) | |
| ! length of ip >= 2+sqrt(n) | |
| ! (n = max(n1/2, n2/2, n3/2)) | |
| ! ip(0),ip(1) are pointers of the cos/sin table. | |
| ! w(0:*) :cos/sin table (real*8) | |
| ! length of w >= max(n1*3/2, n2*3/2, n3*3/2) | |
| ! w(),ip() are initialized if ip(0) = 0. | |
| ! [remark] | |
| ! Inverse of | |
| ! call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) | |
| ! is | |
| ! do j3 = 0, n3 - 1 | |
| ! do j2 = 0, n2 - 1 | |
| ! a(0, j2, j3) = a(0, j2, j3) * 0.5d0 | |
| ! end do | |
| ! do j1 = 0, n1 - 1 | |
| ! a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0 | |
| ! end do | |
| ! end do | |
| ! do j2 = 0, n2 - 1 | |
| ! do j1 = 0, n1 - 1 | |
| ! a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0 | |
| ! end do | |
| ! end do | |
| ! call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) | |
| ! do j3 = 0, n3 - 1 | |
| ! do j2 = 0, n2 - 1 | |
| ! do j1 = 0, n1 - 1 | |
| ! a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3) | |
| ! end do | |
| ! end do | |
| ! end do | |
| ! . | |
| ! | |
| ! | |
| ! -------- DST (Discrete Sine Transform) / Inverse of DST -------- | |
| ! [definition] | |
| ! <case1> IDST (excluding scale) | |
| ! S(k1,k2,k3) = sum_j1=1^n1 sum_j2=1^n2 sum_j3=1^n3 | |
| ! A(j1,j2,j3) * | |
| ! sin(pi*j1*(k1+1/2)/n1) * | |
| ! sin(pi*j2*(k2+1/2)/n2) * | |
| ! sin(pi*j3*(k3+1/2)/n3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! <case2> DST | |
| ! S(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1 | |
| ! a(j1,j2,j3) * | |
| ! sin(pi*(j1+1/2)*k1/n1) * | |
| ! sin(pi*(j2+1/2)*k2/n2) * | |
| ! sin(pi*(j3+1/2)*k3/n3), | |
| ! 0<k1<=n1, 0<k2<=n2, 0<k3<=n3 | |
| ! [usage] | |
| ! <case1> | |
| ! ip(0) = 0 ! first time only | |
| ! call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) | |
| ! <case2> | |
| ! ip(0) = 0 ! first time only | |
| ! call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) | |
| ! [parameters] | |
| ! n1max :row1 size of the 3D array (integer) | |
| ! n2max :row2 size of the 3D array (integer) | |
| ! n1 :data length (integer) | |
| ! n1 >= 2, n1 = power of 2 | |
| ! n2 :data length (integer) | |
| ! n2 >= 2, n2 = power of 2 | |
| ! n3 :data length (integer) | |
| ! n3 >= 2, n3 = power of 2 | |
| ! a(0:n1-1,0:n2-1,0:n3-1) | |
| ! :input/output data (real*8) | |
| ! <case1> | |
| ! input data | |
| ! a(mod(j1,n1),mod(j2,n2),mod(j3,n3)) = A(j1,j2,j3), | |
| ! 0<j1<=n1, 0<j2<=n2, 0<j3<=n3 | |
| ! output data | |
| ! a(k1,k2,k3) = S(k1,k2,k3), | |
| ! 0<=k1<n1, 0<=k2<n2, 0<=k3<n3 | |
| ! <case2> | |
| ! output data | |
| ! a(mod(k1,n1),mod(k2,n2),mod(k3,n3)) = S(k1,k2,k3), | |
| ! 0<k1<=n1, 0<k2<=n2, 0<k3<=n3 | |
| ! t(0:*) :work area (real*8) | |
| ! length of t >= max(4*n2, 4*n3) | |
| ! ip(0:*):work area for bit reversal (integer) | |
| ! length of ip >= 2+sqrt(n) | |
| ! (n = max(n1/2, n2/2, n3/2)) | |
| ! ip(0),ip(1) are pointers of the cos/sin table. | |
| ! w(0:*) :cos/sin table (real*8) | |
| ! length of w >= max(n1*3/2, n2*3/2, n3*3/2) | |
| ! w(),ip() are initialized if ip(0) = 0. | |
| ! [remark] | |
| ! Inverse of | |
| ! call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w) | |
| ! is | |
| ! do j3 = 0, n3 - 1 | |
| ! do j2 = 0, n2 - 1 | |
| ! a(0, j2, j3) = a(0, j2, j3) * 0.5d0 | |
| ! end do | |
| ! do j1 = 0, n1 - 1 | |
| ! a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0 | |
| ! end do | |
| ! end do | |
| ! do j2 = 0, n2 - 1 | |
| ! do j1 = 0, n1 - 1 | |
| ! a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0 | |
| ! end do | |
| ! end do | |
| ! call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w) | |
| ! do j3 = 0, n3 - 1 | |
| ! do j2 = 0, n2 - 1 | |
| ! do j1 = 0, n1 - 1 | |
| ! a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3) | |
| ! end do | |
| ! end do | |
| ! end do | |
| ! . | |
| ! | |
| ! | |
| subroutine cdft3d(n1max, n2max, n1, n2, n3, isgn, a, | |
| & t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), n | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| n = 2 * max(n2, n3) | |
| n = max(n, n1) | |
| if (n .gt. 4 * ip(0)) then | |
| call makewt(n / 4, ip, w) | |
| end if | |
| call xdft3da_sub(n1max, n2max, n1, n2, n3, 0, | |
| & isgn, a, t, ip, w) | |
| call cdft3db_sub(n1max, n2max, n1, n2, n3, | |
| & isgn, a, t, ip, w) | |
| end | |
| ! | |
| subroutine rdft3d(n1max, n2max, n1, n2, n3, isgn, a, | |
| & t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), | |
| & n, nw, nc | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| n = 2 * max(n2, n3) | |
| n = max(n, n1) | |
| nw = ip(0) | |
| if (n .gt. 4 * nw) then | |
| nw = n / 4 | |
| call makewt(nw, ip, w) | |
| end if | |
| nc = ip(1) | |
| if (n1 .gt. 4 * nc) then | |
| nc = n1 / 4 | |
| call makect(nc, ip, w(nw)) | |
| end if | |
| if (isgn .lt. 0) then | |
| call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a) | |
| call cdft3db_sub(n1max, n2max, n1, n2, n3, | |
| & isgn, a, t, ip, w) | |
| call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, | |
| & isgn, a, t, ip, w) | |
| else | |
| call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, | |
| & isgn, a, t, ip, w) | |
| call cdft3db_sub(n1max, n2max, n1, n2, n3, | |
| & isgn, a, t, ip, w) | |
| call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a) | |
| end if | |
| end | |
| ! | |
| subroutine rdft3dsort(n1max, n2max, n1, n2, n3, isgn, a) | |
| integer n1max, n2max, n1, n2, n3, isgn, n2h, n3h, j, k | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), x, y | |
| n2h = n2 / 2 | |
| n3h = n3 / 2 | |
| if (isgn .lt. 0) then | |
| do k = 0, n3 - 1 | |
| do j = n2h + 1, n2 - 1 | |
| a(0, j, k) = a(n1 + 1, j, k) | |
| a(1, j, k) = a(n1, j, k) | |
| end do | |
| end do | |
| do k = n3h + 1, n3 - 1 | |
| a(0, 0, k) = a(n1 + 1, 0, k) | |
| a(1, 0, k) = a(n1, 0, k) | |
| a(0, n2h, k) = a(n1 + 1, n2h, k) | |
| a(1, n2h, k) = a(n1, n2h, k) | |
| end do | |
| a(1, 0, 0) = a(n1, 0, 0) | |
| a(1, n2h, 0) = a(n1, n2h, 0) | |
| a(1, 0, n3h) = a(n1, 0, n3h) | |
| a(1, n2h, n3h) = a(n1, n2h, n3h) | |
| else | |
| do j = n2h + 1, n2 - 1 | |
| y = a(0, j, 0) | |
| x = a(1, j, 0) | |
| a(n1, j, 0) = x | |
| a(n1 + 1, j, 0) = y | |
| a(n1, n2 - j, 0) = x | |
| a(n1 + 1, n2 - j, 0) = -y | |
| a(0, j, 0) = a(0, n2 - j, 0) | |
| a(1, j, 0) = -a(1, n2 - j, 0) | |
| end do | |
| do k = 1, n3 - 1 | |
| do j = n2h + 1, n2 - 1 | |
| y = a(0, j, k) | |
| x = a(1, j, k) | |
| a(n1, j, k) = x | |
| a(n1 + 1, j, k) = y | |
| a(n1, n2 - j, n3 - k) = x | |
| a(n1 + 1, n2 - j, n3 - k) = -y | |
| a(0, j, k) = a(0, n2 - j, n3 - k) | |
| a(1, j, k) = -a(1, n2 - j, n3 - k) | |
| end do | |
| end do | |
| do k = n3h + 1, n3 - 1 | |
| y = a(0, 0, k) | |
| x = a(1, 0, k) | |
| a(n1, 0, k) = x | |
| a(n1 + 1, 0, k) = y | |
| a(n1, 0, n3 - k) = x | |
| a(n1 + 1, 0, n3 - k) = -y | |
| a(0, 0, k) = a(0, 0, n3 - k) | |
| a(1, 0, k) = -a(1, 0, n3 - k) | |
| y = a(0, n2h, k) | |
| x = a(1, n2h, k) | |
| a(n1, n2h, k) = x | |
| a(n1 + 1, n2h, k) = y | |
| a(n1, n2h, n3 - k) = x | |
| a(n1 + 1, n2h, n3 - k) = -y | |
| a(0, n2h, k) = a(0, n2h, n3 - k) | |
| a(1, n2h, k) = -a(1, n2h, n3 - k) | |
| end do | |
| a(n1, 0, 0) = a(1, 0, 0) | |
| a(n1 + 1, 0, 0) = 0 | |
| a(1, 0, 0) = 0 | |
| a(n1, n2h, 0) = a(1, n2h, 0) | |
| a(n1 + 1, n2h, 0) = 0 | |
| a(1, n2h, 0) = 0 | |
| a(n1, 0, n3h) = a(1, 0, n3h) | |
| a(n1 + 1, 0, n3h) = 0 | |
| a(1, 0, n3h) = 0 | |
| a(n1, n2h, n3h) = a(1, n2h, n3h) | |
| a(n1 + 1, n2h, n3h) = 0 | |
| a(1, n2h, n3h) = 0 | |
| end if | |
| end | |
| ! | |
| subroutine ddct3d(n1max, n2max, n1, n2, n3, isgn, a, | |
| & t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), | |
| & n, nw, nc | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| n = max(n2, n3) | |
| n = max(n, n1) | |
| nw = ip(0) | |
| if (n .gt. 4 * nw) then | |
| nw = n / 4 | |
| call makewt(nw, ip, w) | |
| end if | |
| nc = ip(1) | |
| if (n .gt. nc) then | |
| nc = n | |
| call makect(nc, ip, w(nw)) | |
| end if | |
| call ddxt3da_sub(n1max, n2max, n1, n2, n3, 0, | |
| & isgn, a, t, ip, w) | |
| call ddxt3db_sub(n1max, n2max, n1, n2, n3, 0, | |
| & isgn, a, t, ip, w) | |
| end | |
| ! | |
| subroutine ddst3d(n1max, n2max, n1, n2, n3, isgn, a, | |
| & t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), | |
| & n, nw, nc | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| n = max(n2, n3) | |
| n = max(n, n1) | |
| nw = ip(0) | |
| if (n .gt. 4 * nw) then | |
| nw = n / 4 | |
| call makewt(nw, ip, w) | |
| end if | |
| nc = ip(1) | |
| if (n .gt. nc) then | |
| nc = n | |
| call makect(nc, ip, w(nw)) | |
| end if | |
| call ddxt3da_sub(n1max, n2max, n1, n2, n3, 1, | |
| & isgn, a, t, ip, w) | |
| call ddxt3db_sub(n1max, n2max, n1, n2, n3, 1, | |
| & isgn, a, t, ip, w) | |
| end | |
| ! | |
| ! -------- child routines -------- | |
| ! | |
| subroutine xdft3da_sub(n1max, n2max, n1, n2, n3, icr, | |
| & isgn, a, t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, icr, isgn, | |
| & ip(0 : *), i, j, k | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| do k = 0, n3 - 1 | |
| if (icr .eq. 0) then | |
| do j = 0, n2 - 1 | |
| call cdft(n1, isgn, a(0, j, k), ip, w) | |
| end do | |
| else if (isgn .ge. 0) then | |
| do j = 0, n2 - 1 | |
| call rdft(n1, isgn, a(0, j, k), ip, w) | |
| end do | |
| end if | |
| if (n1 .gt. 4) then | |
| do i = 0, n1 - 8, 8 | |
| do j = 0, n2 - 1 | |
| t(2 * j) = a(i, j, k) | |
| t(2 * j + 1) = a(i + 1, j, k) | |
| t(2 * n2 + 2 * j) = a(i + 2, j, k) | |
| t(2 * n2 + 2 * j + 1) = a(i + 3, j, k) | |
| t(4 * n2 + 2 * j) = a(i + 4, j, k) | |
| t(4 * n2 + 2 * j + 1) = a(i + 5, j, k) | |
| t(6 * n2 + 2 * j) = a(i + 6, j, k) | |
| t(6 * n2 + 2 * j + 1) = a(i + 7, j, k) | |
| end do | |
| call cdft(2 * n2, isgn, t, ip, w) | |
| call cdft(2 * n2, isgn, t(2 * n2), ip, w) | |
| call cdft(2 * n2, isgn, t(4 * n2), ip, w) | |
| call cdft(2 * n2, isgn, t(6 * n2), ip, w) | |
| do j = 0, n2 - 1 | |
| a(i, j, k) = t(2 * j) | |
| a(i + 1, j, k) = t(2 * j + 1) | |
| a(i + 2, j, k) = t(2 * n2 + 2 * j) | |
| a(i + 3, j, k) = t(2 * n2 + 2 * j + 1) | |
| a(i + 4, j, k) = t(4 * n2 + 2 * j) | |
| a(i + 5, j, k) = t(4 * n2 + 2 * j + 1) | |
| a(i + 6, j, k) = t(6 * n2 + 2 * j) | |
| a(i + 7, j, k) = t(6 * n2 + 2 * j + 1) | |
| end do | |
| end do | |
| else if (n1 .eq. 4) then | |
| do j = 0, n2 - 1 | |
| t(2 * j) = a(0, j, k) | |
| t(2 * j + 1) = a(1, j, k) | |
| t(2 * n2 + 2 * j) = a(2, j, k) | |
| t(2 * n2 + 2 * j + 1) = a(3, j, k) | |
| end do | |
| call cdft(2 * n2, isgn, t, ip, w) | |
| call cdft(2 * n2, isgn, t(2 * n2), ip, w) | |
| do j = 0, n2 - 1 | |
| a(0, j, k) = t(2 * j) | |
| a(1, j, k) = t(2 * j + 1) | |
| a(2, j, k) = t(2 * n2 + 2 * j) | |
| a(3, j, k) = t(2 * n2 + 2 * j + 1) | |
| end do | |
| else if (n1 .eq. 2) then | |
| do j = 0, n2 - 1 | |
| t(2 * j) = a(0, j, k) | |
| t(2 * j + 1) = a(1, j, k) | |
| end do | |
| call cdft(2 * n2, isgn, t, ip, w) | |
| do j = 0, n2 - 1 | |
| a(0, j, k) = t(2 * j) | |
| a(1, j, k) = t(2 * j + 1) | |
| end do | |
| end if | |
| if (icr .ne. 0 .and. isgn .lt. 0) then | |
| do j = 0, n2 - 1 | |
| call rdft(n1, isgn, a(0, j, k), ip, w) | |
| end do | |
| end if | |
| end do | |
| end | |
| ! | |
| subroutine cdft3db_sub(n1max, n2max, n1, n2, n3, | |
| & isgn, a, t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), | |
| & i, j, k | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| if (n1 .gt. 4) then | |
| do j = 0, n2 - 1 | |
| do i = 0, n1 - 8, 8 | |
| do k = 0, n3 - 1 | |
| t(2 * k) = a(i, j, k) | |
| t(2 * k + 1) = a(i + 1, j, k) | |
| t(2 * n3 + 2 * k) = a(i + 2, j, k) | |
| t(2 * n3 + 2 * k + 1) = a(i + 3, j, k) | |
| t(4 * n3 + 2 * k) = a(i + 4, j, k) | |
| t(4 * n3 + 2 * k + 1) = a(i + 5, j, k) | |
| t(6 * n3 + 2 * k) = a(i + 6, j, k) | |
| t(6 * n3 + 2 * k + 1) = a(i + 7, j, k) | |
| end do | |
| call cdft(2 * n3, isgn, t, ip, w) | |
| call cdft(2 * n3, isgn, t(2 * n3), ip, w) | |
| call cdft(2 * n3, isgn, t(4 * n3), ip, w) | |
| call cdft(2 * n3, isgn, t(6 * n3), ip, w) | |
| do k = 0, n3 - 1 | |
| a(i, j, k) = t(2 * k) | |
| a(i + 1, j, k) = t(2 * k + 1) | |
| a(i + 2, j, k) = t(2 * n3 + 2 * k) | |
| a(i + 3, j, k) = t(2 * n3 + 2 * k + 1) | |
| a(i + 4, j, k) = t(4 * n3 + 2 * k) | |
| a(i + 5, j, k) = t(4 * n3 + 2 * k + 1) | |
| a(i + 6, j, k) = t(6 * n3 + 2 * k) | |
| a(i + 7, j, k) = t(6 * n3 + 2 * k + 1) | |
| end do | |
| end do | |
| end do | |
| else if (n1 .eq. 4) then | |
| do j = 0, n2 - 1 | |
| do k = 0, n3 - 1 | |
| t(2 * k) = a(0, j, k) | |
| t(2 * k + 1) = a(1, j, k) | |
| t(2 * n3 + 2 * k) = a(2, j, k) | |
| t(2 * n3 + 2 * k + 1) = a(3, j, k) | |
| end do | |
| call cdft(2 * n3, isgn, t, ip, w) | |
| call cdft(2 * n3, isgn, t(2 * n3), ip, w) | |
| do k = 0, n3 - 1 | |
| a(0, j, k) = t(2 * k) | |
| a(1, j, k) = t(2 * k + 1) | |
| a(2, j, k) = t(2 * n3 + 2 * k) | |
| a(3, j, k) = t(2 * n3 + 2 * k + 1) | |
| end do | |
| end do | |
| else if (n1 .eq. 2) then | |
| do j = 0, n2 - 1 | |
| do k = 0, n3 - 1 | |
| t(2 * k) = a(0, j, k) | |
| t(2 * k + 1) = a(1, j, k) | |
| end do | |
| call cdft(2 * n3, isgn, t, ip, w) | |
| do k = 0, n3 - 1 | |
| a(0, j, k) = t(2 * k) | |
| a(1, j, k) = t(2 * k + 1) | |
| end do | |
| end do | |
| end if | |
| end | |
| ! | |
| subroutine rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a) | |
| integer n1max, n2max, n1, n2, n3, isgn, | |
| & n2h, n3h, i, j, k, l | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), xi | |
| n2h = n2 / 2 | |
| n3h = n3 / 2 | |
| if (isgn .lt. 0) then | |
| do k = 1, n3h - 1 | |
| l = n3 - k | |
| xi = a(0, 0, k) - a(0, 0, l) | |
| a(0, 0, k) = a(0, 0, k) + a(0, 0, l) | |
| a(0, 0, l) = xi | |
| xi = a(1, 0, l) - a(1, 0, k) | |
| a(1, 0, k) = a(1, 0, k) + a(1, 0, l) | |
| a(1, 0, l) = xi | |
| xi = a(0, n2h, k) - a(0, n2h, l) | |
| a(0, n2h, k) = a(0, n2h, k) + a(0, n2h, l) | |
| a(0, n2h, l) = xi | |
| xi = a(1, n2h, l) - a(1, n2h, k) | |
| a(1, n2h, k) = a(1, n2h, k) + a(1, n2h, l) | |
| a(1, n2h, l) = xi | |
| do i = 1, n2h - 1 | |
| j = n2 - i | |
| xi = a(0, i, k) - a(0, j, l) | |
| a(0, i, k) = a(0, i, k) + a(0, j, l) | |
| a(0, j, l) = xi | |
| xi = a(1, j, l) - a(1, i, k) | |
| a(1, i, k) = a(1, i, k) + a(1, j, l) | |
| a(1, j, l) = xi | |
| xi = a(0, i, l) - a(0, j, k) | |
| a(0, i, l) = a(0, i, l) + a(0, j, k) | |
| a(0, j, k) = xi | |
| xi = a(1, j, k) - a(1, i, l) | |
| a(1, i, l) = a(1, i, l) + a(1, j, k) | |
| a(1, j, k) = xi | |
| end do | |
| end do | |
| do i = 1, n2h - 1 | |
| j = n2 - i | |
| xi = a(0, i, 0) - a(0, j, 0) | |
| a(0, i, 0) = a(0, i, 0) + a(0, j, 0) | |
| a(0, j, 0) = xi | |
| xi = a(1, j, 0) - a(1, i, 0) | |
| a(1, i, 0) = a(1, i, 0) + a(1, j, 0) | |
| a(1, j, 0) = xi | |
| xi = a(0, i, n3h) - a(0, j, n3h) | |
| a(0, i, n3h) = a(0, i, n3h) + a(0, j, n3h) | |
| a(0, j, n3h) = xi | |
| xi = a(1, j, n3h) - a(1, i, n3h) | |
| a(1, i, n3h) = a(1, i, n3h) + a(1, j, n3h) | |
| a(1, j, n3h) = xi | |
| end do | |
| else | |
| do k = 1, n3h - 1 | |
| l = n3 - k | |
| a(0, 0, l) = 0.5d0 * (a(0, 0, k) - a(0, 0, l)) | |
| a(0, 0, k) = a(0, 0, k) - a(0, 0, l) | |
| a(1, 0, l) = 0.5d0 * (a(1, 0, k) + a(1, 0, l)) | |
| a(1, 0, k) = a(1, 0, k) - a(1, 0, l) | |
| a(0, n2h, l) = 0.5d0 * (a(0, n2h, k) - a(0, n2h, l)) | |
| a(0, n2h, k) = a(0, n2h, k) - a(0, n2h, l) | |
| a(1, n2h, l) = 0.5d0 * (a(1, n2h, k) + a(1, n2h, l)) | |
| a(1, n2h, k) = a(1, n2h, k) - a(1, n2h, l) | |
| do i = 1, n2h - 1 | |
| j = n2 - i | |
| a(0, j, l) = 0.5d0 * (a(0, i, k) - a(0, j, l)) | |
| a(0, i, k) = a(0, i, k) - a(0, j, l) | |
| a(1, j, l) = 0.5d0 * (a(1, i, k) + a(1, j, l)) | |
| a(1, i, k) = a(1, i, k) - a(1, j, l) | |
| a(0, j, k) = 0.5d0 * (a(0, i, l) - a(0, j, k)) | |
| a(0, i, l) = a(0, i, l) - a(0, j, k) | |
| a(1, j, k) = 0.5d0 * (a(1, i, l) + a(1, j, k)) | |
| a(1, i, l) = a(1, i, l) - a(1, j, k) | |
| end do | |
| end do | |
| do i = 1, n2h - 1 | |
| j = n2 - i | |
| a(0, j, 0) = 0.5d0 * (a(0, i, 0) - a(0, j, 0)) | |
| a(0, i, 0) = a(0, i, 0) - a(0, j, 0) | |
| a(1, j, 0) = 0.5d0 * (a(1, i, 0) + a(1, j, 0)) | |
| a(1, i, 0) = a(1, i, 0) - a(1, j, 0) | |
| a(0, j, n3h) = 0.5d0 * (a(0, i, n3h) - a(0, j, n3h)) | |
| a(0, i, n3h) = a(0, i, n3h) - a(0, j, n3h) | |
| a(1, j, n3h) = 0.5d0 * (a(1, i, n3h) + a(1, j, n3h)) | |
| a(1, i, n3h) = a(1, i, n3h) - a(1, j, n3h) | |
| end do | |
| end if | |
| end | |
| ! | |
| subroutine ddxt3da_sub(n1max, n2max, n1, n2, n3, ics, | |
| & isgn, a, t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, ics, isgn, | |
| & ip(0 : *), i, j, k | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| do k = 0, n3 - 1 | |
| if (ics .eq. 0) then | |
| do j = 0, n2 - 1 | |
| call ddct(n1, isgn, a(0, j, k), ip, w) | |
| end do | |
| else | |
| do j = 0, n2 - 1 | |
| call ddst(n1, isgn, a(0, j, k), ip, w) | |
| end do | |
| end if | |
| if (n1 .gt. 2) then | |
| do i = 0, n1 - 4, 4 | |
| do j = 0, n2 - 1 | |
| t(j) = a(i, j, k) | |
| t(n2 + j) = a(i + 1, j, k) | |
| t(2 * n2 + j) = a(i + 2, j, k) | |
| t(3 * n2 + j) = a(i + 3, j, k) | |
| end do | |
| if (ics .eq. 0) then | |
| call ddct(n2, isgn, t, ip, w) | |
| call ddct(n2, isgn, t(n2), ip, w) | |
| call ddct(n2, isgn, t(2 * n2), ip, w) | |
| call ddct(n2, isgn, t(3 * n2), ip, w) | |
| else | |
| call ddst(n2, isgn, t, ip, w) | |
| call ddst(n2, isgn, t(n2), ip, w) | |
| call ddst(n2, isgn, t(2 * n2), ip, w) | |
| call ddst(n2, isgn, t(3 * n2), ip, w) | |
| end if | |
| do j = 0, n2 - 1 | |
| a(i, j, k) = t(j) | |
| a(i + 1, j, k) = t(n2 + j) | |
| a(i + 2, j, k) = t(2 * n2 + j) | |
| a(i + 3, j, k) = t(3 * n2 + j) | |
| end do | |
| end do | |
| else if (n1 .eq. 2) then | |
| do j = 0, n2 - 1 | |
| t(j) = a(0, j, k) | |
| t(n2 + j) = a(1, j, k) | |
| end do | |
| if (ics .eq. 0) then | |
| call ddct(n2, isgn, t, ip, w) | |
| call ddct(n2, isgn, t(n2), ip, w) | |
| else | |
| call ddst(n2, isgn, t, ip, w) | |
| call ddst(n2, isgn, t(n2), ip, w) | |
| end if | |
| do j = 0, n2 - 1 | |
| a(0, j, k) = t(j) | |
| a(1, j, k) = t(n2 + j) | |
| end do | |
| end if | |
| end do | |
| end | |
| ! | |
| subroutine ddxt3db_sub(n1max, n2max, n1, n2, n3, ics, | |
| & isgn, a, t, ip, w) | |
| integer n1max, n2max, n1, n2, n3, ics, isgn, | |
| & ip(0 : *), i, j, k | |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), | |
| & t(0 : *), w(0 : *) | |
| if (n1 .gt. 2) then | |
| do j = 0, n2 - 1 | |
| do i = 0, n1 - 4, 4 | |
| do k = 0, n3 - 1 | |
| t(k) = a(i, j, k) | |
| t(n3 + k) = a(i + 1, j, k) | |
| t(2 * n3 + k) = a(i + 2, j, k) | |
| t(3 * n3 + k) = a(i + 3, j, k) | |
| end do | |
| if (ics .eq. 0) then | |
| call ddct(n3, isgn, t, ip, w) | |
| call ddct(n3, isgn, t(n3), ip, w) | |
| call ddct(n3, isgn, t(2 * n3), ip, w) | |
| call ddct(n3, isgn, t(3 * n3), ip, w) | |
| else | |
| call ddst(n3, isgn, t, ip, w) | |
| call ddst(n3, isgn, t(n3), ip, w) | |
| call ddst(n3, isgn, t(2 * n3), ip, w) | |
| call ddst(n3, isgn, t(3 * n3), ip, w) | |
| end if | |
| do k = 0, n3 - 1 | |
| a(i, j, k) = t(k) | |
| a(i + 1, j, k) = t(n3 + k) | |
| a(i + 2, j, k) = t(2 * n3 + k) | |
| a(i + 3, j, k) = t(3 * n3 + k) | |
| end do | |
| end do | |
| end do | |
| else if (n1 .eq. 2) then | |
| do j = 0, n2 - 1 | |
| do k = 0, n3 - 1 | |
| t(k) = a(0, j, k) | |
| t(n3 + k) = a(1, j, k) | |
| end do | |
| if (ics .eq. 0) then | |
| call ddct(n3, isgn, t, ip, w) | |
| call ddct(n3, isgn, t(n3), ip, w) | |
| else | |
| call ddst(n3, isgn, t, ip, w) | |
| call ddst(n3, isgn, t(n3), ip, w) | |
| end if | |
| do k = 0, n3 - 1 | |
| a(0, j, k) = t(k) | |
| a(1, j, k) = t(n3 + k) | |
| end do | |
| end do | |
| end if | |
| end | |
| ! |