! BFFFT (Blazingly Fast Fourier Transform) Mathematical tool Copyright (C) 2004 Nicolas Riesz.
!
! This program is free software; you can redistribute it and/or modify it under the terms of the GNU General
! Public License as published by the Free Software Foundation; either version 2 of the License, or (at 
! you option) any later version. 
!
! BFFFT is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
! for more details.
!
! You should have received a copy of the GNU General Public License
! along with BFFFT.  If not, write to the Free Software Foundation,
! 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

MODULE BFFFT

!
! Real-to-complex 2-D Vectorialized FFT
!
! Sample program to use VFFT2 :
!
! program main
! USE VFFT2
! integer     :: s1(N,N), s2(N,N), s3(N,N), s4(N,N)
! real*8      :: w(N+15)
! real*8      :: a(N,N), b(N,N)
!
! call VRFFTINIT(N,w,s1,s2,s3,s4)
!
! ! 2D-FFT of a
! call FFT2(N,a,b,WSAVE,s1,s2,s3,s4)
! 
! ! 2D-IFFT of a
! call IFFT2(N,a,b,WSAVE,s1,s2,s3,s4)
!
! end program
! 
! FORMAT of <a> after "FFT2" (2D Vectorialized Fast Fourier Transform)
!
!         1           2           3        4       5     ...       N
!1     R(0,0)       R(1,0)     I(1,0)    R(2,0)  I(2,0)  ...   R(N/2,0)
!2     R(0,1)        R11     -I(N-1,1)    R21  -I(N-2,1) ...   R(N/2,1)
!3     I(0,1)        I11      R(N-1,1)    I21   R(N-2,1) ...   I(N/2,1)
!4     R(0,2)        R12     -I(N-1,2)    ...     ...    ...     ...
!5     I(0,2)        I12      R(N-1,2)    ...     ...    ... 
!6      ...          ...         ...      ...     ...   ...      ...
!7      ...          ...         ...      ...     ...   ...      ...
!.      ...          ...         ...      ...     ...   ...      ...
!N-1  I(0,N/2-1)     ...         ...      ...     ...   ...  I(N/2,N/2-1)
!N     R(0,N/2)    R(1,N/2)    I(1,N/2)   ...     ...   ...   R(N/2,N/2)

!  where R(i,j) + sqrt(-1)*I(i,j) = Fourier coefficient on i, j.
! 
! Other coefficients can be deduced by symmetric properties, given that :
! 
!                 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), 
!
! for 0<k1<n1, 0<k2<n2.
!

implicit REAL*8 (a-h), REAL*8 (o-z)
implicit INTEGER(i-n)

integer, parameter :: N=256 ! SIZE OF THE INPUT
integer :: is1(N,N), is2(N,N), is3(N,N), is4(N,N) !flag-vectorized
real    :: dureeFourier2, dureeFourier3

contains

      subroutine FFT2(N,a,b,WSAVE)
      use fastTranspo
      implicit none
      ! 2D vectorizable FFT.
      ! a is the array to be Fourierized
      ! b is a working array
 
      integer :: N, i, j
      real*8  :: a(N,N), b(N,N)
      real*8  :: WSAVE(N+15)
      
      real :: etime, elapsed(2), temps

      call VRFFTF(N, N, a, b, N, WSAVE)

      temps = etime(elapsed)
      call fastTranspose(b, a, N)
      dureeFourier2 = dureeFourier2 + etime(elapsed)-temps

      call VRFFTF(N, N, b, a, N, WSAVE)

!
! matrix A is now of the form :
!
!         1           2           3        4       5    ...       N
!1     R(0,0)       R(1,0)     I(1,0)    R(2,0)  I(2,0) ...    R(N/2,0)
!2     R(0,1)        X11         Y11      X21     Y21   ...    R(N/2,1)
!3     I(0,1)        X12         Y12      X22     Y22   ...    I(N/2,1)
!4     R(0,2)        X13         Y13      X23     Y23   ...      ...
!5     I(0,2)        X14         Y14      X24     Y24   ...      ...
!6      ...          ...         ...      ...     ...   ...      ...
!7      ...          ...         ...      ...     ...   ...      ...
!.      ...          ...         ...      ...     ...   ...      ...
!N-1  I(0,N/2-1)  Y(N/2-1,1)  Y(N/2-1,1)  ...     ...   ...  I(N/2,N/2-1)
!N     R(0,N/2)    R(1,N/2)    I(1,N/2)   ...     ...   ...   R(N/2,N/2)

! with
!      R(k1  , k2) = X(k1, 2*k2-1) - Y(k1, 2*k2  )
!      I(k1  , k2) = X(k1, 2*k2  ) + Y(k1, 2*k2-1)
!      R(N-k1, k2) = X(k1, 2*k2-1) + Y(k1, 2*k2  )
!      I(N-k1, k2) = X(k1, 2*k2  ) - Y(k1, 2*k2-1)
!      0 < k1 < N/2
!      0 < k2 < N/2

! For each interior (2,2) bloc, we do the folowing transform :
!
! | a  b | _\  | a-d  b-c |
! | c  d |  /  | c+b  d+a |
!
      temps = etime(elapsed)

      a(1,:) = b(1,:)
      a(N,:) = b(N,:)
      a(:,1) = b(:,1)
      a(:,N) = b(:,N)

      do j=2,N-1
         do i=2,N-1
               a(i,j) = b(i,j) + is1(i, j) * b(i+1, j+1) &
                               + is2(i, j) * b(i-1, j-1) &
                               + is3(i, j) * b(i-1, j+1) &
                               + is4(i, j) * b(i+1, j-1)
         end do
      end do

      dureeFourier3 = dureeFourier3 + etime(elapsed)-temps

      end subroutine FFT2

      subroutine IFFT2(N,a,b,WSAVE)
      use fastTranspo
      implicit none
      ! 2D vectorizable FFT. 
      ! a is the array to be inverse-Fourierized 
      ! b is a working array

      integer :: N, i, j
      real*8 :: a(N,N), b(N,N)
      real*8 :: WSAVE(N+15)
      
      real :: etime, elapsed(2), temps

      !      Format for a
      !
      !      R(k1  , k2) = X(k1, 2*k2-1) - Y(k1, 2*k2  )
      !      I(k1  , k2) = X(k1, 2*k2  ) + Y(k1, 2*k2-1)
      !      R(N-k1, k2) = X(k1, 2*k2-1) + Y(k1, 2*k2  )
      !      I(N-k1, k2) = X(k1, 2*k2  ) - Y(k1, 2*k2-1)
      !      0 < k1 < N/2
      !      0 < k2 < N/2

      temps = etime(elapsed)

      b(1,:) = a(1,:)
      b(N,:) = a(N,:)
      b(:,1) = a(:,1)
      b(:,N) = a(:,N)

      do j=2,N-1
        do i=2,N-1
             b(i,j) = 0.5d0*( a(i,j) - is1(i, j) * a(i+1, j+1) &
                                     - is2(i, j) * a(i-1, j-1) &
                                     - is3(i, j) * a(i-1, j+1) &
                                     - is4(i, j) * a(i+1, j-1) )
        end do
      end do

      dureeFourier3 = dureeFourier3 + etime(elapsed)-temps

      call VRFFTB(N, N, b, a, N, WSAVE)

      temps = etime(elapsed)
      call fastTranspose(a, b, N)        
      dureeFourier2 = dureeFourier2 + etime(elapsed)-temps

      call VRFFTB(N, N, a, b, N, WSAVE)

      end subroutine IFFT2

      subroutine VRFFTINIT(N,WSAVE,numUnit,verbose)
      use fastTranspo
      implicit none
      
      integer, intent(IN)  :: N, numUnit, verbose
      real*8 , intent(OUT) :: WSAVE(N+15)
      integer :: i1,i2

      call initFastTranspose(N,numUnit,verbose)
      
      dureeFourier2 = 0.0
      dureeFourier3 = 0.0
      is1=0
      is2=0
      is3=0
      is4=0
      do i1=1,N/2-1
       do i2=1,N/2-1
          is1(2*i1  , 2*i2  ) = -1 ! pour faire -d    [cf. VFFT2usuel]
          is2(2*i1+1, 2*i2+1) =  1 ! pour faire +a 
          is3(2*i1+1, 2*i2  ) =  1 ! pour faire +b  
          is4(2*i1  , 2*i2+1) = -1 ! pour faire -c 
       end do
      end do

      call  VRFFTI(N,WSAVE)
   
      end subroutine VRFFTINIT

      SUBROUTINE VCOSQB(M,N,X,XT,MDIMX,WSAVE)
!***BEGIN PROLOGUE  VCOSQB
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, ODD WAVE
!             NUMBERS, MULTIPLE SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Normalized inverse of VCOSQF.
!***DESCRIPTION
!
!  Subroutine VCOSQB computes the backward fast Fourier cosine transform
!  of M quarter wave sequences.  That is, cosine series representations
!  with only odd wave numbers.  The transform is defined below at output
!  parameter X.
!
!  The array WSAVE which is used by subroutine VCOSQB must be
!  initialized by calling subroutine VCOSQI(N,WSAVE).
!
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequences to be transformed.  The method
!          is most efficient when N is a product of small primes.
!
!  X       an array of size at least X(MDIMX,N) which contains the
!          the sequences to be transformed.  The sequences are stored
!          in the ROWS of X.  Thus, the Jth sequence is stored in
!          X(J,I), I=1,..,N.
!
!  XT      a work array of size at least XT(MDIMX,N).
!
!  MDIMX   the first dimension of the array X exactly as it appears in
!          the calling program.
!
!  WSAVE   a work array which must be dimensioned at least 2*N+15
!          in the program that calls VCOSQB.  The WSAVE array must be
!          initialized by calling subroutine VCOSQI(N,WSAVE), and a
!          different WSAVE array must be used for each different
!          value of N.  This initialization does not have to be
!          repeated so long as N remains unchanged.
!
!  Output Parameters
!
!  X       For I=1,...,N and J=1,...,M
!
!               X(I)= the sum from K=1 to K=N of
!
!                 4*X(K)*COS((2*K-1)*(I-1)*PI/(2*N)) /SQRT(4*N)
!
!  WSAVE   contains initialization calculations which must not
!          be destroyed between calls of VCOSQF or VCOSQB.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VCOSQF followed immediately by a call of
!           of VCOSQB will return the original sequences X.  Thus,
!           VCOSQB is the correctly normalized inverse VCOSQF.
!
!  -----------------------------------------------------------------
!
!  VCOSQB is a straightforward extension of the subprogram COSQB to
!  handle M simultaneous sequences.  COSQB was orginally developed
!  by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  VCOSQB
      DIMENSION       X(MDIMX,*), XT(MDIMX,*), WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VCOSQB
      IF (M .LE. 0)  GO TO 900
      IF (N .GT. 2)  GO TO 300
      IF (N .EQ. 2)  GO TO 200
      GO TO 900
!
!  CASE  N = 2
!
  200 CONTINUE
      SCALE = 2.0d0*SQRT(0.50d0)
      DO 210 J=1,M
         X1 = SCALE*(X(J,1)+X(J,2))
         X(J,2) = X(J,1)-X(J,2)
         X(J,1) = X1
  210 CONTINUE
      GO TO 900
!
!  CASE N .GT. 2
!
!     ... PREPROCESSING
!
  300 CONTINUE
      NS2 = (N+1)/2
      NP2 = N+2
      DO 310 I=3,N,2
         DO 310 J=1,M
            XIM1 = X(J,I-1)+X(J,I)
            X(J,I) = X(J,I)-X(J,I-1)
            X(J,I-1) = XIM1
  310 CONTINUE
      DO 320 J=1,M
         X(J,1) = X(J,1)+X(J,1)
  320 CONTINUE
      MODN = MOD(N,2)
      IF (MODN .EQ. 0) THEN
         DO 330 J=1,M
            X(J,N) = X(J,N)+X(J,N)
  330    CONTINUE
      ENDIF
!
!     ... REAL, PERIODIC TRANSFORM
!
      CALL VRFFTB (M,N,X,XT,MDIMX,WSAVE(N+1))
!
!     ... POSTPROCESSING
!
      DO 340 K=2,NS2
         KC = NP2-K
         DO 340 J=1,M
            XT(J,K) = WSAVE(K-1)*X(J,KC)+WSAVE(KC-1)*X(J,K)
            XT(J,KC) = WSAVE(K-1)*X(J,K)-WSAVE(KC-1)*X(J,KC)
  340 CONTINUE
      IF (MODN .EQ. 0) THEN
         DO 350 J=1,M
            X(J,NS2+1) = WSAVE(NS2)*(X(J,NS2+1)+X(J,NS2+1))
  350    CONTINUE
      ENDIF
      DO 360 K=2,NS2
         KC = NP2-K
         DO 360 J=1,M
            X(J,K) = XT(J,K)+XT(J,KC)
            X(J,KC) = XT(J,K)-XT(J,KC)
  360 CONTINUE
      DO 370 J=1,M
         X(J,1) = X(J,1)+X(J,1)
  370 CONTINUE
