third_party/tensorflow_dependencies/fft2d/fftsg2d.f - aom - Git at Google


            
            ! Fast Fourier/Cosine/Sine Transform

            
            !     dimension   :two

            
            !     data length :power of 2

            
            !     decimation  :frequency

            
            !     radix       :split-radix, row-column

            
            !     data        :inplace

            
            !     table       :use

            
            ! subroutines

            
            !     cdft2d: Complex Discrete Fourier Transform

            
            !     rdft2d: Real Discrete Fourier Transform

            
            !     ddct2d: Discrete Cosine Transform

            
            !     ddst2d: Discrete Sine Transform

            
            ! necessary package

            
            !     fftsg.f  : 1D-FFT package

            
            !

            
            !

            
            ! -------- Complex DFT (Discrete Fourier Transform) --------

            
            !     [definition]

            
            !         <case1>

            
            !             X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) *

            
            !                            exp(2*pi*i*j1*k1/n1) *

            
            !                            exp(2*pi*i*j2*k2/n2),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !         <case2>

            
            !             X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) *

            
            !                            exp(-2*pi*i*j1*k1/n1) *

            
            !                            exp(-2*pi*i*j2*k2/n2),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)

            
            !     [usage]

            
            !         <case1>

            
            !             ip(0) = 0  ! first time only

            
            !             call cdft2d(n1max, 2*n1, n2, 1, a, t, ip, w)

            
            !         <case2>

            
            !             ip(0) = 0  ! first time only

            
            !             call cdft2d(n1max, 2*n1, n2, -1, a, t, ip, w)

            
            !     [parameters]

            
            !         n1max  :row size of the 2D array (integer)

            
            !         2*n1   :data length (integer)

            
            !                 n1 >= 1, n1 = power of 2

            
            !         n2     :data length (integer)

            
            !                 n2 >= 1, n2 = power of 2

            
            !         a(0:2*n1-1,0:n2-1)

            
            !                :input/output data (real*8)

            
            !                 input data

            
            !                     a(2*j1,j2) = Re(x(j1,j2)),

            
            !                     a(2*j1+1,j2) = Im(x(j1,j2)),

            
            !                     0<=j1<n1, 0<=j2<n2

            
            !                 output data

            
            !                     a(2*k1,k2) = Re(X(k1,k2)),

            
            !                     a(2*k1+1,k2) = Im(X(k1,k2)),

            
            !                     0<=k1<n1, 0<=k2<n2

            
            !         t(0:8*n2-1)

            
            !                :work area (real*8)

            
            !                 length of t >= 8*n2

            
            !         ip(0:*):work area for bit reversal (integer)

            
            !                 length of ip >= 2+sqrt(n)

            
            !                 (n = max(n1, n2))

            
            !                 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)

            
            !                 w(),ip() are initialized if ip(0) = 0.

            
            !     [remark]

            
            !         Inverse of

            
            !             call cdft2d(n1max, 2*n1, n2, -1, a, t, ip, w)

            
            !         is

            
            !             call cdft2d(n1max, 2*n1, n2, 1, a, t, ip, w)

            
            !             do j2 = 0, n2 - 1

            
            !                 do j1 = 0, 2 * n1 - 1

            
            !                     a(j1, j2) = a(j1, j2) * (1.0d0 / n1 / n2)

            
            !                 end do

            
            !             end do

            
            !         .

            
            !

            
            !

            
            ! -------- Real DFT / Inverse of Real DFT --------

            
            !     [definition]

            
            !         <case1> RDFT

            
            !             R(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *

            
            !                            cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !             I(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *

            
            !                            sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !         <case2> IRDFT (excluding scale)

            
            !             a(k1,k2) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1

            
            !                            (R(j1,j2) *

            
            !                            cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2) +

            
            !                            I(j1,j2) *

            
            !                            sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2)),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !         (notes: R(n1-k1,n2-k2) = R(k1,k2),

            
            !                 I(n1-k1,n2-k2) = -I(k1,k2),

            
            !                 R(n1-k1,0) = R(k1,0),

            
            !                 I(n1-k1,0) = -I(k1,0),

            
            !                 R(0,n2-k2) = R(0,k2),

            
            !                 I(0,n2-k2) = -I(0,k2),

            
            !                 0<k1<n1, 0<k2<n2)

            
            !     [usage]

            
            !         <case1>

            
            !             ip(0) = 0  ! first time only

            
            !             call rdft2d(n1max, n1, n2, 1, a, t, ip, w)

            
            !         <case2>

            
            !             ip(0) = 0  ! first time only

            
            !             call rdft2d(n1max, n1, n2, -1, a, t, ip, w)

            
            !     [parameters]

            
            !         n1max  :row size of the 2D array (integer)

            
            !         n1     :data length (integer)

            
            !                 n1 >= 2, n1 = power of 2

            
            !         n2     :data length (integer)

            
            !                 n2 >= 2, n2 = power of 2

            
            !         a(0:n1-1,0:n2-1)

            
            !                :input/output data (real*8)

            
            !                 <case1>

            
            !                     output data

            
            !                         a(2*k1,k2) = R(k1,k2) = R(n1-k1,n2-k2),

            
            !                         a(2*k1+1,k2) = I(k1,k2) = -I(n1-k1,n2-k2),

            
            !                             0<k1<n1/2, 0<k2<n2,

            
            !                         a(2*k1,0) = R(k1,0) = R(n1-k1,0),

            
            !                         a(2*k1+1,0) = I(k1,0) = -I(n1-k1,0),

            
            !                             0<k1<n1/2,

            
            !                         a(0,k2) = R(0,k2) = R(0,n2-k2),

            
            !                         a(1,k2) = I(0,k2) = -I(0,n2-k2),

            
            !                         a(1,n2-k2) = R(n1/2,k2) = R(n1/2,n2-k2),

            
            !                         a(0,n2-k2) = -I(n1/2,k2) = I(n1/2,n2-k2),

            
            !                             0<k2<n2/2,

            
            !                         a(0,0) = R(0,0),

            
            !                         a(1,0) = R(n1/2,0),

            
            !                         a(0,n2/2) = R(0,n2/2),

            
            !                         a(1,n2/2) = R(n1/2,n2/2)

            
            !                 <case2>

            
            !                     input data

            
            !                         a(2*j1,j2) = R(j1,j2) = R(n1-j1,n2-j2),

            
            !                         a(2*j1+1,j2) = I(j1,j2) = -I(n1-j1,n2-j2),

            
            !                             0<j1<n1/2, 0<j2<n2,

            
            !                         a(2*j1,0) = R(j1,0) = R(n1-j1,0),

            
            !                         a(2*j1+1,0) = I(j1,0) = -I(n1-j1,0),

            
            !                             0<j1<n1/2,

            
            !                         a(0,j2) = R(0,j2) = R(0,n2-j2),

            
            !                         a(1,j2) = I(0,j2) = -I(0,n2-j2),

            
            !                         a(1,n2-j2) = R(n1/2,j2) = R(n1/2,n2-j2),

            
            !                         a(0,n2-j2) = -I(n1/2,j2) = I(n1/2,n2-j2),

            
            !                             0<j2<n2/2,

            
            !                         a(0,0) = R(0,0),

            
            !                         a(1,0) = R(n1/2,0),

            
            !                         a(0,n2/2) = R(0,n2/2),

            
            !                         a(1,n2/2) = R(n1/2,n2/2)

            
            !                 ---- output ordering ----

            
            !                     call rdft2d(n1max, n1, n2, 1, a, t, ip, w)

            
            !                     call rdft2dsort(n1max, n1, n2, 1, a)

            
            !                     ! stored data is a(0:n1-1,0:n2+1):

            
            !                     ! a(2*k1,k2) = R(k1,k2),

            
            !                     ! a(2*k1+1,k2) = I(k1,k2),

            
            !                     ! 0<=k1<=n1/2, 0<=k2<n2.

            
            !                     ! the stored data is larger than the input data!

            
            !                 ---- input ordering ----

            
            !                     call rdft2dsort(n1max, n1, n2, -1, a)

            
            !                     call rdft2d(n1max, n1, n2, -1, a, t, ip, w)

            
            !         t(0:8*n2-1)

            
            !                :work area (real*8)

            
            !                 length of t >= 8*n2

            
            !         ip(0:*):work area for bit reversal (integer)

            
            !                 length of ip >= 2+sqrt(n)

            
            !                 (n = max(n1/2, n2))

            
            !                 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) + n1/4

            
            !                 w(),ip() are initialized if ip(0) = 0.

            
            !     [remark]

            
            !         Inverse of

            
            !             call rdft2d(n1max, n1, n2, 1, a, t, ip, w)

            
            !         is

            
            !             call rdft2d(n1max, n1, n2, -1, a, t, ip, w)

            
            !             do j2 = 0, n2 - 1

            
            !                 do j1 = 0, n1 - 1

            
            !                     a(j1, j2) = a(j1, j2) * (2.0d0 / n1 / n2)

            
            !                 end do

            
            !             end do

            
            !         .

            
            !

            
            !

            
            ! -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------

            
            !     [definition]

            
            !         <case1> IDCT (excluding scale)

            
            !             C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *

            
            !                            cos(pi*j1*(k1+1/2)/n1) *

            
            !                            cos(pi*j2*(k2+1/2)/n2),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !         <case2> DCT

            
            !             C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *

            
            !                            cos(pi*(j1+1/2)*k1/n1) *

            
            !                            cos(pi*(j2+1/2)*k2/n2),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !     [usage]

            
            !         <case1>

            
            !             ip(0) = 0  ! first time only

            
            !             call ddct2d(n1max, n1, n2, 1, a, t, ip, w)

            
            !         <case2>

            
            !             ip(0) = 0  ! first time only

            
            !             call ddct2d(n1max, n1, n2, -1, a, t, ip, w)

            
            !     [parameters]

            
            !         n1max  :row size of the 2D array (integer)

            
            !         n1     :data length (integer)

            
            !                 n1 >= 2, n1 = power of 2

            
            !         n2     :data length (integer)

            
            !                 n2 >= 2, n2 = power of 2

            
            !         a(0:n1-1,0:n2-1)

            
            !                :input/output data (real*8)

            
            !                 output data

            
            !                     a(k1,k2) = C(k1,k2), 0<=k1<n1, 0<=k2<n2

            
            !         t(0:4*n2-1)

            
            !                :work area (real*8)

            
            !                 length of t >= 4*n2

            
            !         ip(0:*):work area for bit reversal (integer)

            
            !                 length of ip >= 2+sqrt(n)

            
            !                 (n = max(n1/2, n2/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)

            
            !                 w(),ip() are initialized if ip(0) = 0.

            
            !     [remark]

            
            !         Inverse of

            
            !             call ddct2d(n1max, n1, n2, -1, a, t, ip, w)

            
            !         is

            
            !             do j1 = 0, n1 - 1

            
            !                 a(j1, 0) = a(j1, 0) * 0.5d0

            
            !             end do

            
            !             do j2 = 0, n2 - 1

            
            !                 a(0, j2) = a(0, j2) * 0.5d0

            
            !             end do

            
            !             call ddct2d(n1max, n1, n2, 1, a, t, ip, w)

            
            !             do j2 = 0, n2 - 1

            
            !                 do j1 = 0, n1 - 1

            
            !                     a(j1, j2) = a(j1, j2) * (4.0d0 / n1 / n2)

            
            !                 end do

            
            !             end do

            
            !         .

            
            !

            
            !

            
            ! -------- DST (Discrete Sine Transform) / Inverse of DST --------

            
            !     [definition]

            
            !         <case1> IDST (excluding scale)

            
            !             S(k1,k2) = sum_j1=1^n1 sum_j2=1^n2 A(j1,j2) *

            
            !                            sin(pi*j1*(k1+1/2)/n1) *

            
            !                            sin(pi*j2*(k2+1/2)/n2),

            
            !                            0<=k1<n1, 0<=k2<n2

            
            !         <case2> DST

            
            !             S(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *

            
            !                            sin(pi*(j1+1/2)*k1/n1) *

            
            !                            sin(pi*(j2+1/2)*k2/n2),

            
            !                            0<k1<=n1, 0<k2<=n2

            
            !     [usage]

            
            !         <case1>

            
            !             ip(0) = 0  ! first time only

            
            !             call ddst2d(n1max, n1, n2, 1, a, t, ip, w)

            
            !         <case2>

            
            !             ip(0) = 0  ! first time only

            
            !             call ddst2d(n1max, n1, n2, -1, a, t, ip, w)

            
            !     [parameters]

            
            !         n1max  :row size of the 2D array (integer)

            
            !         n1     :data length (integer)

            
            !                 n1 >= 2, n1 = power of 2

            
            !         n2     :data length (integer)

            
            !                 n2 >= 2, n2 = power of 2

            
            !         a(0:n1-1,0:n2-1)

            
            !                :input/output data (real*8)

            
            !                 <case1>

            
            !                     input data

            
            !                         a(j1,j2) = A(j1,j2), 0<j1<n1, 0<j2<n2,

            
            !                         a(j1,0) = A(j1,n2), 0<j1<n1,

            
            !                         a(0,j2) = A(n1,j2), 0<j2<n2,

            
            !                         a(0,0) = A(n1,n2)

            
            !                         (i.e. A(j1,j2) = a(mod(j1,n1),mod(j2,n2)))

            
            !                     output data

            
            !                         a(k1,k2) = S(k1,k2), 0<=k1<n1, 0<=k2<n2

            
            !                 <case2>

            
            !                     output data

            
            !                         a(k1,k2) = S(k1,k2), 0<k1<n1, 0<k2<n2,

            
            !                         a(k1,0) = S(k1,n2), 0<k1<n1,

            
            !                         a(0,k2) = S(n1,k2), 0<k2<n2,

            
            !                         a(0,0) = S(n1,n2)

            
            !                         (i.e. S(k1,k2) = a(mod(k1,n1),mod(k2,n2)))

            
            !         t(0:4*n2-1)

            
            !                :work area (real*8)

            
            !                 length of t >= 4*n2

            
            !         ip(0:*):work area for bit reversal (integer)

            
            !                 length of ip >= 2+sqrt(n)

            
            !                 (n = max(n1/2, n2/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)

            
            !                 w(),ip() are initialized if ip(0) = 0.

            
            !     [remark]

            
            !         Inverse of

            
            !             call ddst2d(n1max, n1, n2, -1, a, t, ip, w)

            
            !         is

            
            !             do j1 = 0, n1 - 1

            
            !                 a(j1, 0) = a(j1, 0) * 0.5d0

            
            !             end do

            
            !             do j2 = 0, n2 - 1

            
            !                 a(0, j2) = a(0, j2) * 0.5d0

            
            !             end do

            
            !             call ddst2d(n1max, n1, n2, 1, a, t, ip, w)

            
            !             do j2 = 0, n2 - 1

            
            !                 do j1 = 0, n1 - 1

            
            !                     a(j1, j2) = a(j1, j2) * (4.0d0 / n1 / n2)

            
            !                 end do

            
            !             end do

            
            !         .

            
            !

            
            !

            
                  subroutine cdft2d(n1max, n1, n2, isgn, a, t, ip, w)

            
                  integer n1max, n1, n2, isgn, ip(0 : *), n, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 * n2 - 1),

            
                 &    w(0 : *)

            
                  n = max(n1, 2 * n2)

            
                  if (n .gt. 4 * ip(0)) then

            
                      call makewt(n / 4, ip, w)

            
                  end if

            
                  do j = 0, n2 - 1

            
                      call cdft(n1, isgn, a(0, j), ip, w)

            
                  end do

            
                  call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)

            
                  end

            
            !

            
                  subroutine rdft2d(n1max, n1, n2, isgn, a, t, ip, w)

            
                  integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 * n2 - 1),

            
                 &    w(0 : *)

            
                  n = max(n1, 2 * n2)

            
                  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 rdft2d_sub(n1max, n1, n2, isgn, a)

            
                      call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)

            
                  end if

            
                  do j = 0, n2 - 1

            
                      call rdft(n1, isgn, a(0, j), ip, w)

            
                  end do

            
                  if (isgn .ge. 0) then

            
                      call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)

            
                      call rdft2d_sub(n1max, n1, n2, isgn, a)

            
                  end if

            
                  end

            
            !

            
                  subroutine rdft2dsort(n1max, n1, n2, isgn, a)

            
                  integer n1max, n1, n2, isgn, n2h, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), x, y

            
                  n2h = n2 / 2

            
                  if (isgn .lt. 0) then

            
                      do j = n2h + 1, n2 - 1

            
                          a(0, j) = a(n1 + 1, j)

            
                          a(1, j) = a(n1, j)

            
                      end do

            
                      a(1, 0) = a(n1, 0)

            
                      a(1, n2h) = a(n1, n2h)

            
                  else

            
                      do j = n2h + 1, n2 - 1

            
                          y = a(0, j)

            
                          x = a(1, j)

            
                          a(n1, j) = x

            
                          a(n1 + 1, j) = y

            
                          a(n1, n2 - j) = x

            
                          a(n1 + 1, n2 - j) = -y

            
                          a(0, j) = a(0, n2 - j)

            
                          a(1, j) = -a(1, n2 - j)

            
                      end do

            
                      a(n1, 0) = a(1, 0)

            
                      a(n1 + 1, 0) = 0

            
                      a(1, 0) = 0

            
                      a(n1, n2h) = a(1, n2h)

            
                      a(n1 + 1, n2h) = 0

            
                      a(1, n2h) = 0

            
                  end if

            
                  end

            
            !

            
                  subroutine ddct2d(n1max, n1, n2, isgn, a, t, ip, w)

            
                  integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 * n2 - 1),

            
                 &    w(0 : *)

            
                  n = max(n1, n2)

            
                  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

            
                  do j = 0, n2 - 1

            
                      call ddct(n1, isgn, a(0, j), ip, w)

            
                  end do

            
                  call ddxt2d_sub(n1max, n1, n2, 0, isgn, a, t, ip, w)

            
                  end

            
            !

            
                  subroutine ddst2d(n1max, n1, n2, isgn, a, t, ip, w)

            
                  integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 * n2 - 1),

            
                 &    w(0 : *)

            
                  n = max(n1, n2)

            
                  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

            
                  do j = 0, n2 - 1

            
                      call ddst(n1, isgn, a(0, j), ip, w)

            
                  end do

            
                  call ddxt2d_sub(n1max, n1, n2, 1, isgn, a, t, ip, w)

            
                  end

            
            !

            
            ! -------- child routines --------

            
            !

            
                  subroutine cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)

            
                  integer n1max, n1, n2, isgn, ip(0 : *), i, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 * n2 - 1),

            
                 &    w(0 : *)

            
                  if (n1 .gt. 4) then

            
                      do i = 0, n1 - 8, 8

            
                          do j = 0, n2 - 1

            
                              t(2 * j) = a(i, j)

            
                              t(2 * j + 1) = a(i + 1, j)

            
                              t(2 * n2 + 2 * j) = a(i + 2, j)

            
                              t(2 * n2 + 2 * j + 1) = a(i + 3, j)

            
                              t(4 * n2 + 2 * j) = a(i + 4, j)

            
                              t(4 * n2 + 2 * j + 1) = a(i + 5, j)

            
                              t(6 * n2 + 2 * j) = a(i + 6, j)

            
                              t(6 * n2 + 2 * j + 1) = a(i + 7, j)

            
                          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) = t(2 * j)

            
                              a(i + 1, j) = t(2 * j + 1)

            
                              a(i + 2, j) = t(2 * n2 + 2 * j)

            
                              a(i + 3, j) = t(2 * n2 + 2 * j + 1)

            
                              a(i + 4, j) = t(4 * n2 + 2 * j)

            
                              a(i + 5, j) = t(4 * n2 + 2 * j + 1)

            
                              a(i + 6, j) = t(6 * n2 + 2 * j)

            
                              a(i + 7, j) = 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)

            
                          t(2 * j + 1) = a(1, j)

            
                          t(2 * n2 + 2 * j) = a(2, j)

            
                          t(2 * n2 + 2 * j + 1) = a(3, j)

            
                      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) = t(2 * j)

            
                          a(1, j) = t(2 * j + 1)

            
                          a(2, j) = t(2 * n2 + 2 * j)

            
                          a(3, j) = 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)

            
                          t(2 * j + 1) = a(1, j)

            
                      end do

            
                      call cdft(2 * n2, isgn, t, ip, w)

            
                      do j = 0, n2 - 1

            
                          a(0, j) = t(2 * j)

            
                          a(1, j) = t(2 * j + 1)

            
                      end do

            
                  end if

            
                  end

            
            !

            
                  subroutine rdft2d_sub(n1max, n1, n2, isgn, a)

            
                  integer n1max, n1, n2, isgn, n2h, i, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), xi

            
                  n2h = n2 / 2

            
                  if (isgn .lt. 0) then

            
                      do i = 1, n2h - 1

            
                          j = n2 - i

            
                          xi = a(0, i) - a(0, j)

            
                          a(0, i) = a(0, i) + a(0, j)

            
                          a(0, j) = xi

            
                          xi = a(1, j) - a(1, i)

            
                          a(1, i) = a(1, i) + a(1, j)

            
                          a(1, j) = xi

            
                      end do

            
                  else

            
                      do i = 1, n2h - 1

            
                          j = n2 - i

            
                          a(0, j) = 0.5d0 * (a(0, i) - a(0, j))

            
                          a(0, i) = a(0, i) - a(0, j)

            
                          a(1, j) = 0.5d0 * (a(1, i) + a(1, j))

            
                          a(1, i) = a(1, i) - a(1, j)

            
                      end do

            
                  end if

            
                  end

            
            !

            
                  subroutine ddxt2d_sub(n1max, n1, n2, ics, isgn, a, t,

            
                 &    ip, w)

            
                  integer n1max, n1, n2, ics, isgn, ip(0 : *), i, j

            
                  real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 * n2 - 1),

            
                 &    w(0 : *)

            
                  if (n1 .gt. 2) then

            
                      do i = 0, n1 - 4, 4

            
                          do j = 0, n2 - 1

            
                              t(j) = a(i, j)

            
                              t(n2 + j) = a(i + 1, j)

            
                              t(2 * n2 + j) = a(i + 2, j)

            
                              t(3 * n2 + j) = a(i + 3, j)

            
                          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) = t(j)

            
                              a(i + 1, j) = t(n2 + j)

            
                              a(i + 2, j) = t(2 * n2 + j)

            
                              a(i + 3, j) = t(3 * n2 + j)

            
                          end do

            
                      end do

            
                  else if (n1 .eq. 2) then

            
                      do j = 0, n2 - 1

            
                          t(j) = a(0, j)

            
                          t(n2 + j) = a(1, j)

            
                      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) = t(j)

            
                          a(1, j) = t(n2 + j)

            
                      end do

            
                  end if

            
                  end

            
            !
	! Fast Fourier/Cosine/Sine Transform
	! dimension :two
	! data length :power of 2
	! decimation :frequency
	! radix :split-radix, row-column
	! data :inplace
	! table :use
	! subroutines
	! cdft2d: Complex Discrete Fourier Transform
	! rdft2d: Real Discrete Fourier Transform
	! ddct2d: Discrete Cosine Transform
	! ddst2d: Discrete Sine Transform
	! necessary package
	! fftsg.f : 1D-FFT package
	!
	!
	! -------- Complex DFT (Discrete Fourier Transform) --------
	! [definition]
	! <case1>
	! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) *
	! exp(2piij1k1/n1) *
	! exp(2piij2k2/n2),
	! 0<=k1<n1, 0<=k2<n2
	! <case2>
	! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) *
	! exp(-2piij1k1/n1) *
	! exp(-2piij2k2/n2),
	! 0<=k1<n1, 0<=k2<n2
	! (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
	! [usage]
	! <case1>
	! ip(0) = 0 ! first time only
	! call cdft2d(n1max, 2*n1, n2, 1, a, t, ip, w)
	! <case2>
	! ip(0) = 0 ! first time only
	! call cdft2d(n1max, 2*n1, n2, -1, a, t, ip, w)
	! [parameters]
	! n1max :row size of the 2D array (integer)
	! 2*n1 :data length (integer)
	! n1 >= 1, n1 = power of 2
	! n2 :data length (integer)
	! n2 >= 1, n2 = power of 2
	! a(0:2*n1-1,0:n2-1)
	! :input/output data (real*8)
	! input data
	! a(2*j1,j2) = Re(x(j1,j2)),
	! a(2*j1+1,j2) = Im(x(j1,j2)),
	! 0<=j1<n1, 0<=j2<n2
	! output data
	! a(2*k1,k2) = Re(X(k1,k2)),
	! a(2*k1+1,k2) = Im(X(k1,k2)),
	! 0<=k1<n1, 0<=k2<n2
	! t(0:8*n2-1)
	! :work area (real*8)
	! length of t >= 8*n2
	! ip(0:*):work area for bit reversal (integer)
	! length of ip >= 2+sqrt(n)
	! (n = max(n1, n2))
	! ip(0),ip(1) are pointers of the cos/sin table.
	! w(0:) :cos/sin table (real8)
	! length of w >= max(n1/2, n2/2)
	! w(),ip() are initialized if ip(0) = 0.
	! [remark]
	! Inverse of
	! call cdft2d(n1max, 2*n1, n2, -1, a, t, ip, w)
	! is
	! call cdft2d(n1max, 2*n1, n2, 1, a, t, ip, w)
	! do j2 = 0, n2 - 1
	! do j1 = 0, 2 * n1 - 1
	! a(j1, j2) = a(j1, j2) * (1.0d0 / n1 / n2)
	! end do
	! end do
	! .
	!
	!
	! -------- Real DFT / Inverse of Real DFT --------
	! [definition]
	! <case1> RDFT
	! R(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *
	! cos(2pij1k1/n1 + 2pij2k2/n2),
	! 0<=k1<n1, 0<=k2<n2
	! I(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *
	! sin(2pij1k1/n1 + 2pij2k2/n2),
	! 0<=k1<n1, 0<=k2<n2
	! <case2> IRDFT (excluding scale)
	! a(k1,k2) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1
	! (R(j1,j2) *
	! cos(2pij1k1/n1 + 2pij2k2/n2) +
	! I(j1,j2) *
	! sin(2pij1k1/n1 + 2pij2k2/n2)),
	! 0<=k1<n1, 0<=k2<n2
	! (notes: R(n1-k1,n2-k2) = R(k1,k2),
	! I(n1-k1,n2-k2) = -I(k1,k2),
	! R(n1-k1,0) = R(k1,0),
	! I(n1-k1,0) = -I(k1,0),
	! R(0,n2-k2) = R(0,k2),
	! I(0,n2-k2) = -I(0,k2),
	! 0<k1<n1, 0<k2<n2)
	! [usage]
	! <case1>
	! ip(0) = 0 ! first time only
	! call rdft2d(n1max, n1, n2, 1, a, t, ip, w)
	! <case2>
	! ip(0) = 0 ! first time only
	! call rdft2d(n1max, n1, n2, -1, a, t, ip, w)
	! [parameters]
	! n1max :row size of the 2D array (integer)
	! n1 :data length (integer)
	! n1 >= 2, n1 = power of 2
	! n2 :data length (integer)
	! n2 >= 2, n2 = power of 2
	! a(0:n1-1,0:n2-1)
	! :input/output data (real*8)
	! <case1>
	! output data
	! a(2*k1,k2) = R(k1,k2) = R(n1-k1,n2-k2),
	! a(2*k1+1,k2) = I(k1,k2) = -I(n1-k1,n2-k2),
	! 0<k1<n1/2, 0<k2<n2,
	! a(2*k1,0) = R(k1,0) = R(n1-k1,0),
	! a(2*k1+1,0) = I(k1,0) = -I(n1-k1,0),
	! 0<k1<n1/2,
	! a(0,k2) = R(0,k2) = R(0,n2-k2),
	! a(1,k2) = I(0,k2) = -I(0,n2-k2),
	! a(1,n2-k2) = R(n1/2,k2) = R(n1/2,n2-k2),
	! a(0,n2-k2) = -I(n1/2,k2) = I(n1/2,n2-k2),
	! 0<k2<n2/2,
	! a(0,0) = R(0,0),
	! a(1,0) = R(n1/2,0),
	! a(0,n2/2) = R(0,n2/2),
	! a(1,n2/2) = R(n1/2,n2/2)
	! <case2>
	! input data
	! a(2*j1,j2) = R(j1,j2) = R(n1-j1,n2-j2),
	! a(2*j1+1,j2) = I(j1,j2) = -I(n1-j1,n2-j2),
	! 0<j1<n1/2, 0<j2<n2,
	! a(2*j1,0) = R(j1,0) = R(n1-j1,0),
	! a(2*j1+1,0) = I(j1,0) = -I(n1-j1,0),
	! 0<j1<n1/2,
	! a(0,j2) = R(0,j2) = R(0,n2-j2),
	! a(1,j2) = I(0,j2) = -I(0,n2-j2),
	! a(1,n2-j2) = R(n1/2,j2) = R(n1/2,n2-j2),
	! a(0,n2-j2) = -I(n1/2,j2) = I(n1/2,n2-j2),
	! 0<j2<n2/2,
	! a(0,0) = R(0,0),
	! a(1,0) = R(n1/2,0),
	! a(0,n2/2) = R(0,n2/2),
	! a(1,n2/2) = R(n1/2,n2/2)
	! ---- output ordering ----
	! call rdft2d(n1max, n1, n2, 1, a, t, ip, w)
	! call rdft2dsort(n1max, n1, n2, 1, a)
	! ! stored data is a(0:n1-1,0:n2+1):
	! ! a(2*k1,k2) = R(k1,k2),
	! ! a(2*k1+1,k2) = I(k1,k2),
	! ! 0<=k1<=n1/2, 0<=k2<n2.
	! ! the stored data is larger than the input data!
	! ---- input ordering ----
	! call rdft2dsort(n1max, n1, n2, -1, a)
	! call rdft2d(n1max, n1, n2, -1, a, t, ip, w)
	! t(0:8*n2-1)
	! :work area (real*8)
	! length of t >= 8*n2
	! ip(0:*):work area for bit reversal (integer)
	! length of ip >= 2+sqrt(n)
	! (n = max(n1/2, n2))
	! ip(0),ip(1) are pointers of the cos/sin table.
	! w(0:) :cos/sin table (real8)
	! length of w >= max(n1/4, n2/2) + n1/4
	! w(),ip() are initialized if ip(0) = 0.
	! [remark]
	! Inverse of
	! call rdft2d(n1max, n1, n2, 1, a, t, ip, w)
	! is
	! call rdft2d(n1max, n1, n2, -1, a, t, ip, w)
	! do j2 = 0, n2 - 1
	! do j1 = 0, n1 - 1
	! a(j1, j2) = a(j1, j2) * (2.0d0 / n1 / n2)
	! end do
	! end do
	! .
	!
	!
	! -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
	! [definition]
	! <case1> IDCT (excluding scale)
	! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *
	! cos(pij1(k1+1/2)/n1) *
	! cos(pij2(k2+1/2)/n2),
	! 0<=k1<n1, 0<=k2<n2
	! <case2> DCT
	! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *
	! cos(pi(j1+1/2)k1/n1) *
	! cos(pi(j2+1/2)k2/n2),
	! 0<=k1<n1, 0<=k2<n2
	! [usage]
	! <case1>
	! ip(0) = 0 ! first time only
	! call ddct2d(n1max, n1, n2, 1, a, t, ip, w)
	! <case2>
	! ip(0) = 0 ! first time only
	! call ddct2d(n1max, n1, n2, -1, a, t, ip, w)
	! [parameters]
	! n1max :row size of the 2D array (integer)
	! n1 :data length (integer)
	! n1 >= 2, n1 = power of 2
	! n2 :data length (integer)
	! n2 >= 2, n2 = power of 2
	! a(0:n1-1,0:n2-1)
	! :input/output data (real*8)
	! output data
	! a(k1,k2) = C(k1,k2), 0<=k1<n1, 0<=k2<n2
	! t(0:4*n2-1)
	! :work area (real*8)
	! length of t >= 4*n2
	! ip(0:*):work area for bit reversal (integer)
	! length of ip >= 2+sqrt(n)
	! (n = max(n1/2, n2/2))
	! ip(0),ip(1) are pointers of the cos/sin table.
	! w(0:) :cos/sin table (real8)
	! length of w >= max(n13/2, n23/2)
	! w(),ip() are initialized if ip(0) = 0.
	! [remark]
	! Inverse of
	! call ddct2d(n1max, n1, n2, -1, a, t, ip, w)
	! is
	! do j1 = 0, n1 - 1
	! a(j1, 0) = a(j1, 0) * 0.5d0
	! end do
	! do j2 = 0, n2 - 1
	! a(0, j2) = a(0, j2) * 0.5d0
	! end do
	! call ddct2d(n1max, n1, n2, 1, a, t, ip, w)
	! do j2 = 0, n2 - 1
	! do j1 = 0, n1 - 1
	! a(j1, j2) = a(j1, j2) * (4.0d0 / n1 / n2)
	! end do
	! end do
	! .
	!
	!
	! -------- DST (Discrete Sine Transform) / Inverse of DST --------
	! [definition]
	! <case1> IDST (excluding scale)
	! S(k1,k2) = sum_j1=1^n1 sum_j2=1^n2 A(j1,j2) *
	! sin(pij1(k1+1/2)/n1) *
	! sin(pij2(k2+1/2)/n2),
	! 0<=k1<n1, 0<=k2<n2
	! <case2> DST
	! S(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) *
	! sin(pi(j1+1/2)k1/n1) *
	! sin(pi(j2+1/2)k2/n2),
	! 0<k1<=n1, 0<k2<=n2
	! [usage]
	! <case1>
	! ip(0) = 0 ! first time only
	! call ddst2d(n1max, n1, n2, 1, a, t, ip, w)
	! <case2>
	! ip(0) = 0 ! first time only
	! call ddst2d(n1max, n1, n2, -1, a, t, ip, w)
	! [parameters]
	! n1max :row size of the 2D array (integer)
	! n1 :data length (integer)
	! n1 >= 2, n1 = power of 2
	! n2 :data length (integer)
	! n2 >= 2, n2 = power of 2
	! a(0:n1-1,0:n2-1)
	! :input/output data (real*8)
	! <case1>
	! input data
	! a(j1,j2) = A(j1,j2), 0<j1<n1, 0<j2<n2,
	! a(j1,0) = A(j1,n2), 0<j1<n1,
	! a(0,j2) = A(n1,j2), 0<j2<n2,
	! a(0,0) = A(n1,n2)
	! (i.e. A(j1,j2) = a(mod(j1,n1),mod(j2,n2)))
	! output data
	! a(k1,k2) = S(k1,k2), 0<=k1<n1, 0<=k2<n2
	! <case2>
	! output data
	! a(k1,k2) = S(k1,k2), 0<k1<n1, 0<k2<n2,
	! a(k1,0) = S(k1,n2), 0<k1<n1,
	! a(0,k2) = S(n1,k2), 0<k2<n2,
	! a(0,0) = S(n1,n2)
	! (i.e. S(k1,k2) = a(mod(k1,n1),mod(k2,n2)))
	! t(0:4*n2-1)
	! :work area (real*8)
	! length of t >= 4*n2
	! ip(0:*):work area for bit reversal (integer)
	! length of ip >= 2+sqrt(n)
	! (n = max(n1/2, n2/2))
	! ip(0),ip(1) are pointers of the cos/sin table.
	! w(0:) :cos/sin table (real8)
	! length of w >= max(n13/2, n23/2)
	! w(),ip() are initialized if ip(0) = 0.
	! [remark]
	! Inverse of
	! call ddst2d(n1max, n1, n2, -1, a, t, ip, w)
	! is
	! do j1 = 0, n1 - 1
	! a(j1, 0) = a(j1, 0) * 0.5d0
	! end do
	! do j2 = 0, n2 - 1
	! a(0, j2) = a(0, j2) * 0.5d0
	! end do
	! call ddst2d(n1max, n1, n2, 1, a, t, ip, w)
	! do j2 = 0, n2 - 1
	! do j1 = 0, n1 - 1
	! a(j1, j2) = a(j1, j2) * (4.0d0 / n1 / n2)
	! end do
	! end do
	! .
	!
	!
	subroutine cdft2d(n1max, n1, n2, isgn, a, t, ip, w)
	integer n1max, n1, n2, isgn, ip(0 : *), n, j
	real8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 n2 - 1),
	& w(0 : *)
	n = max(n1, 2 * n2)
	if (n .gt. 4 * ip(0)) then
	call makewt(n / 4, ip, w)
	end if
	do j = 0, n2 - 1
	call cdft(n1, isgn, a(0, j), ip, w)
	end do
	call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)
	end
	!
	subroutine rdft2d(n1max, n1, n2, isgn, a, t, ip, w)
	integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j
	real8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 n2 - 1),
	& w(0 : *)
	n = max(n1, 2 * n2)
	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 rdft2d_sub(n1max, n1, n2, isgn, a)
	call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)
	end if
	do j = 0, n2 - 1
	call rdft(n1, isgn, a(0, j), ip, w)
	end do
	if (isgn .ge. 0) then
	call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)
	call rdft2d_sub(n1max, n1, n2, isgn, a)
	end if
	end
	!
	subroutine rdft2dsort(n1max, n1, n2, isgn, a)
	integer n1max, n1, n2, isgn, n2h, j
	real*8 a(0 : n1max - 1, 0 : n2 - 1), x, y
	n2h = n2 / 2
	if (isgn .lt. 0) then
	do j = n2h + 1, n2 - 1
	a(0, j) = a(n1 + 1, j)
	a(1, j) = a(n1, j)
	end do
	a(1, 0) = a(n1, 0)
	a(1, n2h) = a(n1, n2h)
	else
	do j = n2h + 1, n2 - 1
	y = a(0, j)
	x = a(1, j)
	a(n1, j) = x
	a(n1 + 1, j) = y
	a(n1, n2 - j) = x
	a(n1 + 1, n2 - j) = -y
	a(0, j) = a(0, n2 - j)
	a(1, j) = -a(1, n2 - j)
	end do
	a(n1, 0) = a(1, 0)
	a(n1 + 1, 0) = 0
	a(1, 0) = 0
	a(n1, n2h) = a(1, n2h)
	a(n1 + 1, n2h) = 0
	a(1, n2h) = 0
	end if
	end
	!
	subroutine ddct2d(n1max, n1, n2, isgn, a, t, ip, w)
	integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j
	real8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 n2 - 1),
	& w(0 : *)
	n = max(n1, n2)
	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
	do j = 0, n2 - 1
	call ddct(n1, isgn, a(0, j), ip, w)
	end do
	call ddxt2d_sub(n1max, n1, n2, 0, isgn, a, t, ip, w)
	end
	!
	subroutine ddst2d(n1max, n1, n2, isgn, a, t, ip, w)
	integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j
	real8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 n2 - 1),
	& w(0 : *)
	n = max(n1, n2)
	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
	do j = 0, n2 - 1
	call ddst(n1, isgn, a(0, j), ip, w)
	end do
	call ddxt2d_sub(n1max, n1, n2, 1, isgn, a, t, ip, w)
	end
	!
	! -------- child routines --------
	!
	subroutine cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w)
	integer n1max, n1, n2, isgn, ip(0 : *), i, j
	real8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 n2 - 1),
	& w(0 : *)
	if (n1 .gt. 4) then
	do i = 0, n1 - 8, 8
	do j = 0, n2 - 1
	t(2 * j) = a(i, j)
	t(2 * j + 1) = a(i + 1, j)
	t(2 * n2 + 2 * j) = a(i + 2, j)
	t(2 * n2 + 2 * j + 1) = a(i + 3, j)
	t(4 * n2 + 2 * j) = a(i + 4, j)
	t(4 * n2 + 2 * j + 1) = a(i + 5, j)
	t(6 * n2 + 2 * j) = a(i + 6, j)
	t(6 * n2 + 2 * j + 1) = a(i + 7, j)
	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) = t(2 * j)
	a(i + 1, j) = t(2 * j + 1)
	a(i + 2, j) = t(2 * n2 + 2 * j)
	a(i + 3, j) = t(2 * n2 + 2 * j + 1)
	a(i + 4, j) = t(4 * n2 + 2 * j)
	a(i + 5, j) = t(4 * n2 + 2 * j + 1)
	a(i + 6, j) = t(6 * n2 + 2 * j)
	a(i + 7, j) = 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)
	t(2 * j + 1) = a(1, j)
	t(2 * n2 + 2 * j) = a(2, j)
	t(2 * n2 + 2 * j + 1) = a(3, j)
	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) = t(2 * j)
	a(1, j) = t(2 * j + 1)
	a(2, j) = t(2 * n2 + 2 * j)
	a(3, j) = 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)
	t(2 * j + 1) = a(1, j)
	end do
	call cdft(2 * n2, isgn, t, ip, w)
	do j = 0, n2 - 1
	a(0, j) = t(2 * j)
	a(1, j) = t(2 * j + 1)
	end do
	end if
	end
	!
	subroutine rdft2d_sub(n1max, n1, n2, isgn, a)
	integer n1max, n1, n2, isgn, n2h, i, j
	real*8 a(0 : n1max - 1, 0 : n2 - 1), xi
	n2h = n2 / 2
	if (isgn .lt. 0) then
	do i = 1, n2h - 1
	j = n2 - i
	xi = a(0, i) - a(0, j)
	a(0, i) = a(0, i) + a(0, j)
	a(0, j) = xi
	xi = a(1, j) - a(1, i)
	a(1, i) = a(1, i) + a(1, j)
	a(1, j) = xi
	end do
	else
	do i = 1, n2h - 1
	j = n2 - i
	a(0, j) = 0.5d0 * (a(0, i) - a(0, j))
	a(0, i) = a(0, i) - a(0, j)
	a(1, j) = 0.5d0 * (a(1, i) + a(1, j))
	a(1, i) = a(1, i) - a(1, j)
	end do
	end if
	end
	!
	subroutine ddxt2d_sub(n1max, n1, n2, ics, isgn, a, t,
	& ip, w)
	integer n1max, n1, n2, ics, isgn, ip(0 : *), i, j
	real8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 n2 - 1),
	& w(0 : *)
	if (n1 .gt. 2) then
	do i = 0, n1 - 4, 4
	do j = 0, n2 - 1
	t(j) = a(i, j)
	t(n2 + j) = a(i + 1, j)
	t(2 * n2 + j) = a(i + 2, j)
	t(3 * n2 + j) = a(i + 3, j)
	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) = t(j)
	a(i + 1, j) = t(n2 + j)
	a(i + 2, j) = t(2 * n2 + j)
	a(i + 3, j) = t(3 * n2 + j)
	end do
	end do
	else if (n1 .eq. 2) then
	do j = 0, n2 - 1
	t(j) = a(0, j)
	t(n2 + j) = a(1, j)
	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) = t(j)
	a(1, j) = t(n2 + j)
	end do
	end if
	end
	!