!
!     ... NORMALIZATION
!
      SCALE = 0.5d0
      DO 380 I=1,N
         DO 380 J=1,M
            X(J,I) = SCALE*X(J,I)
  380 CONTINUE
!
!  EXIT
!
  900 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VCOSQF(M,N,X,XT,MDIMX,WSAVE)
!***BEGIN PROLOGUE  VCOSQF
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, ODD WAVE
!             NUMBERS, MULTIPLE SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Forward cosine transform, odd wave numbers, M sequences.
!***DESCRIPTION
!
!  Subroutine VCOSQF computes the forward fast Fourier cosine transform
!  of M quarter wave sequences.  That is, cosine series representations
!  with only odd wave numbers.  The transform is defined below at output
!  parameter X.
!
!  The array WSAVE which is used by subroutine VCOSQF must be
!  initialized by calling subroutine VCOSQI(N,WSAVE).
!
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequences to be transformed.  The method
!          is most efficient when N is a product of small primes.
!
!  X       an array of size at least X(MDIMX,N) which contains the
!          the sequences to be transformed.  The sequences are stored
!          in the ROWS of X.  Thus, the Jth sequence is stored in
!          X(J,I), I=1,..,N.
!
!  XT      a work array of size at least XT(MDIMX,N).
!
!  MDIMX   the first dimension of the array X exactly as it appears in
!          the calling program.
!
!  WSAVE   a work array which must be dimensioned at least 2*N+15
!          in the program that calls VCOSQF.  The WSAVE array must be
!          initialized by calling subroutine VCOSQI(N,WSAVE), and a
!          different WSAVE array must be used for each different
!          value of N.  This initialization does not have to be
!          repeated so long as N remains unchanged.
!
!  Output Parameters
!
!  X       For I=1,...,N and J=1,...,M
!
!               X(I) = ( X(1) + the sum from K=2 to K=N of
!
!                  2*X(K)*COS((2*I-1)*(K-1)*PI/(2*N)) )/SQRT(4*N)
!
!  WSAVE   contains initialization calculations which must not
!          be destroyed between calls of VCOSQF or VCOSQB.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VCOSQF followed immediately by a call of
!           of VCOSQB will return the original sequences X.  Thus,
!           VCOSQB is the correctly normalized inverse VCOSQF.
!
!  -----------------------------------------------------------------
!
!  VCOSQF is a straightforward extension of the subprogram COSQF to
!  handle M simultaneous sequences.  COSQF was orginally developed
!  by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  VCOSQF
      DIMENSION       X(MDIMX,*), XT(MDIMX,*), WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VCOSQF
      IF (M .LE. 0)  GO TO 900
      IF (N .GT. 2)  GO TO 300
      IF (N .LT. 2)  GO TO 900
!
!  CASE  N = 2
!
  200 CONTINUE
      SQRT2 = SQRT(2.0d0)
      SCALE = 0.50d0/SQRT2
      DO 210 J=1,M
         TSQX = SQRT2*X(J,2)
         X(J,2) = SCALE*(X(J,1)-TSQX)
         X(J,1) = SCALE*(X(J,1)+TSQX)
  210 CONTINUE
      GO TO 900
!
!  CASE N .GT. 2
!
  300 CONTINUE
!
!     ... PREPROCESSING
!
      NS2 = (N+1)/2
      NP2 = N+2
      DO 310 K=2,NS2
         KC = NP2-K
         DO 310 J=1,M
            XT(J,K) = X(J,K)+X(J,KC)
            XT(J,KC) = X(J,K)-X(J,KC)
  310 CONTINUE
      MODN = MOD(N,2)
      IF (MODN .EQ. 0) THEN
         DO 320 J=1,M
            XT(J,NS2+1) = X(J,NS2+1)+X(J,NS2+1)
  320    CONTINUE
      ENDIF
      DO 330 K=2,NS2
         KC = NP2-K
         DO 330 J=1,M
            X(J,K) = WSAVE(K-1)*XT(J,KC)+WSAVE(KC-1)*XT(J,K)
            X(J,KC) = WSAVE(K-1)*XT(J,K)-WSAVE(KC-1)*XT(J,KC)
  330 CONTINUE
      IF (MODN .EQ. 0) THEN
         DO 340 J=1,M
            X(J,NS2+1) = WSAVE(NS2)*XT(J,NS2+1)
  340    CONTINUE
      ENDIF
!
!     ... REAL, PERIODIC TRANSFORM
!
      CALL VRFFTF (M,N,X,XT,MDIMX,WSAVE(N+1))
!
!     ... POSTPROCESSING
!
      DO 350 I=3,N,2
         DO 350 J=1,M
            XIM1 = X(J,I-1)-X(J,I)
            X(J,I) = X(J,I-1)+X(J,I)
            X(J,I-1) = XIM1
  350 CONTINUE
!
!     ... NORMALIZATION
!
      SCALE = 0.5d0
      DO 360 I=1,N
         DO 360 J=1,M
            X(J,I) = SCALE*X(J,I)
  360 CONTINUE
!
!  EXIT
!
  900 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VCOSQI(N,WSAVE)
!***BEGIN PROLOGUE  VCOSQI
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, ODD WAVE
!             NUMBERS, MULTIPLE SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Initialize for VCOSQF and VCOSQB.
!***DESCRIPTION
!
!  Subroutine VCOSQI initializes the array WSAVE which is used in
!  both VCOSQF and VCOSQB.  The prime factorization of N together with
!  a tabulation of the trigonometric functions are computed and
!  stored in WSAVE.
!
!  Input Parameter
!
!  N       the length of the array to be transformed.  The method
!          is most efficient when N is a product of small primes.
!
!  Output Parameter
!
!  WSAVE   a work array which must be dimensioned at least 2*N+15.
!          The same work array can be used for both VCOSQF and VCOSQB
!          as long as N remains unchanged.  Different WSAVE arrays
!          are required for different values of N.  The contents of
!          WSAVE must not be changed between calls of VCOSQF or VCOSQB.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  VRFFTI
!***END PROLOGUE  VCOSQI
      DIMENSION       WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VCOSQI
      PIH = 2.0d0*ATAN(1.0d0)
      DT = PIH/REAL(N)
      FK = 0.0d0
      DO 101 K=1,N
         FK = FK+1.0d0
         WSAVE(K) = COS(FK*DT)
  101 CONTINUE
      CALL VRFFTI (N,WSAVE(N+1))
      RETURN
      END SUBROUTINE
      SUBROUTINE VCOST(M,N,X,XT,MDIMX,WSAVE)
!***BEGIN PROLOGUE  VCOST
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, MULTIPLE
!             SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Cosine transform of one or more real, even sequences.
!***DESCRIPTION
!
!  Subroutine VCOST computes the discrete Fourier cosine transform
!  of M even sequences X(J,I), J=1,...,M.  The transform is defined
!  below at output parameter X.
!
!  The array WSAVE which is used by subroutine VCOST must be
!  initialized by calling subroutine VCOSTI(N,WSAVE).
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequence to be transformed.  N must be
!          greater than 1.  The method is most efficient when N-1 is
!          is a product of small primes.
!
!  X       an array of size at least X(MDIMX,N) which contains the
!          the sequences to be transformed.  The sequences are stored
!          in the ROWS of X.  Thus, the Jth sequence is stored in
!          X(J,I), I=1,..,N.
!
!  XT      a work array of size at least XT(MDIMX,N-1).
!
!  MDIMX   the first dimension of the array X exactly as it appears in
!          the calling program.
!
!  WSAVE   a work array which must be dimensioned at least 3*N+15
!          in the program that calls VCOST.  The WSAVE array must be
!          initialized by calling subroutine VCOSTI(N,WSAVE), and a
!          different WSAVE array must be used for each different
!          value of N.  This initialization does not have to be
!          repeated so long as N remains unchanged.  Thus subsequent
!          transforms can be obtained faster than the first.
!
!  Output Parameters
!
!  X       For I=1,...,N and J=1,...,M
!
!             X(J,I) = ( X(J,1)+(-1)**(I-1)*X(J,N)
!
!               + the sum from K=2 to K=N-1
!
!                 2*X(J,K)*COS((K-1)*(I-1)*PI/(N-1)) )/SQRT(2*(N-1))
!
!  WSAVE   contains initialization calculations which must not be
!          destroyed between calls of VCOST.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VCOST followed immediately by another call
!           of VCOST will return the original sequences X.  Thus,
!           VCOST is the correctly normalized inverse of itself.
!
!  -----------------------------------------------------------------
!
!  VCOST is a straightforward extension of the subprogram COST to
!  handle M simultaneous sequences.  The scaling of the sequences
!  computed by VCOST is different than that of COST.  COST was
!  orginally developed by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  VRFFTF
!***END PROLOGUE  VCOST
      DIMENSION       X(MDIMX,*), XT(MDIMX,*), WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VCOST
      IF (M .LE. 0)  GO TO 900
      IF (N .LE. 1)  GO TO 900
      IF (N .GT. 3)  GO TO 400
      IF (N .EQ. 3)  GO TO 300
!
!  CASE  N = 2
!
  200 CONTINUE
      SCALE = SQRT(0.50d0)
      DO 210 J=1,M
         X1H = SCALE*(X(J,1)+X(J,2))
         X(J,2) = SCALE*(X(J,1)-X(J,2))
         X(J,1) = X1H
  210 CONTINUE
      GO TO 900
!
!  CASE  N = 3
!
  300 CONTINUE
      SCALE = 0.50d0
      DO 310 J=1,M
         X1P3 = X(J,1)+X(J,3)
         TX2 = X(J,2)+X(J,2)
         X(J,2) = SCALE*(X(J,1)-X(J,3))
         X(J,1) = SCALE*(X1P3+TX2)
         X(J,3) = SCALE*(X1P3-TX2)
  310 CONTINUE
      GO TO 900
!
!  CASE  N .GT. 3
!
!     ... PREPROCESSING
!
  400 CONTINUE
      NM1 = N-1
      NP1 = N+1
      NS2 = N/2
      DO 410 J=1,M
         XT(J,1) = X(J,1)-X(J,N)
         X(J,1) = X(J,1)+X(J,N)
  410 CONTINUE
      DO 420 K=2,NS2
         KC = NP1-K
         DO 420 J=1,M
            T1 = X(J,K)+X(J,KC)
            T2 = X(J,K)-X(J,KC)
            XT(J,1) = XT(J,1)+WSAVE(KC)*T2
            T2 = WSAVE(K)*T2
            X(J,K) = T1-T2
            X(J,KC) = T1+T2
  420 CONTINUE
      MODN = MOD(N,2)
      IF (MODN .NE. 0) THEN
         DO 430 J=1,M
            X(J,NS2+1) = X(J,NS2+1)+X(J,NS2+1)
  430    CONTINUE
      ENDIF
      DO 435 J=1,M
         X(J,N) = XT(J,1)
  435 CONTINUE
!
!     ... REAL*8 PERIODIC TRANSFORM
!
      CALL VRFFTF (M,NM1,X,XT,MDIMX,WSAVE(NP1))
!
!     ... POSTPROCESSING
!
      FACTOR = 1.0d0/SQRT(REAL(NM1))
      DO 440 J=1,M
         XT(J,1) = X(J,2)
         X(J,2) = FACTOR*X(J,N)
  440 CONTINUE
      DO 450 I=4,N,2
         DO 450 J=1,M
            XI = X(J,I)
            X(J,I) = X(J,I-2)-X(J,I-1)
            X(J,I-1) = XT(J,1)
            XT(J,1) = XI
  450 CONTINUE
      IF (MODN .NE. 0) THEN
         DO 460 J=1,M
            X(J,N) = XT(J,1)
  460    CONTINUE
      ENDIF
!
!     ... NORMALIZATION
!
      SCALE = SQRT(0.5d0)
      DO 490 I=1,N
         DO 490 J=1,M
            X(J,I) = SCALE*X(J,I)
  490 CONTINUE
!
!  EXIT
!
  900 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VCOSTI(N,WSAVE)
!***BEGIN PROLOGUE  VCOSTI
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, MULTIPLE
!             SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Initialize for VCOST.
!***DESCRIPTION
!
!  Subroutine VCOSTI initializes the array WSAVE which is used in
!  subroutine VCOST.  The prime factorization of N together with
!  a tabulation of the trigonometric functions are computed and
!  stored in WSAVE.
!
!  Input Parameter
!
!  N       the length of the sequence to be transformed.  The method
!          is most efficient when N-1 is a product of small primes.
!
!  Output Parameter
!
!  WSAVE   a work array which must be dimensioned at least 3*N+15.
!          Different WSAVE arrays are required for different values
!          of N.  The contents of WSAVE must not be changed between
!          calls of VCOST.
!
!  -----------------------------------------------------------------
!
!  VCOSTI is a straightforward extension of the subprogram COSTI to
!  handle M simultaneous sequences.  COSTI was orginally developed
!  by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  VRFFTI
!***END PROLOGUE  VCOSTI
      DIMENSION       WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VCOSTI
      PI = 4.0d0*ATAN(1.0d0)
      IF (N .LE. 3) RETURN
      NM1 = N-1
      NP1 = N+1
      NS2 = N/2
      DT = PI/REAL(NM1)
      FK = 0.0d0
      DO 101 K=2,NS2
         KC = NP1-K
         FK = FK+1.0d0
         WSAVE(K) = 2.0d0*SIN(FK*DT)
         WSAVE(KC) = 2.0d0*COS(FK*DT)
  101 CONTINUE
      CALL VRFFTI (NM1,WSAVE(N+1))
      RETURN
      END SUBROUTINE
      SUBROUTINE VRADB2 (MP,IDO,L1,CC,CH,MDIMC,WA1)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION  CC(MDIMC,IDO,2,L1)    ,CH(MDIMC,IDO,L1,2),&
                     WA1(IDO)
      DO 101 K=1,L1
          DO 1001 M=1,MP
         CH(M,1,K,1) = CC(M,1,1,K)+CC(M,IDO,2,K)
         CH(M,1,K,2) = CC(M,1,1,K)-CC(M,IDO,2,K)
 1001     CONTINUE
  101 CONTINUE
      IF (IDO-2) 107,105,102
  102 IDP2 = IDO+2
      DO 104 K=1,L1
         DO 103 I=3,IDO,2
            IC = IDP2-I
               DO 1002 M=1,MP
            CH(M,I-1,K,1) = CC(M,I-1,1,K)+CC(M,IC-1,2,K)
            CH(M,I,K,1) = CC(M,I,1,K)-CC(M,IC,2,K)
            CH(M,I-1,K,2) = WA1(I-2)*(CC(M,I-1,1,K)-CC(M,IC-1,2,K))&
           -WA1(I-1)*(CC(M,I,1,K)+CC(M,IC,2,K))
            CH(M,I,K,2) = WA1(I-2)*(CC(M,I,1,K)+CC(M,IC,2,K))+WA1(I-1)&
           *(CC(M,I-1,1,K)-CC(M,IC-1,2,K))
 1002          CONTINUE
  103    CONTINUE
  104 CONTINUE
      IF (MOD(IDO,2) .EQ. 1) RETURN
  105 DO 106 K=1,L1
          DO 1003 M=1,MP
         CH(M,IDO,K,1) = CC(M,IDO,1,K)+CC(M,IDO,1,K)
         CH(M,IDO,K,2) = -(CC(M,1,2,K)+CC(M,1,2,K))
 1003     CONTINUE
  106 CONTINUE
  107 RETURN
      END SUBROUTINE
      SUBROUTINE VRADB3 (MP,IDO,L1,CC,CH,MDIMC,WA1,WA2)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION  CC(MDIMC,IDO,3,L1)    ,CH(MDIMC,IDO,L1,3),&
                     WA1(IDO)   ,WA2(IDO)
      ARG=2.0d0*3.141592653589793238462643/3.
      TAUR=COS(ARG)
      TAUI=SIN(ARG)
      DO 101 K=1,L1
          DO 1001 M=1,MP
         CH(M,1,K,1) = CC(M,1,1,K)+2.0d0*CC(M,IDO,2,K)
         CH(M,1,K,2) = CC(M,1,1,K)+(2.0d0*TAUR)*CC(M,IDO,2,K)&
        -(2.0d0*TAUI)*CC(M,1,3,K)
         CH(M,1,K,3) = CC(M,1,1,K)+(2.0d0*TAUR)*CC(M,IDO,2,K)&
        +2.0d0*TAUI*CC(M,1,3,K)
 1001     CONTINUE
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      DO 103 K=1,L1
         DO 102 I=3,IDO,2
            IC = IDP2-I
               DO 1002 M=1,MP
            CH(M,I-1,K,1) = CC(M,I-1,1,K)+(CC(M,I-1,3,K)+CC(M,IC-1,2,K))
            CH(M,I,K,1) = CC(M,I,1,K)+(CC(M,I,3,K)-CC(M,IC,2,K))
            CH(M,I-1,K,2) = WA1(I-2)*&
      ((CC(M,I-1,1,K)+TAUR*(CC(M,I-1,3,K)+CC(M,IC-1,2,K)))-&
       (TAUI*(CC(M,I,3,K)+CC(M,IC,2,K))))&
                         -WA1(I-1)*&
       ((CC(M,I,1,K)+TAUR*(CC(M,I,3,K)-CC(M,IC,2,K)))+&
       (TAUI*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))))
            CH(M,I,K,2) = WA1(I-2)*&
       ((CC(M,I,1,K)+TAUR*(CC(M,I,3,K)-CC(M,IC,2,K)))+&
       (TAUI*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))))&
                        +WA1(I-1)*&
       ((CC(M,I-1,1,K)+TAUR*(CC(M,I-1,3,K)+CC(M,IC-1,2,K)))-&
       (TAUI*(CC(M,I,3,K)+CC(M,IC,2,K))))
              CH(M,I-1,K,3) = WA2(I-2)*&
       ((CC(M,I-1,1,K)+TAUR*(CC(M,I-1,3,K)+CC(M,IC-1,2,K)))+&
       (TAUI*(CC(M,I,3,K)+CC(M,IC,2,K))))&
                           -WA2(I-1)*&
       ((CC(M,I,1,K)+TAUR*(CC(M,I,3,K)-CC(M,IC,2,K)))-&
       (TAUI*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))))
            CH(M,I,K,3) = WA2(I-2)*&
      ((CC(M,I,1,K)+TAUR*(CC(M,I,3,K)-CC(M,IC,2,K)))-&
       (TAUI*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))))&
                       +WA2(I-1)*&
       ((CC(M,I-1,1,K)+TAUR*(CC(M,I-1,3,K)+CC(M,IC-1,2,K)))+&
       (TAUI*(CC(M,I,3,K)+CC(M,IC,2,K))))
 1002          CONTINUE
  102    CONTINUE
  103 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VRADB4 (MP,IDO,L1,CC,CH,MDIMC,WA1,WA2,WA3)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION  CC(MDIMC,IDO,4,L1)  ,CH(MDIMC,IDO,L1,4)    ,&
                     WA1(IDO)  ,WA2(IDO)  ,WA3(IDO)
      SQRT2=SQRT(2.0d0)
      DO 101 K=1,L1
          DO 1001 M=1,MP
         CH(M,1,K,3) = (CC(M,1,1,K)+CC(M,IDO,4,K))&
        -(CC(M,IDO,2,K)+CC(M,IDO,2,K))
         CH(M,1,K,1) = (CC(M,1,1,K)+CC(M,IDO,4,K))&
        +(CC(M,IDO,2,K)+CC(M,IDO,2,K))
         CH(M,1,K,4) = (CC(M,1,1,K)-CC(M,IDO,4,K))&
        +(CC(M,1,3,K)+CC(M,1,3,K))
         CH(M,1,K,2) = (CC(M,1,1,K)-CC(M,IDO,4,K))&
        -(CC(M,1,3,K)+CC(M,1,3,K))
 1001     CONTINUE
  101 CONTINUE
      IF (IDO-2) 107,105,102
  102 IDP2 = IDO+2
      DO 104 K=1,L1
         DO 103 I=3,IDO,2
            IC = IDP2-I
               DO 1002 M=1,MP
            CH(M,I-1,K,1) = (CC(M,I-1,1,K)+CC(M,IC-1,4,K))&
           +(CC(M,I-1,3,K)+CC(M,IC-1,2,K))
            CH(M,I,K,1) = (CC(M,I,1,K)-CC(M,IC,4,K))&
           +(CC(M,I,3,K)-CC(M,IC,2,K))
            CH(M,I-1,K,2)=WA1(I-2)*((CC(M,I-1,1,K)-CC(M,IC-1,4,K))&
           -(CC(M,I,3,K)+CC(M,IC,2,K)))-WA1(I-1)&
           *((CC(M,I,1,K)+CC(M,IC,4,K))+(CC(M,I-1,3,K)-CC(M,IC-1,2,K)))
            CH(M,I,K,2)=WA1(I-2)*((CC(M,I,1,K)+CC(M,IC,4,K))&
           +(CC(M,I-1,3,K)-CC(M,IC-1,2,K)))+WA1(I-1)&
           *((CC(M,I-1,1,K)-CC(M,IC-1,4,K))-(CC(M,I,3,K)+CC(M,IC,2,K)))
            CH(M,I-1,K,3)=WA2(I-2)*((CC(M,I-1,1,K)+CC(M,IC-1,4,K))&
           -(CC(M,I-1,3,K)+CC(M,IC-1,2,K)))-WA2(I-1)&
           *((CC(M,I,1,K)-CC(M,IC,4,K))-(CC(M,I,3,K)-CC(M,IC,2,K)))
            CH(M,I,K,3)=WA2(I-2)*((CC(M,I,1,K)-CC(M,IC,4,K))&
           -(CC(M,I,3,K)-CC(M,IC,2,K)))+WA2(I-1)&
           *((CC(M,I-1,1,K)+CC(M,IC-1,4,K))-(CC(M,I-1,3,K)&
           +CC(M,IC-1,2,K)))
            CH(M,I-1,K,4)=WA3(I-2)*((CC(M,I-1,1,K)-CC(M,IC-1,4,K))&
           +(CC(M,I,3,K)+CC(M,IC,2,K)))-WA3(I-1)&
          *((CC(M,I,1,K)+CC(M,IC,4,K))-(CC(M,I-1,3,K)-CC(M,IC-1,2,K)))
            CH(M,I,K,4)=WA3(I-2)*((CC(M,I,1,K)+CC(M,IC,4,K))&
           -(CC(M,I-1,3,K)-CC(M,IC-1,2,K)))+WA3(I-1)&
           *((CC(M,I-1,1,K)-CC(M,IC-1,4,K))+(CC(M,I,3,K)+CC(M,IC,2,K)))
 1002          CONTINUE
  103    CONTINUE
  104 CONTINUE
      IF (MOD(IDO,2) .EQ. 1) RETURN
  105 CONTINUE
      DO 106 K=1,L1
               DO 1003 M=1,MP
         CH(M,IDO,K,1) = (CC(M,IDO,1,K)+CC(M,IDO,3,K))&
        +(CC(M,IDO,1,K)+CC(M,IDO,3,K))
         CH(M,IDO,K,2) = SQRT2*((CC(M,IDO,1,K)-CC(M,IDO,3,K))&
        -(CC(M,1,2,K)+CC(M,1,4,K)))
         CH(M,IDO,K,3) = (CC(M,1,4,K)-CC(M,1,2,K))&
        +(CC(M,1,4,K)-CC(M,1,2,K))
         CH(M,IDO,K,4) = -SQRT2*((CC(M,IDO,1,K)-CC(M,IDO,3,K))&
        +(CC(M,1,2,K)+CC(M,1,4,K)))
 1003          CONTINUE
  106 CONTINUE
  107 RETURN
      END SUBROUTINE
      SUBROUTINE VRADB5 (MP,IDO,L1,CC,CH,MDIMC,WA1,WA2,WA3,WA4)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION  CC(MDIMC,IDO,5,L1)    ,CH(MDIMC,IDO,L1,5),&
                  WA1(IDO)     ,WA2(IDO)     ,WA3(IDO)     ,WA4(IDO)
      ARG=2.0d0*3.141592653589793238462643/5.
      TR11=COS(ARG)
      TI11=SIN(ARG)
      TR12=COS(2.0d0*ARG)
      TI12=SIN(2.0d0*ARG)
      DO 101 K=1,L1
      DO 1001 M=1,MP
         CH(M,1,K,1) = CC(M,1,1,K)+2.0d0*CC(M,IDO,2,K)+2.0d0*CC(M,IDO,4,K)
         CH(M,1,K,2) = (CC(M,1,1,K)+TR11*2.0d0*CC(M,IDO,2,K)&
        +TR12*2.0d0*CC(M,IDO,4,K))-(TI11*2.0d0*CC(M,1,3,K)&
        +TI12*2.0d0*CC(M,1,5,K))
         CH(M,1,K,3) = (CC(M,1,1,K)+TR12*2.0d0*CC(M,IDO,2,K)&
        +TR11*2.0d0*CC(M,IDO,4,K))-(TI12*2.0d0*CC(M,1,3,K)&
        -TI11*2.0d0*CC(M,1,5,K))
         CH(M,1,K,4) = (CC(M,1,1,K)+TR12*2.0d0*CC(M,IDO,2,K)&
        +TR11*2.0d0*CC(M,IDO,4,K))+(TI12*2.0d0*CC(M,1,3,K)&
        -TI11*2.0d0*CC(M,1,5,K))
         CH(M,1,K,5) = (CC(M,1,1,K)+TR11*2.0d0*CC(M,IDO,2,K)&
        +TR12*2.0d0*CC(M,IDO,4,K))+(TI11*2.0d0*CC(M,1,3,K)&
        +TI12*2.0d0*CC(M,1,5,K))
 1001          CONTINUE
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      DO 103 K=1,L1
         DO 102 I=3,IDO,2
            IC = IDP2-I
      DO 1002 M=1,MP
            CH(M,I-1,K,1) = CC(M,I-1,1,K)+(CC(M,I-1,3,K)+CC(M,IC-1,2,K))&
           +(CC(M,I-1,5,K)+CC(M,IC-1,4,K))
            CH(M,I,K,1) = CC(M,I,1,K)+(CC(M,I,3,K)-CC(M,IC,2,K))&
           +(CC(M,I,5,K)-CC(M,IC,4,K))
            CH(M,I-1,K,2) = WA1(I-2)*((CC(M,I-1,1,K)+TR11*&
           (CC(M,I-1,3,K)+CC(M,IC-1,2,K))+TR12&
           *(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))-(TI11*(CC(M,I,3,K)&
           +CC(M,IC,2,K))+TI12*(CC(M,I,5,K)+CC(M,IC,4,K))))&
           -WA1(I-1)*((CC(M,I,1,K)+TR11*(CC(M,I,3,K)-CC(M,IC,2,K))&
           +TR12*(CC(M,I,5,K)-CC(M,IC,4,K)))+(TI11*(CC(M,I-1,3,K)&
           -CC(M,IC-1,2,K))+TI12*(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))
            CH(M,I,K,2) = WA1(I-2)*((CC(M,I,1,K)+TR11*(CC(M,I,3,K)&
           -CC(M,IC,2,K))+TR12*(CC(M,I,5,K)-CC(M,IC,4,K)))&
           +(TI11*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))+TI12&
           *(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))+WA1(I-1)&
           *((CC(M,I-1,1,K)+TR11*(CC(M,I-1,3,K)&
           +CC(M,IC-1,2,K))+TR12*(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))&
           -(TI11*(CC(M,I,3,K)+CC(M,IC,2,K))+TI12&
           *(CC(M,I,5,K)+CC(M,IC,4,K))))
            CH(M,I-1,K,3) = WA2(I-2)&
           *((CC(M,I-1,1,K)+TR12*(CC(M,I-1,3,K)+CC(M,IC-1,2,K))&
           +TR11*(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))-(TI12*(CC(M,I,3,K)&
           +CC(M,IC,2,K))-TI11*(CC(M,I,5,K)+CC(M,IC,4,K))))&
          -WA2(I-1)&
          *((CC(M,I,1,K)+TR12*(CC(M,I,3,K)-&
           CC(M,IC,2,K))+TR11*(CC(M,I,5,K)-CC(M,IC,4,K)))&
           +(TI12*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))-TI11&
           *(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))
            CH(M,I,K,3) = WA2(I-2)&
          *((CC(M,I,1,K)+TR12*(CC(M,I,3,K)-&
           CC(M,IC,2,K))+TR11*(CC(M,I,5,K)-CC(M,IC,4,K)))&
           +(TI12*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))-TI11&
           *(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))&
           +WA2(I-1)&
           *((CC(M,I-1,1,K)+TR12*(CC(M,I-1,3,K)+CC(M,IC-1,2,K))&
           +TR11*(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))-(TI12*(CC(M,I,3,K)&
           +CC(M,IC,2,K))-TI11*(CC(M,I,5,K)+CC(M,IC,4,K))))
            CH(M,I-1,K,4) = WA3(I-2)&
           *((CC(M,I-1,1,K)+TR12*(CC(M,I-1,3,K)+CC(M,IC-1,2,K))&
           +TR11*(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))+(TI12*(CC(M,I,3,K)&
           +CC(M,IC,2,K))-TI11*(CC(M,I,5,K)+CC(M,IC,4,K))))&
           -WA3(I-1)&
          *((CC(M,I,1,K)+TR12*(CC(M,I,3,K)-&
           CC(M,IC,2,K))+TR11*(CC(M,I,5,K)-CC(M,IC,4,K)))&
           -(TI12*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))-TI11&
           *(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))
            CH(M,I,K,4) = WA3(I-2)&
          *((CC(M,I,1,K)+TR12*(CC(M,I,3,K)-&
           CC(M,IC,2,K))+TR11*(CC(M,I,5,K)-CC(M,IC,4,K)))&
           -(TI12*(CC(M,I-1,3,K)-CC(M,IC-1,2,K))-TI11&
           *(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))&
           +WA3(I-1)&
           *((CC(M,I-1,1,K)+TR12*(CC(M,I-1,3,K)+CC(M,IC-1,2,K))&
           +TR11*(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))+(TI12*(CC(M,I,3,K)&
           +CC(M,IC,2,K))-TI11*(CC(M,I,5,K)+CC(M,IC,4,K))))
            CH(M,I-1,K,5) = WA4(I-2)&
           *((CC(M,I-1,1,K)+TR11*(CC(M,I-1,3,K)+CC(M,IC-1,2,K))&
           +TR12*(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))+(TI11*(CC(M,I,3,K)&
           +CC(M,IC,2,K))+TI12*(CC(M,I,5,K)+CC(M,IC,4,K))))&
           -WA4(I-1)&
           *((CC(M,I,1,K)+TR11*(CC(M,I,3,K)-CC(M,IC,2,K))&
           +TR12*(CC(M,I,5,K)-CC(M,IC,4,K)))-(TI11*(CC(M,I-1,3,K)&
           -CC(M,IC-1,2,K))+TI12*(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))
            CH(M,I,K,5) = WA4(I-2)&
           *((CC(M,I,1,K)+TR11*(CC(M,I,3,K)-CC(M,IC,2,K))&
           +TR12*(CC(M,I,5,K)-CC(M,IC,4,K)))-(TI11*(CC(M,I-1,3,K)&
           -CC(M,IC-1,2,K))+TI12*(CC(M,I-1,5,K)-CC(M,IC-1,4,K))))&
           +WA4(I-1)&
           *((CC(M,I-1,1,K)+TR11*(CC(M,I-1,3,K)+CC(M,IC-1,2,K))&
           +TR12*(CC(M,I-1,5,K)+CC(M,IC-1,4,K)))+(TI11*(CC(M,I,3,K)&
           +CC(M,IC,2,K))+TI12*(CC(M,I,5,K)+CC(M,IC,4,K))))
 1002          CONTINUE
  102    CONTINUE
  103 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VRADBG (MP,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,&
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
                       MDIMC,WA)
      DIMENSION    CH(MDIMC,IDO,L1,IP)    ,CC(MDIMC,IDO,IP,L1) ,&
                C1(MDIMC,IDO,L1,IP)     ,C2(MDIMC,IDL1,IP),&
                      CH2(MDIMC,IDL1,IP)       ,WA(IDO)
      TPI=2.0d0*3.141592653589793238462643
      ARG = TPI/FLOAT(IP)
      DCP = COS(ARG)
      DSP = SIN(ARG)
      IDP2 = IDO+2
      NBD = (IDO-1)/2
      IPP2 = IP+2
      IPPH = (IP+1)/2
      IF (IDO .LT. L1) GO TO 103
      DO 102 K=1,L1
         DO 101 I=1,IDO
            DO 1001 M=1,MP
            CH(M,I,K,1) = CC(M,I,1,K)
 1001       CONTINUE
  101    CONTINUE
  102 CONTINUE
      GO TO 106
  103 DO 105 I=1,IDO
         DO 104 K=1,L1
            DO 1004 M=1,MP
            CH(M,I,K,1) = CC(M,I,1,K)
 1004       CONTINUE
  104    CONTINUE
  105 CONTINUE
  106 DO 108 J=2,IPPH
         JC = IPP2-J
         J2 = J+J
         DO 107 K=1,L1
            DO 1007 M=1,MP
            CH(M,1,K,J) = CC(M,IDO,J2-2,K)+CC(M,IDO,J2-2,K)
            CH(M,1,K,JC) = CC(M,1,J2-1,K)+CC(M,1,J2-1,K)
 1007       CONTINUE
  107    CONTINUE
  108 CONTINUE
      IF (IDO .EQ. 1) GO TO 116
      IF (NBD .LT. L1) GO TO 112
      DO 111 J=2,IPPH
         JC = IPP2-J
         DO 110 K=1,L1
            DO 109 I=3,IDO,2
               IC = IDP2-I
               DO 1009 M=1,MP
               CH(M,I-1,K,J) = CC(M,I-1,2*J-1,K)+CC(M,IC-1,2*J-2,K)
               CH(M,I-1,K,JC) = CC(M,I-1,2*J-1,K)-CC(M,IC-1,2*J-2,K)
               CH(M,I,K,J) = CC(M,I,2*J-1,K)-CC(M,IC,2*J-2,K)
               CH(M,I,K,JC) = CC(M,I,2*J-1,K)+CC(M,IC,2*J-2,K)
 1009          CONTINUE
  109       CONTINUE
  110    CONTINUE
  111 CONTINUE
      GO TO 116
  112 DO 115 J=2,IPPH
         JC = IPP2-J
         DO 114 I=3,IDO,2
            IC = IDP2-I
            DO 113 K=1,L1
               DO 1013 M=1,MP
               CH(M,I-1,K,J) = CC(M,I-1,2*J-1,K)+CC(M,IC-1,2*J-2,K)
               CH(M,I-1,K,JC) = CC(M,I-1,2*J-1,K)-CC(M,IC-1,2*J-2,K)
               CH(M,I,K,J) = CC(M,I,2*J-1,K)-CC(M,IC,2*J-2,K)
               CH(M,I,K,JC) = CC(M,I,2*J-1,K)+CC(M,IC,2*J-2,K)
 1013          CONTINUE
  113       CONTINUE
  114    CONTINUE
  115 CONTINUE
  116 AR1 = 1.0d0
      AI1 = 0.0d0
      DO 120 L=2,IPPH
         LC = IPP2-L
         AR1H = DCP*AR1-DSP*AI1
         AI1 = DCP*AI1+DSP*AR1
         AR1 = AR1H
         DO 117 IK=1,IDL1
            DO 1017 M=1,MP
            C2(M,IK,L) = CH2(M,IK,1)+AR1*CH2(M,IK,2)
            C2(M,IK,LC) = AI1*CH2(M,IK,IP)
 1017       CONTINUE
  117    CONTINUE
         DC2 = AR1
         DS2 = AI1
         AR2 = AR1
         AI2 = AI1
         DO 119 J=3,IPPH
            JC = IPP2-J
            AR2H = DC2*AR2-DS2*AI2
            AI2 = DC2*AI2+DS2*AR2
            AR2 = AR2H
            DO 118 IK=1,IDL1
               DO 1018 M=1,MP
               C2(M,IK,L) = C2(M,IK,L)+AR2*CH2(M,IK,J)
               C2(M,IK,LC) = C2(M,IK,LC)+AI2*CH2(M,IK,JC)
 1018          CONTINUE
  118       CONTINUE
  119    CONTINUE
  120 CONTINUE
      DO 122 J=2,IPPH
         DO 121 IK=1,IDL1
            DO 1021 M=1,MP
            CH2(M,IK,1) = CH2(M,IK,1)+CH2(M,IK,J)
 1021       CONTINUE
  121    CONTINUE
  122 CONTINUE
      DO 124 J=2,IPPH
         JC = IPP2-J
         DO 123 K=1,L1
            DO 1023 M=1,MP
            CH(M,1,K,J) = C1(M,1,K,J)-C1(M,1,K,JC)
            CH(M,1,K,JC) = C1(M,1,K,J)+C1(M,1,K,JC)
 1023       CONTINUE
  123    CONTINUE
  124 CONTINUE
      IF (IDO .EQ. 1) GO TO 132
      IF (NBD .LT. L1) GO TO 128
      DO 127 J=2,IPPH
         JC = IPP2-J
         DO 126 K=1,L1
            DO 125 I=3,IDO,2
               DO 1025 M=1,MP
               CH(M,I-1,K,J) = C1(M,I-1,K,J)-C1(M,I,K,JC)
               CH(M,I-1,K,JC) = C1(M,I-1,K,J)+C1(M,I,K,JC)
               CH(M,I,K,J) = C1(M,I,K,J)+C1(M,I-1,K,JC)
               CH(M,I,K,JC) = C1(M,I,K,J)-C1(M,I-1,K,JC)
 1025          CONTINUE
  125       CONTINUE
  126    CONTINUE
  127 CONTINUE
      GO TO 132
  128 DO 131 J=2,IPPH
         JC = IPP2-J
         DO 130 I=3,IDO,2
            DO 129 K=1,L1
               DO 1029 M=1,MP
               CH(M,I-1,K,J) = C1(M,I-1,K,J)-C1(M,I,K,JC)
               CH(M,I-1,K,JC) = C1(M,I-1,K,J)+C1(M,I,K,JC)
               CH(M,I,K,J) = C1(M,I,K,J)+C1(M,I-1,K,JC)
               CH(M,I,K,JC) = C1(M,I,K,J)-C1(M,I-1,K,JC)
 1029          CONTINUE
  129       CONTINUE
  130    CONTINUE
  131 CONTINUE
  132 CONTINUE
      IF (IDO .EQ. 1) RETURN
      DO 133 IK=1,IDL1
         DO 1033 M=1,MP
         C2(M,IK,1) = CH2(M,IK,1)
 1033    CONTINUE
  133 CONTINUE
      DO 135 J=2,IP
         DO 134 K=1,L1
            DO 1034 M=1,MP
            C1(M,1,K,J) = CH(M,1,K,J)
 1034       CONTINUE
  134    CONTINUE
  135 CONTINUE
      IF (NBD .GT. L1) GO TO 139
      IS = -IDO
      DO 138 J=2,IP
         IS = IS+IDO
         IDIJ = IS
         DO 137 I=3,IDO,2
            IDIJ = IDIJ+2
            DO 136 K=1,L1
               DO 1036 M=1,MP
               C1(M,I-1,K,J) = WA(IDIJ-1)*CH(M,I-1,K,J)-WA(IDIJ)*&
               CH(M,I,K,J)
               C1(M,I,K,J) = WA(IDIJ-1)*CH(M,I,K,J)+WA(IDIJ)*&
               CH(M,I-1,K,J)
 1036          CONTINUE
  136       CONTINUE
  137    CONTINUE
  138 CONTINUE
      GO TO 143
  139 IS = -IDO
      DO 142 J=2,IP
         IS = IS+IDO
         DO 141 K=1,L1
            IDIJ = IS
            DO 140 I=3,IDO,2
               IDIJ = IDIJ+2
               DO 1040 M=1,MP
               C1(M,I-1,K,J) = WA(IDIJ-1)*CH(M,I-1,K,J)-WA(IDIJ)*&
               CH(M,I,K,J)
               C1(M,I,K,J) = WA(IDIJ-1)*CH(M,I,K,J)+WA(IDIJ)*&
               CH(M,I-1,K,J)
 1040          CONTINUE
  140       CONTINUE
  141    CONTINUE
  142 CONTINUE
  143 RETURN
      END SUBROUTINE
      SUBROUTINE VRADF2 (MP,IDO,L1,CC,CH,MDIMC,WA1)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION   CH(MDIMC,IDO,2,L1)  ,CC(MDIMC,IDO,L1,2)     ,&
                     WA1(IDO)
      DO 101 K=1,L1
         DO 1001 M=1,MP
         CH(M,1,1,K) = CC(M,1,K,1)+CC(M,1,K,2)
         CH(M,IDO,2,K) = CC(M,1,K,1)-CC(M,1,K,2)
 1001    CONTINUE
  101 CONTINUE
      IF (IDO-2) 107,105,102
  102 IDP2 = IDO+2
      DO 104 K=1,L1
         DO 103 I=3,IDO,2
            IC = IDP2-I
            DO 1003 M=1,MP
            CH(M,I,1,K) = CC(M,I,K,1)+(WA1(I-2)*CC(M,I,K,2)-&
            WA1(I-1)*CC(M,I-1,K,2))
            CH(M,IC,2,K) = (WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*&
            CC(M,I-1,K,2))-CC(M,I,K,1)
            CH(M,I-1,1,K) = CC(M,I-1,K,1)+(WA1(I-2)*CC(M,I-1,K,2)+&
            WA1(I-1)*CC(M,I,K,2))
            CH(M,IC-1,2,K) = CC(M,I-1,K,1)-(WA1(I-2)*CC(M,I-1,K,2)+&
            WA1(I-1)*CC(M,I,K,2))
 1003       CONTINUE
  103    CONTINUE
  104 CONTINUE
      IF (MOD(IDO,2) .EQ. 1) RETURN
  105 DO 106 K=1,L1
         DO 1006 M=1,MP
         CH(M,1,2,K) = -CC(M,IDO,K,2)
         CH(M,IDO,1,K) = CC(M,IDO,K,1)
 1006    CONTINUE
  106 CONTINUE
  107 RETURN
      END SUBROUTINE
      SUBROUTINE VRADF3 (MP,IDO,L1,CC,CH,MDIMC,WA1,WA2)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION   CH(MDIMC,IDO,3,L1)  ,CC(MDIMC,IDO,L1,3)     ,&
                     WA1(IDO)     ,WA2(IDO)
      ARG=2.0d0*3.141592653589793238462643/3.
      TAUR=COS(ARG)
      TAUI=SIN(ARG)
      DO 101 K=1,L1
         DO 1001 M=1,MP
         CH(M,1,1,K) = CC(M,1,K,1)+(CC(M,1,K,2)+CC(M,1,K,3))
         CH(M,1,3,K) = TAUI*(CC(M,1,K,3)-CC(M,1,K,2))
         CH(M,IDO,2,K) = CC(M,1,K,1)+TAUR*&
           (CC(M,1,K,2)+CC(M,1,K,3))
 1001    CONTINUE
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      DO 103 K=1,L1
         DO 102 I=3,IDO,2
            IC = IDP2-I
            DO 1002 M=1,MP
            CH(M,I-1,1,K) = CC(M,I-1,K,1)+((WA1(I-2)*CC(M,I-1,K,2)+&
            WA1(I-1)*CC(M,I,K,2))+(WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3)))
            CH(M,I,1,K) = CC(M,I,K,1)+((WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*&
            CC(M,I-1,K,2))+(WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3)))
            CH(M,I-1,3,K) = (CC(M,I-1,K,1)+TAUR*((WA1(I-2)*&
            CC(M,I-1,K,2)+WA1(I-1)*CC(M,I,K,2))+(WA2(I-2)*&
            CC(M,I-1,K,3)+WA2(I-1)*CC(M,I,K,3))))+(TAUI*((WA1(I-2)*&
            CC(M,I,K,2)-WA1(I-1)*CC(M,I-1,K,2))-(WA2(I-2)*&
            CC(M,I,K,3)-WA2(I-1)*CC(M,I-1,K,3))))
            CH(M,IC-1,2,K) = (CC(M,I-1,K,1)+TAUR*((WA1(I-2)*&
            CC(M,I-1,K,2)+WA1(I-1)*CC(M,I,K,2))+(WA2(I-2)*&
            CC(M,I-1,K,3)+WA2(I-1)*CC(M,I,K,3))))-(TAUI*((WA1(I-2)*&
            CC(M,I,K,2)-WA1(I-1)*CC(M,I-1,K,2))-(WA2(I-2)*&
            CC(M,I,K,3)-WA2(I-1)*CC(M,I-1,K,3))))
            CH(M,I,3,K) = (CC(M,I,K,1)+TAUR*((WA1(I-2)*CC(M,I,K,2)-&
            WA1(I-1)*CC(M,I-1,K,2))+(WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3))))+(TAUI*((WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2))))
            CH(M,IC,2,K) = (TAUI*((WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2))))-(CC(M,I,K,1)+TAUR*((WA1(I-2)*CC(M,I,K,2)-&
            WA1(I-1)*CC(M,I-1,K,2))+(WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3))))
 1002       CONTINUE
  102    CONTINUE
  103 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VRADF4 (MP,IDO,L1,CC,CH,MDIMC,WA1,WA2,WA3)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION    CC(MDIMC,IDO,L1,4)   ,CH(MDIMC,IDO,4,L1)     ,&
                     WA1(IDO)     ,WA2(IDO)     ,WA3(IDO)
      HSQT2=SQRT(2.0d0)/2.0d0
      DO 101 K=1,L1
         DO 1001 M=1,MP
         CH(M,1,1,K) = (CC(M,1,K,2)+CC(M,1,K,4))&
           +(CC(M,1,K,1)+CC(M,1,K,3))
         CH(M,IDO,4,K) = (CC(M,1,K,1)+CC(M,1,K,3))&
           -(CC(M,1,K,2)+CC(M,1,K,4))
         CH(M,IDO,2,K) = CC(M,1,K,1)-CC(M,1,K,3)
         CH(M,1,3,K) = CC(M,1,K,4)-CC(M,1,K,2)
 1001    CONTINUE
  101 CONTINUE
      IF (IDO-2) 107,105,102
  102 IDP2 = IDO+2
      DO 104 K=1,L1
         DO 103 I=3,IDO,2
            IC = IDP2-I
            DO 1003 M=1,MP
            CH(M,I-1,1,K) = ((WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2))+(WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*&
            CC(M,I,K,4)))+(CC(M,I-1,K,1)+(WA2(I-2)*CC(M,I-1,K,3)+&
            WA2(I-1)*CC(M,I,K,3)))
            CH(M,IC-1,4,K) = (CC(M,I-1,K,1)+(WA2(I-2)*CC(M,I-1,K,3)+&
            WA2(I-1)*CC(M,I,K,3)))-((WA1(I-2)*CC(M,I-1,K,2)+&
            WA1(I-1)*CC(M,I,K,2))+(WA3(I-2)*CC(M,I-1,K,4)+&
            WA3(I-1)*CC(M,I,K,4)))
            CH(M,I,1,K) = ((WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*&
            CC(M,I-1,K,2))+(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4)))+(CC(M,I,K,1)+(WA2(I-2)*CC(M,I,K,3)-&
            WA2(I-1)*CC(M,I-1,K,3)))
            CH(M,IC,4,K) = ((WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*&
            CC(M,I-1,K,2))+(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4)))-(CC(M,I,K,1)+(WA2(I-2)*CC(M,I,K,3)-&
            WA2(I-1)*CC(M,I-1,K,3)))
            CH(M,I-1,3,K) = ((WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*&
            CC(M,I-1,K,2))-(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4)))+(CC(M,I-1,K,1)-(WA2(I-2)*CC(M,I-1,K,3)+&
            WA2(I-1)*CC(M,I,K,3)))
            CH(M,IC-1,2,K) = (CC(M,I-1,K,1)-(WA2(I-2)*CC(M,I-1,K,3)+&
            WA2(I-1)*CC(M,I,K,3)))-((WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*&
            CC(M,I-1,K,2))-(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4)))
            CH(M,I,3,K) = ((WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*&
            CC(M,I,K,4))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2)))+(CC(M,I,K,1)-(WA2(I-2)*CC(M,I,K,3)-&
            WA2(I-1)*CC(M,I-1,K,3)))
            CH(M,IC,2,K) = ((WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*&
            CC(M,I,K,4))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2)))-(CC(M,I,K,1)-(WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3)))
 1003       CONTINUE
  103    CONTINUE
  104 CONTINUE
      IF (MOD(IDO,2) .EQ. 1) RETURN
  105 CONTINUE
      DO 106 K=1,L1
         DO 1006 M=1,MP
            CH(M,IDO,1,K) = (HSQT2*(CC(M,IDO,K,2)-CC(M,IDO,K,4)))+&
            CC(M,IDO,K,1)
            CH(M,IDO,3,K) = CC(M,IDO,K,1)-(HSQT2*(CC(M,IDO,K,2)-&
            CC(M,IDO,K,4)))
            CH(M,1,2,K) = (-HSQT2*(CC(M,IDO,K,2)+CC(M,IDO,K,4)))-&
            CC(M,IDO,K,3)
            CH(M,1,4,K) = (-HSQT2*(CC(M,IDO,K,2)+CC(M,IDO,K,4)))+&
            CC(M,IDO,K,3)
 1006    CONTINUE
  106 CONTINUE
  107 RETURN
      END SUBROUTINE
      SUBROUTINE VRADF5 (MP,IDO,L1,CC,CH,MDIMC,WA1,WA2,WA3,WA4)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION  CC(MDIMC,IDO,L1,5)    ,CH(MDIMC,IDO,5,L1)     ,&
                WA1(IDO)     ,WA2(IDO)     ,WA3(IDO)     ,WA4(IDO)
      ARG=2.0d0*3.141592653589793238462643/5.
      TR11=COS(ARG)
      TI11=SIN(ARG)
      TR12=COS(2.0d0*ARG)
      TI12=SIN(2.0d0*ARG)
      DO 101 K=1,L1
         DO 1001 M=1,MP
         CH(M,1,1,K) = CC(M,1,K,1)+(CC(M,1,K,5)+CC(M,1,K,2))+&
         (CC(M,1,K,4)+CC(M,1,K,3))
         CH(M,IDO,2,K) = CC(M,1,K,1)+TR11*(CC(M,1,K,5)+CC(M,1,K,2))+&
         TR12*(CC(M,1,K,4)+CC(M,1,K,3))
         CH(M,1,3,K) = TI11*(CC(M,1,K,5)-CC(M,1,K,2))+TI12*&
         (CC(M,1,K,4)-CC(M,1,K,3))
         CH(M,IDO,4,K) = CC(M,1,K,1)+TR12*(CC(M,1,K,5)+CC(M,1,K,2))+&
         TR11*(CC(M,1,K,4)+CC(M,1,K,3))
         CH(M,1,5,K) = TI12*(CC(M,1,K,5)-CC(M,1,K,2))-TI11*&
         (CC(M,1,K,4)-CC(M,1,K,3))
 1001    CONTINUE
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      DO 103 K=1,L1
         DO 102 I=3,IDO,2
            IC = IDP2-I
            DO 1002 M=1,MP
            CH(M,I-1,1,K) = CC(M,I-1,K,1)+((WA1(I-2)*CC(M,I-1,K,2)+&
            WA1(I-1)*CC(M,I,K,2))+(WA4(I-2)*CC(M,I-1,K,5)+WA4(I-1)*&
            CC(M,I,K,5)))+((WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3))+(WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*CC(M,I,K,4)))
            CH(M,I,1,K) = CC(M,I,K,1)+((WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*&
            CC(M,I-1,K,2))+(WA4(I-2)*CC(M,I,K,5)-WA4(I-1)*&
            CC(M,I-1,K,5)))+((WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3))+(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4)))
            CH(M,I-1,3,K) = CC(M,I-1,K,1)+TR11*&
           ( WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*CC(M,I,K,2)&
            +WA4(I-2)*CC(M,I-1,K,5)+WA4(I-1)*CC(M,I,K,5))+TR12*&
           ( WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*CC(M,I,K,3)&
            +WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*CC(M,I,K,4))+TI11*&
           ( WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*CC(M,I-1,K,2)&
            -(WA4(I-2)*CC(M,I,K,5)-WA4(I-1)*CC(M,I-1,K,5)))+TI12*&
           ( WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*CC(M,I-1,K,3)&
            -(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*CC(M,I-1,K,4)))
            CH(M,IC-1,2,K) = CC(M,I-1,K,1)+TR11*&
           ( WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*CC(M,I,K,2)&
            +WA4(I-2)*CC(M,I-1,K,5)+WA4(I-1)*CC(M,I,K,5))+TR12*&
          ( WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*CC(M,I,K,3)&
           +WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*CC(M,I,K,4))-(TI11*&
           ( WA1(I-2)*CC(M,I,K,2)-WA1(I-1)*CC(M,I-1,K,2)&
              -(WA4(I-2)*CC(M,I,K,5)-WA4(I-1)*CC(M,I-1,K,5)))+TI12*&
           ( WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*CC(M,I-1,K,3)&
            -(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*CC(M,I-1,K,4))))
            CH(M,I,3,K) = (CC(M,I,K,1)+TR11*((WA1(I-2)*CC(M,I,K,2)-&
            WA1(I-1)*CC(M,I-1,K,2))+(WA4(I-2)*CC(M,I,K,5)-WA4(I-1)*&
            CC(M,I-1,K,5)))+TR12*((WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3))+(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4))))+(TI11*((WA4(I-2)*CC(M,I-1,K,5)+&
            WA4(I-1)*CC(M,I,K,5))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2)))+TI12*((WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*&
            CC(M,I,K,4))-(WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3))))
            CH(M,IC,2,K) = (TI11*((WA4(I-2)*CC(M,I-1,K,5)+WA4(I-1)*&
            CC(M,I,K,5))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2)))+TI12*((WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*&
            CC(M,I,K,4))-(WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3))))-(CC(M,I,K,1)+TR11*((WA1(I-2)*CC(M,I,K,2)-&
            WA1(I-1)*CC(M,I-1,K,2))+(WA4(I-2)*CC(M,I,K,5)-WA4(I-1)*&
            CC(M,I-1,K,5)))+TR12*((WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3))+(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4))))
            CH(M,I-1,5,K) = (CC(M,I-1,K,1)+TR12*((WA1(I-2)*&
            CC(M,I-1,K,2)+WA1(I-1)*CC(M,I,K,2))+(WA4(I-2)*&
            CC(M,I-1,K,5)+WA4(I-1)*CC(M,I,K,5)))+TR11*((WA2(I-2)*&
            CC(M,I-1,K,3)+WA2(I-1)*CC(M,I,K,3))+(WA3(I-2)*&
            CC(M,I-1,K,4)+WA3(I-1)*CC(M,I,K,4))))+(TI12*((WA1(I-2)*&
            CC(M,I,K,2)-WA1(I-1)*CC(M,I-1,K,2))-(WA4(I-2)*CC(M,I,K,5)-&
            WA4(I-1)*CC(M,I-1,K,5)))-TI11*((WA2(I-2)*CC(M,I,K,3)-&
            WA2(I-1)*CC(M,I-1,K,3))-(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4))))
            CH(M,IC-1,4,K) = (CC(M,I-1,K,1)+TR12*((WA1(I-2)*&
            CC(M,I-1,K,2)+WA1(I-1)*CC(M,I,K,2))+(WA4(I-2)*&
            CC(M,I-1,K,5)+WA4(I-1)*CC(M,I,K,5)))+TR11*((WA2(I-2)*&
            CC(M,I-1,K,3)+WA2(I-1)*CC(M,I,K,3))+(WA3(I-2)*&
            CC(M,I-1,K,4)+WA3(I-1)*CC(M,I,K,4))))-(TI12*((WA1(I-2)*&
            CC(M,I,K,2)-WA1(I-1)*CC(M,I-1,K,2))-(WA4(I-2)*CC(M,I,K,5)-&
            WA4(I-1)*CC(M,I-1,K,5)))-TI11*((WA2(I-2)*CC(M,I,K,3)-&
            WA2(I-1)*CC(M,I-1,K,3))-(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4))))
            CH(M,I,5,K) = (CC(M,I,K,1)+TR12*((WA1(I-2)*CC(M,I,K,2)-&
            WA1(I-1)*CC(M,I-1,K,2))+(WA4(I-2)*CC(M,I,K,5)-WA4(I-1)*&
            CC(M,I-1,K,5)))+TR11*((WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3))+(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4))))+(TI12*((WA4(I-2)*CC(M,I-1,K,5)+&
            WA4(I-1)*CC(M,I,K,5))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2)))-TI11*((WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*&
            CC(M,I,K,4))-(WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3))))
            CH(M,IC,4,K) = (TI12*((WA4(I-2)*CC(M,I-1,K,5)+WA4(I-1)*&
            CC(M,I,K,5))-(WA1(I-2)*CC(M,I-1,K,2)+WA1(I-1)*&
            CC(M,I,K,2)))-TI11*((WA3(I-2)*CC(M,I-1,K,4)+WA3(I-1)*&
            CC(M,I,K,4))-(WA2(I-2)*CC(M,I-1,K,3)+WA2(I-1)*&
            CC(M,I,K,3))))-(CC(M,I,K,1)+TR12*((WA1(I-2)*CC(M,I,K,2)-&
            WA1(I-1)*CC(M,I-1,K,2))+(WA4(I-2)*CC(M,I,K,5)-WA4(I-1)*&
            CC(M,I-1,K,5)))+TR11*((WA2(I-2)*CC(M,I,K,3)-WA2(I-1)*&
            CC(M,I-1,K,3))+(WA3(I-2)*CC(M,I,K,4)-WA3(I-1)*&
            CC(M,I-1,K,4))))
 1002       CONTINUE
  102    CONTINUE
  103 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VRADFG (MP,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,MDIMC,WA)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION     CH(MDIMC,IDO,L1,IP)   ,CC(MDIMC,IDO,IP,L1)  ,&
                 C1(MDIMC,IDO,L1,IP)    ,C2(MDIMC,IDL1,IP),&
                      CH2(MDIMC,IDL1,IP)           ,WA(IDO)
      TPI=2.0d0*3.141592653589793238462643
      ARG = TPI/FLOAT(IP)
      DCP = COS(ARG)
      DSP = SIN(ARG)
      IPPH = (IP+1)/2
      IPP2 = IP+2
      IDP2 = IDO+2
      NBD = (IDO-1)/2
      IF (IDO .EQ. 1) GO TO 119
      DO 101 IK=1,IDL1
         DO 1001 M=1,MP
         CH2(M,IK,1) = C2(M,IK,1)
 1001    CONTINUE
  101 CONTINUE
      DO 103 J=2,IP
         DO 102 K=1,L1
            DO 1002 M=1,MP
            CH(M,1,K,J) = C1(M,1,K,J)
 1002       CONTINUE
  102    CONTINUE
  103 CONTINUE
      IF (NBD .GT. L1) GO TO 107
      IS = -IDO
      DO 106 J=2,IP
         IS = IS+IDO
         IDIJ = IS
         DO 105 I=3,IDO,2
            IDIJ = IDIJ+2
            DO 104 K=1,L1
               DO 1004 M=1,MP
               CH(M,I-1,K,J) = WA(IDIJ-1)*C1(M,I-1,K,J)+WA(IDIJ)&
                *C1(M,I,K,J)
               CH(M,I,K,J) = WA(IDIJ-1)*C1(M,I,K,J)-WA(IDIJ)&
                *C1(M,I-1,K,J)
 1004          CONTINUE
  104       CONTINUE
  105    CONTINUE
  106 CONTINUE
      GO TO 111
  107 IS = -IDO
      DO 110 J=2,IP
         IS = IS+IDO
         DO 109 K=1,L1
            IDIJ = IS
            DO 108 I=3,IDO,2
               IDIJ = IDIJ+2
               DO 1008 M=1,MP
               CH(M,I-1,K,J) = WA(IDIJ-1)*C1(M,I-1,K,J)+WA(IDIJ)&
                *C1(M,I,K,J)
               CH(M,I,K,J) = WA(IDIJ-1)*C1(M,I,K,J)-WA(IDIJ)&
                *C1(M,I-1,K,J)
 1008          CONTINUE
  108       CONTINUE
  109    CONTINUE
  110 CONTINUE
  111 IF (NBD .LT. L1) GO TO 115
      DO 114 J=2,IPPH
         JC = IPP2-J
         DO 113 K=1,L1
            DO 112 I=3,IDO,2
               DO 1012 M=1,MP
               C1(M,I-1,K,J) = CH(M,I-1,K,J)+CH(M,I-1,K,JC)
               C1(M,I-1,K,JC) = CH(M,I,K,J)-CH(M,I,K,JC)
               C1(M,I,K,J) = CH(M,I,K,J)+CH(M,I,K,JC)
               C1(M,I,K,JC) = CH(M,I-1,K,JC)-CH(M,I-1,K,J)
 1012          CONTINUE
  112       CONTINUE
  113    CONTINUE
  114 CONTINUE
      GO TO 121
  115 DO 118 J=2,IPPH
         JC = IPP2-J
         DO 117 I=3,IDO,2
            DO 116 K=1,L1
               DO 1016 M=1,MP
               C1(M,I-1,K,J) = CH(M,I-1,K,J)+CH(M,I-1,K,JC)
               C1(M,I-1,K,JC) = CH(M,I,K,J)-CH(M,I,K,JC)
               C1(M,I,K,J) = CH(M,I,K,J)+CH(M,I,K,JC)
               C1(M,I,K,JC) = CH(M,I-1,K,JC)-CH(M,I-1,K,J)
 1016          CONTINUE
  116       CONTINUE
  117    CONTINUE
  118 CONTINUE
      GO TO 121
  119 DO 120 IK=1,IDL1
         DO 1020 M=1,MP
         C2(M,IK,1) = CH2(M,IK,1)
 1020    CONTINUE
  120 CONTINUE
  121 DO 123 J=2,IPPH
         JC = IPP2-J
         DO 122 K=1,L1
            DO 1022 M=1,MP
            C1(M,1,K,J) = CH(M,1,K,J)+CH(M,1,K,JC)
            C1(M,1,K,JC) = CH(M,1,K,JC)-CH(M,1,K,J)
 1022       CONTINUE
  122    CONTINUE
  123 CONTINUE
!
      AR1 = 1.0d0
      AI1 = 0.0d0
      DO 127 L=2,IPPH
         LC = IPP2-L
         AR1H = DCP*AR1-DSP*AI1
         AI1 = DCP*AI1+DSP*AR1
         AR1 = AR1H
         DO 124 IK=1,IDL1
            DO 1024 M=1,MP
            CH2(M,IK,L) = C2(M,IK,1)+AR1*C2(M,IK,2)
            CH2(M,IK,LC) = AI1*C2(M,IK,IP)
 1024       CONTINUE
  124    CONTINUE
         DC2 = AR1
         DS2 = AI1
         AR2 = AR1
         AI2 = AI1
         DO 126 J=3,IPPH
            JC = IPP2-J
            AR2H = DC2*AR2-DS2*AI2
            AI2 = DC2*AI2+DS2*AR2
            AR2 = AR2H
            DO 125 IK=1,IDL1
               DO 1025 M=1,MP
               CH2(M,IK,L) = CH2(M,IK,L)+AR2*C2(M,IK,J)
               CH2(M,IK,LC) = CH2(M,IK,LC)+AI2*C2(M,IK,JC)
 1025          CONTINUE
  125       CONTINUE
  126    CONTINUE
  127 CONTINUE
      DO 129 J=2,IPPH
         DO 128 IK=1,IDL1
            DO 1028 M=1,MP
            CH2(M,IK,1) = CH2(M,IK,1)+C2(M,IK,J)
 1028       CONTINUE
  128    CONTINUE
  129 CONTINUE
!
      IF (IDO .LT. L1) GO TO 132
      DO 131 K=1,L1
         DO 130 I=1,IDO
            DO 1030 M=1,MP
            CC(M,I,1,K) = CH(M,I,K,1)
 1030       CONTINUE
  130    CONTINUE
  131 CONTINUE
      GO TO 135
  132 DO 134 I=1,IDO
         DO 133 K=1,L1
            DO 1033 M=1,MP
            CC(M,I,1,K) = CH(M,I,K,1)
 1033       CONTINUE
  133    CONTINUE
  134 CONTINUE
  135 DO 137 J=2,IPPH
         JC = IPP2-J
         J2 = J+J
         DO 136 K=1,L1
            DO 1036 M=1,MP
            CC(M,IDO,J2-2,K) = CH(M,1,K,J)
            CC(M,1,J2-1,K) = CH(M,1,K,JC)
 1036       CONTINUE
  136    CONTINUE
  137 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IF (NBD .LT. L1) GO TO 141
      DO 140 J=2,IPPH
         JC = IPP2-J
         J2 = J+J
         DO 139 K=1,L1
            DO 138 I=3,IDO,2
               IC = IDP2-I
               DO 1038 M=1,MP
               CC(M,I-1,J2-1,K) = CH(M,I-1,K,J)+CH(M,I-1,K,JC)
               CC(M,IC-1,J2-2,K) = CH(M,I-1,K,J)-CH(M,I-1,K,JC)
               CC(M,I,J2-1,K) = CH(M,I,K,J)+CH(M,I,K,JC)
               CC(M,IC,J2-2,K) = CH(M,I,K,JC)-CH(M,I,K,J)
 1038          CONTINUE
  138       CONTINUE
  139    CONTINUE
  140 CONTINUE
      RETURN
  141 DO 144 J=2,IPPH
         JC = IPP2-J
         J2 = J+J
         DO 143 I=3,IDO,2
            IC = IDP2-I
            DO 142 K=1,L1
               DO 1042 M=1,MP
               CC(M,I-1,J2-1,K) = CH(M,I-1,K,J)+CH(M,I-1,K,JC)
               CC(M,IC-1,J2-2,K) = CH(M,I-1,K,J)-CH(M,I-1,K,JC)
               CC(M,I,J2-1,K) = CH(M,I,K,J)+CH(M,I,K,JC)
               CC(M,IC,J2-2,K) = CH(M,I,K,JC)-CH(M,I,K,J)
 1042          CONTINUE
  142       CONTINUE
  143    CONTINUE
  144 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VRFFTB(M,N,R,RT,MDIMR,WSAVE)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION     R(MDIMR,N),RT(MDIMR,N),WSAVE(N+15)
      IF (N .EQ. 1) RETURN
      CALL VRFTB1 (M,N,R,RT,MDIMR,WSAVE(1),WSAVE(N+1))
      RETURN
      END SUBROUTINE
      SUBROUTINE VRFFTF (M,N,R,RT,MDIMR,WSAVE)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION       R(MDIMR,N)  ,RT(MDIMR,N)    ,WSAVE(N+15)
      IF (N .EQ. 1) RETURN
      CALL VRFTF1 (M,N,R,RT,MDIMR,WSAVE(1),WSAVE(N+1))
      RETURN
      END SUBROUTINE
      SUBROUTINE VRFFTI (N,WSAVE)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION       WSAVE(N+15)
      IF (N .EQ. 1) RETURN
      CALL VRFTI1 (N,WSAVE(1),WSAVE(N+1))
      RETURN
      END SUBROUTINE
      SUBROUTINE VRFTB1 (M,N,C,CH,MDIMC,WA,FAC)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION       CH(MDIMC,N), C(MDIMC,N), WA(N) ,FAC(15)
      NF = FAC(2)
      NA = 0
      L1 = 1
      IW = 1
      DO 116 K1=1,NF
         IP = FAC(K1+2)
         L2 = IP*L1
         IDO = N/L2
         IDL1 = IDO*L1
         IF (IP .NE. 4) GO TO 103
         IX2 = IW+IDO
         IX3 = IX2+IDO
         IF (NA .NE. 0) GO TO 101
         CALL VRADB4 (M,IDO,L1,C,CH,MDIMC,WA(IW),WA(IX2),WA(IX3))
         GO TO 102
  101    CALL VRADB4 (M,IDO,L1,CH,C,MDIMC,WA(IW),WA(IX2),WA(IX3))
  102    NA = 1-NA
         GO TO 115
  103    IF (IP .NE. 2) GO TO 106
         IF (NA .NE. 0) GO TO 104
         CALL VRADB2 (M,IDO,L1,C,CH,MDIMC,WA(IW))
         GO TO 105
  104    CALL VRADB2 (M,IDO,L1,CH,C,MDIMC,WA(IW))
  105    NA = 1-NA
         GO TO 115
  106    IF (IP .NE. 3) GO TO 109
         IX2 = IW+IDO
         IF (NA .NE. 0) GO TO 107
         CALL VRADB3 (M,IDO,L1,C,CH,MDIMC,WA(IW),WA(IX2))
         GO TO 108
  107    CALL VRADB3 (M,IDO,L1,CH,C,MDIMC,WA(IW),WA(IX2))
  108    NA = 1-NA
         GO TO 115
  109    IF (IP .NE. 5) GO TO 112
         IX2 = IW+IDO
         IX3 = IX2+IDO
         IX4 = IX3+IDO
         IF (NA .NE. 0) GO TO 110
      CALL VRADB5 (M,IDO,L1,C,CH,MDIMC,WA(IW),WA(IX2),WA(IX3),WA(IX4))
         GO TO 111
  110 CALL VRADB5 (M,IDO,L1,CH,C,MDIMC,WA(IW),WA(IX2),WA(IX3),WA(IX4))
  111    NA = 1-NA
         GO TO 115
  112    IF (NA .NE. 0) GO TO 113
         CALL VRADBG (M,IDO,IP,L1,IDL1,C,C,C,CH,CH,MDIMC,WA(IW))
         GO TO 114
  113    CALL VRADBG (M,IDO,IP,L1,IDL1,CH,CH,CH,C,C,MDIMC,WA(IW))
  114    IF (IDO .EQ. 1) NA = 1-NA
  115    L1 = L2
         IW = IW+(IP-1)*IDO
  116 CONTINUE
      SCALE=SQRT(1.0d0/N)
      IF (NA .EQ. 0) GO TO 118
      DO 117 J=1,N
      DO 117 I=1,M
         C(I,J) = SCALE*CH(I,J)
  117 CONTINUE
      RETURN
  118 DO 119 J=1,N
      DO 119 I=1,M
         C(I,J)=SCALE*C(I,J)
  119 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VRFTF1 (M,N,C,CH,MDIMC,WA,FAC)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION       CH(MDIMC,N) ,C(MDIMC,N)  ,WA(N)   ,FAC(15)
      NF = FAC(2)
      NA = 1
      L2 = N
      IW = N
      DO 111 K1=1,NF
         KH = NF-K1
         IP = FAC(KH+3)
         L1 = L2/IP
         IDO = N/L2
         IDL1 = IDO*L1
         IW = IW-(IP-1)*IDO
         NA = 1-NA
         IF (IP .NE. 4) GO TO 102
         IX2 = IW+IDO
         IX3 = IX2+IDO
         IF (NA .NE. 0) GO TO 101
         CALL VRADF4 (M,IDO,L1,C,CH,MDIMC,WA(IW),WA(IX2),WA(IX3))
         GO TO 110
  101    CALL VRADF4 (M,IDO,L1,CH,C,MDIMC,WA(IW),WA(IX2),WA(IX3))
         GO TO 110
  102    IF (IP .NE. 2) GO TO 104
         IF (NA .NE. 0) GO TO 103
         CALL VRADF2 (M,IDO,L1,C,CH,MDIMC,WA(IW))
         GO TO 110
  103    CALL VRADF2 (M,IDO,L1,CH,C,MDIMC,WA(IW))
         GO TO 110
  104    IF (IP .NE. 3) GO TO 106
         IX2 = IW+IDO
         IF (NA .NE. 0) GO TO 105
         CALL VRADF3 (M,IDO,L1,C,CH,MDIMC,WA(IW),WA(IX2))
         GO TO 110
  105    CALL VRADF3 (M,IDO,L1,CH,C,MDIMC,WA(IW),WA(IX2))
         GO TO 110
  106    IF (IP .NE. 5) GO TO 108
         IX2 = IW+IDO
         IX3 = IX2+IDO
         IX4 = IX3+IDO
         IF (NA .NE. 0) GO TO 107
      CALL VRADF5(M,IDO,L1,C,CH,MDIMC,WA(IW),WA(IX2),WA(IX3),WA(IX4))
         GO TO 110
  107 CALL VRADF5 (M,IDO,L1,CH,C,MDIMC,WA(IW),WA(IX2),WA(IX3),WA(IX4))
         GO TO 110
  108    IF (IDO .EQ. 1) NA = 1-NA
         IF (NA .NE. 0) GO TO 109
         CALL VRADFG (M,IDO,IP,L1,IDL1,C,C,C,CH,CH,MDIMC,WA(IW))
         NA = 1
         GO TO 110
  109    CALL VRADFG (M,IDO,IP,L1,IDL1,CH,CH,CH,C,C,MDIMC,WA(IW))
         NA = 0
  110    L2 = L1
  111 CONTINUE
      SCALE=SQRT(1.0d0/N)
      IF (NA .EQ. 1) GO TO 113
      DO 112 J=1,N
      DO 112 I=1,M
         C(I,J) = SCALE*CH(I,J)
  112 CONTINUE
      RETURN
  113 DO 114 J=1,N
      DO 114 I=1,M
         C(I,J)=SCALE*C(I,J)
  114 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VRFTI1 (N,WA,FAC)
!
!     VRFFTPK, VERSION 1, AUGUST 1985
!
      DIMENSION       WA(N)      ,FAC(15)    ,NTRYH(4)
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
      NL = N
      NF = 0
      J = 0
  101 J = J+1
      IF (J-4) 102,102,103
  102 NTRY = NTRYH(J)
      GO TO 104
  103 NTRY = NTRY+2
  104 NQ = NL/NTRY
      NR = NL-NTRY*NQ
      IF (NR) 101,105,101
  105 NF = NF+1
      FAC(NF+2) = NTRY
      NL = NQ
      IF (NTRY .NE. 2) GO TO 107
      IF (NF .EQ. 1) GO TO 107
      DO 106 I=2,NF
         IB = NF-I+2
         FAC(IB+2) = FAC(IB+1)
  106 CONTINUE
      FAC(3) = 2
  107 IF (NL .NE. 1) GO TO 104
      FAC(1) = N
      FAC(2) = NF
      TPI = 2.0d0*3.141592653589793238462643
      ARGH = TPI/FLOAT(N)
      IS = 0
      NFM1 = NF-1
      L1 = 1
      IF (NFM1 .EQ. 0) RETURN
      DO 110 K1=1,NFM1
         IP = FAC(K1+2)
         LD = 0
         L2 = L1*IP
         IDO = N/L2
         IPM = IP-1
         DO 109 J=1,IPM
            LD = LD+L1
            I = IS
            ARGLD = FLOAT(LD)*ARGH
            FI = 0.0d0
            DO 108 II=3,IDO,2
               I = I+2
               FI = FI+1.0d0
               ARG = FI*ARGLD
               WA(I-1) = COS(ARG)
               WA(I) = SIN(ARG)
  108       CONTINUE
            IS = IS+IDO
  109    CONTINUE
         L1 = L2
  110 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VSINQB(M,N,X,XT,MDIMX,WSAVE)
!***BEGIN PROLOGUE  VSINQB
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, ODD WAVE
!             NUMBERS, MULTIPLE SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Normalized inverse of VSINQF.
!***DESCRIPTION
!
!  Subroutine VSINQB computes the backward fast Fourier sine transform
!  of M quarter wave sequences.  That is, sine series representations
!  with only odd wave numbers.  The transform is defined below at output
!  parameter X.
!
!  The array WSAVE which is used by subroutine VSINQB must be
!  initialized by calling subroutine VSINQI(N,WSAVE).
!
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequences to be transformed.  The method
!          is most efficient when N is a product of small primes.
!
!  X       an array of size at least X(MDIMX,N) which contains the
!          the sequences to be transformed.  The sequences are stored
!          in the ROWS of X.  Thus, the Jth sequence is stored in
!          X(J,I), I=1,..,N.
!
!  XT      a work array of size at least XT(MDIMX,N).
!
!  MDIMX   the first dimension of the array X exactly as it appears in
!          the calling program.
!
!  WSAVE   a work array which must be dimensioned at least 2*N+15
!          in the program that calls VSINQB.  The WSAVE array must be
!          initialized by calling subroutine VSINQI(N,WSAVE), and a
!          different WSAVE array must be used for each different
!          value of N.  This initialization does not have to be
!          repeated so long as N remains unchanged.
!
!  Output Parameters
!
!  X       for I=1,...,N and J=1,...,M
!
!               X(I)= the sum from K=1 to K=N of
!
!                 4*X(K)*SIN((2k-1)*I*PI/(2*N)) /SQRT(4*N)
!
!
!  WSAVE   contains initialization calculations which must not
!          be destroyed between calls of VSINQB or VSINQF.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VSINQF followed immediately by a call of
!           of VSINQB will return the original sequences X.  Thus,
!           VSINQB is the correctly normalized inverse VSINQF.
!
!  -----------------------------------------------------------------
!
!  VSINQB is a straightforward extension of the subprogram SINQB to
!  handle M simultaneous sequences.  SINQB was orginally developed
!  by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  VCOSQB
!***END PROLOGUE  VSINQB
      DIMENSION       X(MDIMX,*), XT(MDIMX,*), WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VSINQB
      IF (M .LE. 0)  GO TO 900
      IF (N .LE. 1)  GO TO 900
!
!  CASE  N .GT. 1
!
!     ... PREPROCESSING
!
      NS2 = N/2
      DO 210 K=2,N,2
         DO 210 J=1,M
            X(J,K) = -X(J,K)
  210 CONTINUE
!
!     ... COSINE QUARTER WAVE TRANSFORM
!
      CALL VCOSQB (M,N,X,XT,MDIMX,WSAVE)
!
!     ... POSTPROCESSING
!
      DO 220 K=1,NS2
         KC = N-K
         DO 220 J=1,M
            XHOLD = X(J,K)
            X(J,K) = X(J,KC+1)
            X(J,KC+1) = XHOLD
  220 CONTINUE
!
!  EXIT
!
  900 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VSINQF(M,N,X,XT,MDIMX,WSAVE)
!***BEGIN PROLOGUE  VSINQF
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, ODD WAVE
!             NUMBERS, MULTIPLE SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Forward sine transform, odd wave numbers, M sequences.
!***DESCRIPTION
!
!  Subroutine VSINQF computes the forward fast Fourier sine transform
!  of M quarter wave sequences.  That is, sine series representations
!  with only odd wave numbers.  The transform is defined below at output
!  parameter X.
!
!  The array WSAVE which is used by subroutine VSINQF must be
!  initialized by calling subroutine VSINQI(N,WSAVE).
!
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequences to be transformed.  The method
!          is most efficient when N is a product of small primes.
!
!  X       an array of size at least X(MDIMX,N) which contains the
!          the sequences to be transformed.  The sequences are stored
!          in the ROWS of X.  Thus, the Jth sequence is stored in
!          X(J,I), I=1,..,N.
!
!  XT      a work array of size at least XT(MDIMX,N).
!
!  MDIMX   the first dimension of the array X exactly as it appears in
!          the calling program.
!
!  WSAVE   a work array which must be dimensioned at least 2*N+15
!          in the program that calls VSINQF.  The WSAVE array must be
!          initialized by calling subroutine VSINQI(N,WSAVE), and a
!          different WSAVE array must be used for each different
!          value of N.  This initialization does not have to be
!          repeated so long as N remains unchanged.
!
!  Output Parameters
!
!  X       For I=1,...,N and J=1,...,M
!
!               X(I) = ( (-1)**(I-1)*X(N)
!
!                        + the sum from K=1 to K=N-1 of
!
!                        2*X(K)*SIN((2*I-1)*K*PI/(2*N)) )/SQRT(4*N)
!
!  WSAVE   contains initialization calculations which must not
!          be destroyed between calls of VSINQF or VSINQB.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VSINQF followed immediately by a call of
!           of VSINQB will return the original sequences X.  Thus,
!           VSINQB is the correctly normalized inverse VSINQF.
!
!  -----------------------------------------------------------------
!
!  VSINQF is a straightforward extension of the subprogram SINQF to
!  handle M simultaneous sequences.  SINQF was orginally developed
!  by P. N. Swarztrauber of NCAR.
!
 
!***ROUTINES CALLED  VCOSQF
!***END PROLOGUE  VSINQF
      DIMENSION       X(MDIMX,*), XT(MDIMX,*), WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VSINQF
      IF (M .LE. 0)  GO TO 900
      IF (N .LE. 1)  GO TO 900
!
!     ... PREPROCESSING
!
      NS2 = N/2
      DO 110 K=1,NS2
         KC = N-K
         DO 110 J=1,M
            XHOLD = X(J,K)
            X(J,K) = X(J,KC+1)
            X(J,KC+1) = XHOLD
  110 CONTINUE
!
!     ... COSINE QUARTER WAVE TRANSFORM
!
      CALL VCOSQF (M,N,X,XT,MDIMX,WSAVE)
!
!     ... POSTPROCESSING
!
      DO 120 K=2,N,2
         DO 120 J=1,M
            X(J,K) = -X(J,K)
  120 CONTINUE
!
!  EXIT
!
  900 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VSINQI(N,WSAVE)
!***BEGIN PROLOGUE  VSINQI
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, ODD WAVE
!             NUMBERS, MULTIPLE SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Initialize for VSINQF and VSINQB.
!***DESCRIPTION
!
!  Subroutine VSINQI initializes the array WSAVE which is used in
!  both VSINQF and VSINQB.  The prime factorization of N together with
!  a tabulation of the trigonometric functions are computed and
!  stored in WSAVE.
!
!  Input Parameter
!
!  N       the length of the sequence to be transformed.  The method
!          is most efficient when N is a product of small primes.
!
!  Output Parameter
!
!  WSAVE   a work array which must be dimensioned at least 3*N+15.
!          The same work array can be used for both VSINQF and VSINQB
!          as long as N remains unchanged.  Different WSAVE arrays
!          are required for different values of N.  The contents of
!
!          WSAVE must not be changed between calls of VSINQF or VSINQB.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  VCOSQI
!***END PROLOGUE  VSINQI
      DIMENSION       WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  VSINQI
      CALL VCOSQI (N,WSAVE)
      RETURN
      END SUBROUTINE
 
      SUBROUTINE VSINT(M,N,X,XT,MDIMX,WSAVE)
!***BEGIN PROLOGUE  VSINT
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, MULTIPLE
!             SEQUENCES
!***AUTHOR  BOISVERT, R. F., (NBS)
!***PURPOSE  Sine transform of one or more real, odd sequences.
!***DESCRIPTION
!
!  Subroutine VSINT computes the discrete Fourier sine transform
!  of M odd sequences X(J,I), J=1,...,M.  The transform is defined
!  below at output parameter X.
!
!  The array WSAVE which is used by subroutine VSINT must be
!  initialized by calling subroutine VSINTI(N,WSAVE).
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequence to be transformed.  The method
!          is most efficient when N+1 is the product of small primes.
!
!  X       an array of size at least X(MDIMX,N+1) which contains the
!          the sequences to be transformed.  The sequences are stored
!          in the ROWS of X.  Thus, the Jth sequence is stored in
!          X(J,I), I=1,..,N.  The extra column of X is used as work
!          storage.
!
!  XT      a work array of size at least XT(MDIMX,N+1).
!
!  MDIMX   the first dimension of the array X exactly as it appears in
!          the calling program.
!
!  WSAVE   a work array with dimension at least INT(2.5*N+15)
!          in the program that calls VSINT.  The WSAVE array must be
!          initialized by calling subroutine VSINTI(N,WSAVE), and a
!          different WSAVE array must be used for each different
!          value of N.  This initialization does not have to be
!          repeated so long as N remains unchanged.
!
!  Output Parameters
!
!  X       for I=1,...,N and J=1,...,M
!
!               X(J,I)= the sum from k=1 to k=N
!
!                    2*X(J,K)*SIN(K*I*PI/(N+1))/SQRT(2*(N+1))
!
!  WSAVE   contains initialization calculations which must not be
!          destroyed between calls of VSINT.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VSINT followed immediately by another call
!           of VSINT will return the original sequences X.  Thus,
!           VSINT is the correctly normalized inverse of itself.
!
!  -----------------------------------------------------------------
!
!  VSINT is a straightforward extension of the subprogram SINT to
!  handle M simultaneous sequences.  The scaling of the sequences
!  computed by VSINT is different than that of SINT.  SINT was
!  orginally developed by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  VRFFTF
!***END PROLOGUE  VSINT
      DIMENSION       X(MDIMX,*), XT(MDIMX,*), WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  SINT
      IF (M .LE. 0)  GO TO 900
      IF (N .LE. 1)  GO TO 900
      IF (N .GT. 2)  GO TO 300
!
!  CASE   N = 2
!
      SQRTH = SQRT(0.50d0)
      DO 201 J=1,M
         XH = SQRTH*(X(J,1)+X(J,2))
         X(J,2) = SQRTH*(X(J,1)-X(J,2))
         X(J,1) = XH
  201 CONTINUE
      GO TO 900
!
!  CASE   N .GT. 2
!
!     ... PREPROCESSING
!
  300 CONTINUE
      NP1 = N+1
      NS2 = N/2
      DO 301 J=1,M
         XT(J,1) = 0.0d0
  301 CONTINUE
      DO 310 K=1,NS2
         KC = NP1-K
         DO 310 J=1,M
            T1 = X(J,K)-X(J,KC)
            T2 = WSAVE(K)*(X(J,K)+X(J,KC))
            XT(J,K+1) = T1+T2
            XT(J,KC+1) = T2-T1
  310 CONTINUE
      MODN = MOD(N,2)
      IF (MODN .NE. 0) THEN
         DO 320 J=1,M
            XT(J,NS2+2) = 4.0d0*X(J,NS2+1)
  320    CONTINUE
      ENDIF
!
!     ... REAL*8 PERIODIC TRANSFORM
!
      NF = NS2+1
      CALL VRFFTF(M,NP1,XT,X,MDIMX,WSAVE(NF))
!
!     ... POSTPROCESSING
!
      DO 330 J=1,M
         X(J,1) = 0.5d0*XT(J,1)
  330 CONTINUE
      DO 350 I=3,N,2
         DO 340 J=1,M
            X(J,I-1) = -XT(J,I)
  340    CONTINUE
         DO 345 J=1,M
            X(J,I) = X(J,I-2)+XT(J,I-1)
  345    CONTINUE
  350 CONTINUE
      IF (MODN .EQ. 0) THEN
         DO 360 J=1,M
            X(J,N) = -XT(J,N+1)
  360    CONTINUE
      ENDIF
!
!     ... NORMALIZATION
!
      SCALE = SQRT(0.5d0)
      DO 370 I=1,N
         DO 370 J=1,M
            X(J,I) = SCALE*X(J,I)
  370 CONTINUE
!
!  EXIT
!
  900 CONTINUE
      RETURN
      END SUBROUTINE
      SUBROUTINE VSINTI(N,WSAVE)
!***BEGIN PROLOGUE  VSINTI
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  860701   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, SINE TRANSFORM, MULTIPLE
!             SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NBS)
!***PURPOSE  Initialize for VSINT.
!***DESCRIPTION
!
!  Subroutine VSINTI initializes the array WSAVE which is used in
!  subroutine SINT.  The prime factorization of N together with
!  a tabulation of the trigonometric functions are computed and
!  stored in WSAVE.
!
!  Input Parameter
!
!  N       the length of the sequence to be transformed.  The method
!          is most efficient when N+1 is a product of small primes.
!
!  Output Parameter
!
!  WSAVE   a work array with at least INT(2.5*N+15) locations.
!          Different WSAVE arrays are required for different values
!          of N.  The contents of WSAVE must not be changed between
!          calls of VSINT.
!
!  -----------------------------------------------------------------
!
!  VSINTI is a straightforward extension of the subprogram SINTI to
!  handle M simultaneous sequences.  SINTI was originally developed
!  P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!               pp. 51-83.
!***ROUTINES CALLED  VRFFTI
!***END PROLOGUE  VSINTI
      DIMENSION       WSAVE(*)
!***FIRST EXECUTABLE STATEMENT  SINTI
      PI = 4.0d0*ATAN(1.0d0)
      IF (N .LE. 1) RETURN
      NP1 = N+1
      NS2 = N/2
      DT = PI/REAL(NP1)
      KS = 1
      KF = KS+NS2-1
      FK = 0.0d0
      DO 101 K=KS,KF
         FK = FK+1.0d0
         WSAVE(K) = 2.0d0*SIN(FK*DT)
  101 CONTINUE
      CALL VRFFTI (NP1,WSAVE(KF+1))
      RETURN
      END SUBROUTINE

END MODULE BFFFT