LCOV - code coverage report
Current view: top level - math - rfft.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 469 861 54.5 %
Date: 2019-09-08 04:53:50 Functions: 12 14 85.7 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_rfft
       8             :    use m_juDFT
       9             :    PRIVATE
      10             : ! module also contains routines vrffti,vrfftf&vrfftb as private routines below
      11             :    PUBLIC rfft
      12             : CONTAINS
      13        9458 :    SUBROUTINE rfft( &
      14             :       isn,n1d,n2d,n3d,n1,n2,n3, &
      15       18916 :       nw1,nw2,nw3,wsave,b, &
      16        9458 :       a)
      17             : ! **********************************************************************
      18             : 
      19             : ! isn = +1 : FFT a real, 3dim array "a" of dimensions n1*n2*n3 to a complex
      20             : !            array of size n1*n2*(n3+1)/2 [if n3 is odd] or n1*n2*(n3/2+1)
      21             : !            [n3 is even].
      22             : ! isn = -1 : back-FFT of a complex array "a" of size n1*n2*(n3+1)/2 [odd]
      23             : !            or n1*n2*(n3/2+1) [even] to a real, 3dim array
      24             : 
      25             : ! the actual array is assumed to be located between
      26             : !                                  1 ... nw1+1    &     n1-nw1+1 ... n1
      27             : !                                  1 ... nw2+1    &     n2-nw2+1 ... n2
      28             : !                                  1 ... nw3+1    &     n3-nw3+1 ... n3
      29             : ! and padded with zeros in between.
      30             : ! if nw1 >= (n1-1)/2, no padding is assumed (-> aliasing errors!)
      31             : 
      32             : !                                                  G.Bihlmayer (UniWien)
      33             : 
      34             : ! **********************************************************************
      35             :       USE m_cfft
      36             :       IMPLICIT NONE
      37             : 
      38             :       INTEGER :: n1d,n2d,n3d,n1,n2,n3,nw1,nw2,nw3,isn
      39             :       REAL :: a(n1d,n2d,0:n3d),b(n1d,n2d,n3d),wsave(n3d+15)
      40             : 
      41             :       INTEGER :: i1,i2,i3,nup
      42             :       REAL :: factor
      43             :       LOGICAL :: l_nopad
      44             : 
      45             : ! a ... array for FFT
      46             : ! b ... work array
      47             : ! wsave ... stores tables for r-2-c FFT
      48             : ! n1,n2,n3 ... dimensions to be transformed
      49             : ! nw1,nw2,nw3 ... actual dimensions of a before FFT
      50             : ! n1d,n2d,n3d ... dimensions of a,b
      51             : 
      52             : ! check for input errors
      53             : 
      54        9458 :       IF ((isn/=-1) .AND. (isn /= 1))  CALL juDFT_error("choose isn=+/- 1" &
      55           0 :                                                         ,calledby ="rfft")
      56        9458 :       IF ((n1d < n1) .OR. (n2d < n2) .OR. (n3d < n3)) THEN
      57           0 :          WRITE (6,*) 'n1d,n2d,n3d =',n1d,n2d,n3d
      58           0 :          WRITE (6,*) 'n1 ,n2 ,n3  =',n1 ,n2 ,n3
      59           0 :          CALL juDFT_error("n(i) > n(i)d",calledby ="rfft")
      60             :       ENDIF
      61             :       IF ((n1 <= 2*nw1+1) .OR. &
      62        9458 :           (n2 <= 2*nw2+1) .OR. &
      63             :           (n3 <= 2*nw3+1)) THEN
      64             :          !        WRITE (6,*) 'n1 ,n2 ,n3  =',n1 ,n2 ,n3
      65             :          !        WRITE (6,*) 'nw1,nw2,nw3 =',nw1,nw2,nw3
      66             :          l_nopad= .TRUE.
      67             :       ELSE
      68             :          l_nopad= .FALSE.
      69             :       ENDIF
      70             : 
      71             : ! ******************** f o r w a r d - t r a n s f o r m *******************
      72             : 
      73        9458 :       IF (isn == 1) THEN
      74             : 
      75             :          ! array a is assumed to be zero from (1,1,0) to (n1d,n2d,0) and the array
      76             :          ! to be FFT'd starts at (1,1,1) as n1*n2*n3 real numbers.
      77             :          ! first transform n1*n2 real sequences of lenghth n3 to n3/2 complex values
      78             : 
      79        9108 :          CALL vrffti(n3,wsave)
      80        9108 :          IF (l_nopad) THEN
      81           0 :             CALL vrfftf(n1*n2,n3,a(1,1,1),b,n1d*n2d,wsave)
      82             :          ELSE
      83       52571 :             DO i2=1,nw2+1
      84       86926 :                CALL vrfftf(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
      85       96034 :                CALL vrfftf(nw1  ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
      86             :             ENDDO
      87       43463 :             DO i2=n2-nw2+1,n2
      88       68710 :                CALL vrfftf(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
      89       77818 :                CALL vrfftf(nw1  ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
      90             :             ENDDO
      91             :          ENDIF
      92             : 
      93             :          ! now we have the FFT'd array stored as described in vrfftf
      94             :          ! (mixed real & compex data)
      95             :          ! remove the norm 1/sqrt(n3) (to be compatible with cfft)
      96             : 
      97        9108 :          factor = sqrt(1.0*n3)
      98             :          ! ALL CPP_BLAS_sscal(n1d*n2d*n3,factor,a(1,1,1),1)
      99      492302 :          a=a*factor
     100             : 
     101             :          ! now, the real part of f(0) has to be moved to a(n1,n2,0) to get a purely
     102             :          ! complex array starting at a(1,1,0)
     103             : 
     104      146158 :          DO i1=1,n1
     105     2237042 :             DO i2=1,n2
     106     2090884 :                a(i1,i2,0)=a(i1,i2,1)
     107     2227934 :                a(i1,i2,1)=0.0
     108             :             ENDDO
     109             :          ENDDO
     110             : 
     111             :       ENDIF
     112             : 
     113             : ! ******************** g e n e r a l  p a r t *******************
     114             : 
     115             : ! now perform n2*n3/2 and n1*n3/2 complex FFT's; a is assumed to be
     116             : ! complex and starting at (1,1,0)
     117             : 
     118        9458 :       IF (ABS((n3/2.)-NINT(n3/2.)) > 0.1) THEN
     119             :          nup = n3
     120             :       ELSE
     121        9458 :          nup = n3+1
     122        9458 :          IF (n3+1>n3d)  CALL juDFT_error("n3 even & n3+1 > n3d" ,calledby &
     123           0 :                                          ="rfft")
     124        9458 :          a(:,:,n3+1)=0.0
     125             :       ENDIF
     126             : 
     127        9458 :       IF (l_nopad) THEN
     128        8016 :          DO i3=1,nup,2
     129        7666 :             CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n1,n1,isn)
     130        8016 :             CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n2,n1*n2,isn)
     131             :          ENDDO
     132             :       ELSE
     133      250705 :          DO i3=1,nup,2
     134      241597 :             CALL cfft(a(1,1,i3-1),a(1,1,i3),(nw2+1)*n1,n1,n1,isn)
     135             :             CALL cfft(a(1,n2-nw2+1,i3-1),a(1,n2-nw2+1,i3), &
     136      241597 :                       nw2*n1,n1,n1,isn)
     137      250705 :             CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n2,n1*n2,isn)
     138             :          ENDDO
     139             :       ENDIF
     140             : 
     141             : ! ******************** b a c k w a r d - t r a n s f o r m *******************
     142             : 
     143        9458 :       IF (isn == -1) THEN
     144             : 
     145             :          ! the real part of f(0) has to be moved to a(n1,n2,1) for a correct
     146             :          ! setup for vrfftb (see comments therein)
     147             : 
     148        5438 :          DO i1=1,n1
     149       80942 :             DO i2=1,n2
     150       75504 :                a(i1,i2,1)=a(i1,i2,0)
     151       80592 :                a(i1,i2,0)=0.0
     152             :             ENDDO
     153             :          ENDDO
     154             : 
     155             :          ! transform n1*n2 mixed real and complex sequences of lenghth n3/2
     156             :          ! to n3 real values
     157             : 
     158         350 :          CALL vrffti(n3,wsave)
     159         350 :          IF (l_nopad) THEN
     160         700 :             CALL vrfftb(n1*n2,n3,a(1,1,1),b,n1d*n2d,wsave)
     161             :          ELSE
     162           0 :             DO i2=1,nw2+1
     163           0 :                CALL vrfftb(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
     164           0 :                CALL vrfftb(nw1  ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
     165             :             ENDDO
     166           0 :             DO i2=n2-nw2+1,n2
     167           0 :                CALL vrfftb(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
     168           0 :                CALL vrfftb(nw1  ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
     169             :             ENDDO
     170             :          ENDIF
     171             : 
     172             :          ! remove the norm 1/sqrt(n3) (compatibility with cfft)
     173             : 
     174         350 :          factor = sqrt(1.0*n3)
     175             :          ! ALL CPP_BLAS_sscal(n1d*n2d*n3,factor,a(1,1,1),1)
     176       15682 :          a=a*factor
     177             : 
     178             :       ENDIF
     179             : 
     180        9458 :    END SUBROUTINE rfft
     181             : 
     182        9458 :    SUBROUTINE vrffti(n,wsave)
     183             : !***BEGIN PROLOGUE  VRFFTI
     184             : !***DATE WRITTEN   860701   (YYMMDD)
     185             : !***REVISION DATE  900509   (YYMMDD)
     186             : !***CATEGORY NO.  J1A1
     187             : !***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
     188             : !             MULTIPLE SEQUENCES
     189             : !***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
     190             : !***PURPOSE  Initialization for VRFFTF and VRFFTB.
     191             : !***DESCRIPTION
     192             : 
     193             : !  Subroutine VRFFTI initializes the array WSAVE which is used in
     194             : !  both VRFFTF and VRFFTB.  The prime factorization of N together with
     195             : !  a tabulation of certain trigonometric functions are computed and
     196             : !  stored in the array WSAVE.
     197             : 
     198             : !  Input Parameter
     199             : 
     200             : !  N       the length of the sequence to be transformed.  There is no
     201             : !          restriction on N.
     202             : 
     203             : !  Output Parameter
     204             : 
     205             : !  WSAVE   a work array which must be dimensioned at least N+15.
     206             : !          The same work array can be used for both VRFFTF and VRFFTB
     207             : !          as long as N remains unchanged.  Different WSAVE arrays
     208             : !          are required for different values of N.  The contents of
     209             : !          WSAVE must not be changed between calls of VRFFTF or VRFFTB.
     210             : 
     211             : !              * * * * * * * * * * * * * * * * * * * * *
     212             : !              *                                       *
     213             : !              *         PROGRAM SPECIFICATIONS        *
     214             : !              *                                       *
     215             : !              * * * * * * * * * * * * * * * * * * * * *
     216             : 
     217             : !     Dimension of    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
     218             : !     Arguments
     219             : 
     220             : !     Latest          AUGUST 1, 1985
     221             : !     Revision
     222             : 
     223             : !     Subprograms     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
     224             : !     Required        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
     225             : !                     VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
     226             : 
     227             : !     Special         NONE
     228             : !     Conditions
     229             : 
     230             : !     Common          NONE
     231             : !     blocks
     232             : 
     233             : !     I/O             NONE
     234             : 
     235             : !     Precision       SINGLE
     236             : 
     237             : !     Specialist      ROLAND SWEET
     238             : 
     239             : !     Language        FORTRAN
     240             : 
     241             : !     History         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
     242             : !                     NATIONAL BUREAU OF STANDARDS (BOULDER).
     243             : 
     244             : !     Algorithm       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
     245             : !                     OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
     246             : 
     247             : !     Portability     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
     248             : !                     THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
     249             : !                     THE FUNCTION PIMACH.
     250             : 
     251             : !     Required        COS,SIN
     252             : !     resident
     253             : !     Routines
     254             : 
     255             : !***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
     256             : !               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
     257             : !               pp. 51-83.
     258             : !***ROUTINES CALLED  VRFTI1
     259             : !***END PROLOGUE  VRFFTI
     260             : 
     261             : !     VRFFTPK, VERSION 1, AUGUST 1985
     262             : 
     263             :       DIMENSION wsave(n+15)
     264             : !***FIRST EXECUTABLE STATEMENT  VRFFTI
     265        9458 :       IF (n == 1) RETURN
     266        9458 :       CALL vrfti1(n,wsave(1),wsave(n+1))
     267             :       RETURN
     268             :    END subroutine
     269        9458 :    SUBROUTINE vrfti1(n,wa,fac)
     270             : 
     271             : !     VRFFTPK, VERSION 1, AUGUST 1985
     272             : 
     273             :       USE m_constants, ONLY : pimach
     274             :       DIMENSION wa(n),fac(15),ntryh(4)
     275             :       DATA ntryh(1),ntryh(2),ntryh(3),ntryh(4)/4,2,3,5/
     276             : 
     277        9458 :       nl = n
     278        9458 :       nf = 0
     279       10049 :       j = 0
     280       10049 : 10    j = j + 1
     281       10049 :       IF ( j <=4 ) THEN
     282       10049 :          ntry = ntryh(j)
     283             :       ELSE
     284           0 :          ntry = ntry + 2
     285             :       ENDIF
     286       26565 : 40    nq = nl/ntry
     287       26565 :       nr = nl - ntry*nq
     288       26565 :       IF ( nr /= 0 ) THEN
     289             :          GOTO 10
     290             :       ENDIF
     291       25974 :       nf = nf + 1
     292       25974 :       fac(nf+2) = ntry
     293       25974 :       nl = nq
     294       25974 :       IF (ntry /= 2) GO TO 70
     295         271 :       IF (nf == 1) GO TO 70
     296         198 :       DO 60 i = 2,nf
     297          66 :          ib = nf - i + 2
     298          66 :          fac(ib+2) = fac(ib+1)
     299          66 : 60    END DO
     300       25974 :       fac(3) = 2
     301       25974 : 70    IF (nl /= 1) GO TO 40
     302        9458 :       fac(1) = n
     303        9458 :       fac(2) = nf
     304        9458 :       tpi = 2.*pimach()
     305        9458 :       argh = tpi/real(n)
     306        9458 :       is = 0
     307        9458 :       nfm1 = nf - 1
     308        9458 :       l1 = 1
     309        9458 :       IF (nfm1 == 0) RETURN
     310       42490 :       DO 100 k1 = 1,nfm1
     311       16516 :          ip = fac(k1+2)
     312       16516 :          ld = 0
     313       16516 :          l2 = l1*ip
     314       16516 :          ido = n/l2
     315       16516 :          ipm = ip - 1
     316       65415 :          DO 90 j = 1,ipm
     317       48899 :             ld = ld + l1
     318       48899 :             i = is
     319       48899 :             argld = real(ld)*argh
     320       48899 :             fi = 0.
     321      221431 :             DO 80 ii = 3,ido,2
     322      172532 :                i = i + 2
     323      172532 :                fi = fi + 1.
     324      172532 :                arg = fi*argld
     325      172532 :                wa(i-1) = cos(arg)
     326      172532 :                wa(i) = sin(arg)
     327       48899 : 80          END DO
     328       48899 :             is = is + ido
     329       16516 : 90       END DO
     330       16516 :          l1 = l2
     331        9458 : 100   END DO
     332             :       RETURN
     333             :    END subroutine
     334      155636 :    SUBROUTINE vrfftf(m,n,r,rt,mdimr,wsave)
     335             : 
     336             : !***BEGIN PROLOGUE  VRFFTF
     337             : !***DATE WRITTEN   850801   (YYMMDD)
     338             : !***REVISION DATE  900509   (YYMMDD)
     339             : !***CATEGORY NO.  J1A1
     340             : !***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
     341             : !             FOURIER ANALYSIS, FORWARD TRANSFORM, MULTIPLE SEQUENCES
     342             : !***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
     343             : !***PURPOSE  Forward real periodic transform, M sequences.
     344             : !***DESCRIPTION
     345             : 
     346             : !  Subroutine VRFFTF computes the Fourier coefficients (forward
     347             : !  transform) of a number of real periodic sequences.  Specifically,
     348             : !  for each sequence the subroutine claculates the independent
     349             : !  Fourier coefficients described below at output parameter R.
     350             : 
     351             : !  The array WSAVE which is used by subroutine VRFFTF must be
     352             : !  initialized by calling subroutine VRFFTI(N,WSAVE).
     353             : 
     354             : !  Input Parameters
     355             : 
     356             : !  M       the number of sequences to be transformed.
     357             : 
     358             : !  N       the length of the sequences to be transformed.  The method
     359             : !          is most efficient when N is a product of small primes,
     360             : !          however n may be any positive integer.
     361             : 
     362             : !  R       areal two-dimensional array of size MDIMX x N containing the
     363             : !          the sequences to be transformed.  The sequences are stored
     364             : !          in the ROWS of R.  Thus, the I-th sequence to be transformed,
     365             : !          X(I,J), J=0,1,...,N-1, is stored as
     366             : 
     367             : !               R(I,J) = X(I,J-1) , J=1, 2, . . . , N.
     368             : 
     369             : !  RT      a real two-dimensional work array of size MDIMX x N.
     370             : 
     371             : !  MDIMR   the row (or first) dimension of the arrays R and RT exactly
     372             : !          as they appear in the calling program.  This parameter is
     373             : !          used to specify the variable dimension of these arrays.
     374             : 
     375             : !  WSAVE   a real one-dimensional work array which must be dimensioned
     376             : !          at least N+15.  The WSAVE array must be initialized by
     377             : !          calling subroutine VRFFTI.  A different WSAVE array must be
     378             : !          used for each different value of N.  This initialization does
     379             : !          not have to be repeated so long as N remains unchanged.  The
     380             : !          same WSAVE array may be used by VRFFTF and VRFFTB.
     381             : 
     382             : !  Output Parameters
     383             : 
     384             : !  R       contains the Fourier coefficients F(K) for each of the M
     385             : !          input sequences.  Specifically, row I of R, R(I,J),
     386             : !          J=1,2,..,N, contains the independent Fourier coefficients
     387             : !          F(I,K), for the I-th input sequence stored as
     388             : 
     389             : !             R(I,1) = REAL( F(I,0) ),
     390             : !                    = SQRT(1/N)*SUM(J=0,N-1)[ X(I,J) ],
     391             : 
     392             : !             R(I,2*K) = REAL( F(I,K) )
     393             : !                      = SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*COS(2J*K*PI/N)]
     394             : 
     395             : !             R(I,2*K+1) = IMAG( F(I,K) )
     396             : !                        =-SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*SIN(2J*K*PI/N)]
     397             : 
     398             : !                   for K = 1, 2, . . . , M-1,
     399             : 
     400             : !              and, when N is even,
     401             : 
     402             : !              R(I,N) = REAL( F(I,N/2) ).
     403             : !                     = SQRT(1/N)*SUM(J=0,N-1)[ (-1)**J*X(I,J) ].
     404             : 
     405             : !  WSAVE   contains results which must not be destroyed between calls
     406             : !          to VRFFTF or VRFFTB.
     407             : 
     408             : !  -----------------------------------------------------------------
     409             : 
     410             : !  NOTE  -  A call of VRFFTF followed immediately by a call of
     411             : !           of VRFFTB will return the original sequences R.  Thus,
     412             : !           VRFFTB is the correctly normalized inverse of VRFFTF.
     413             : 
     414             : !  -----------------------------------------------------------------
     415             : 
     416             : !  VRFFTF is a straightforward extension of the subprogram RFFTF to
     417             : !  handle M simultaneous sequences.  RFFTF was originally developed
     418             : !  by P. N. Swarztrauber of NCAR.
     419             : 
     420             : !              * * * * * * * * * * * * * * * * * * * * *
     421             : !              *                                       *
     422             : !              *         PROGRAM SPECIFICATIONS        *
     423             : !              *                                       *
     424             : !              * * * * * * * * * * * * * * * * * * * * *
     425             : 
     426             : !     Dimension of    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
     427             : !     Arguments
     428             : 
     429             : !     Latest          AUGUST 1, 1985
     430             : !     Revision
     431             : 
     432             : !     Subprograms     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
     433             : !     Required        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
     434             : !                     VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
     435             : 
     436             : !     Special         NONE
     437             : !     Conditions
     438             : 
     439             : !     Common          NONE
     440             : !     blocks
     441             : 
     442             : !     I/O             NONE
     443             : 
     444             : !     Precision       SINGLE
     445             : 
     446             : !     Specialist      ROLAND SWEET
     447             : 
     448             : !     Language        FORTRAN
     449             : 
     450             : !     History         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
     451             : !                     NATIONAL BUREAU OF STANDARDS (BOULDER).
     452             : 
     453             : !     Algorithm       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
     454             : !                     OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
     455             : 
     456             : !     Portability     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
     457             : !                     THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
     458             : !                     THE FUNCTION PIMACH.
     459             : 
     460             : !     Required        COS,SIN
     461             : !     resident
     462             : !     Routines
     463             : 
     464             : !***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
     465             : !               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
     466             : !               pp. 51-83.
     467             : !***ROUTINES CALLED  VRFTF1
     468             : !***END PROLOGUE  VRFFTF
     469             : 
     470             : !     VRFFTPK, VERSION 1, AUGUST 1985
     471             : 
     472             :       DIMENSION r(mdimr,n),rt(mdimr,n),wsave(n+15)
     473             : !***FIRST EXECUTABLE STATEMENT  VRFFTF
     474      155636 :       IF (n == 1) RETURN
     475      155636 :       CALL vrftf1(m,n,r,rt,mdimr,wsave(1),wsave(n+1))
     476             :       RETURN
     477             :    END subroutine
     478        5042 :    SUBROUTINE vradf2(mp,ido,l1,cc,ch,mdimc,wa1)
     479             : 
     480             : !     VRFFTPK, VERSION 1, AUGUST 1985
     481             : 
     482             :       DIMENSION ch(mdimc,ido,2,l1),cc(mdimc,ido,l1,2),wa1(ido)
     483             : 
     484       10084 :       DO 20 k = 1,l1
     485       31801 :          DO 10 m = 1,mp
     486       26759 :             ch(m,1,1,k) = cc(m,1,k,1) + cc(m,1,k,2)
     487       26759 :             ch(m,ido,2,k) = cc(m,1,k,1) - cc(m,1,k,2)
     488        5042 : 10       END DO
     489        5042 : 20    END DO
     490        5042 :       IF ( ido ==  2 ) THEN
     491             :          GOTO 70
     492        5042 :       ELSEIF ( ido < 2 ) THEN
     493             :          GOTO 100
     494             :       ENDIF
     495        5042 :       idp2 = ido + 2
     496       10084 :       DO 60 k = 1,l1
     497       33398 :          DO 50 i = 3,ido,2
     498       28356 :             ic = idp2 - i
     499      179454 :             DO 40 m = 1,mp
     500             :                ch(m,i,1,k) = cc(m,i,k,1) + &
     501             :                              (wa1(i-2)*cc(m,i,k,2)-wa1(i-1)* &
     502      151098 :                               cc(m,i-1,k,2))
     503             :                ch(m,ic,2,k) = (wa1(i-2)*cc(m,i,k,2)- &
     504      151098 :                                wa1(i-1)*cc(m,i-1,k,2)) - cc(m,i,k,1)
     505             :                ch(m,i-1,1,k) = cc(m,i-1,k,1) + &
     506             :                                (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)* &
     507      151098 :                                 cc(m,i,k,2))
     508             :                ch(m,ic-1,2,k) = cc(m,i-1,k,1) - &
     509             :                                 (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)* &
     510      151098 :                                  cc(m,i,k,2))
     511       28356 : 40          END DO
     512        5042 : 50       END DO
     513        5042 : 60    END DO
     514        5042 :       IF (mod(ido,2) == 1) RETURN
     515        2916 : 70    DO 90 k = 1,l1
     516        5346 :          DO 80 m = 1,mp
     517        4374 :             ch(m,1,2,k) = -cc(m,ido,k,2)
     518        4374 :             ch(m,ido,1,k) = cc(m,ido,k,1)
     519         972 : 80       END DO
     520         972 : 90    END DO
     521             : 100   RETURN
     522             :    END subroutine
     523        7308 :    SUBROUTINE vradf3(mp,ido,l1,cc,ch,mdimc,wa1,wa2)
     524             : 
     525             : !     VRFFTPK, VERSION 1, AUGUST 1985
     526             : 
     527             :       USE m_constants, ONLY : pimach
     528             :       DIMENSION ch(mdimc,ido,3,l1),cc(mdimc,ido,l1,3),wa1(ido),wa2(ido)
     529             : 
     530        7308 :       arg = 2.*pimach()/3.
     531        7308 :       taur = cos(arg)
     532        7308 :       taui = sin(arg)
     533       40428 :       DO 20 k = 1,l1
     534      207504 :          DO 10 m = 1,mp
     535      174384 :             ch(m,1,1,k) = cc(m,1,k,1) + (cc(m,1,k,2)+cc(m,1,k,3))
     536      174384 :             ch(m,1,3,k) = taui* (cc(m,1,k,3)-cc(m,1,k,2))
     537             :             ch(m,ido,2,k) = cc(m,1,k,1) + &
     538      174384 :                             taur* (cc(m,1,k,2)+cc(m,1,k,3))
     539       33120 : 10       END DO
     540        7308 : 20    END DO
     541        7308 :       IF (ido == 1) RETURN
     542        3168 :       idp2 = ido + 2
     543        9504 :       DO 50 k = 1,l1
     544       12672 :          DO 40 i = 3,ido,2
     545        6336 :             ic = idp2 - i
     546       41184 :             DO 30 m = 1,mp
     547             :                ch(m,i-1,1,k) = cc(m,i-1,k,1) + &
     548             :                ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i, &
     549             :                k,2))+ (wa2(i-2)*cc(m,i-1,k, &
     550       34848 :                &                          3)+wa2(i-1)*cc(m,i,k,3)))
     551             :                ch(m,i,1,k) = cc(m,i,k,1) + &
     552             :                ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
     553             :                &                        2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
     554       34848 :                i-1,k,3)))
     555             :                ch(m,i-1,3,k) = (cc(m,i-1,k,1)+ &
     556             :                taur* ((wa1(i-2)*cc(m,i-1,k, &
     557             :                &                          2)+wa1(i-1)*cc(m,i,k,2))+ (wa2(i-2)*cc(m, &
     558             :                i-1,k,3)+wa2(i-1)*cc(m,i,k,3)))) + &
     559             :                (taui* ((wa1(i-2)*cc(m,i,k, &
     560             :                &                          2)-wa1(i-1)*cc(m,i-1,k, &
     561             :                &                          2))- (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
     562       34848 :                i-1,k,3))))
     563             :                ch(m,ic-1,2,k) = (cc(m,i-1,k,1)+ &
     564             :                taur* ((wa1(i-2)*cc(m,i-1,k, &
     565             :                &                           2)+wa1(i-1)*cc(m,i,k, &
     566             :                &                           2))+ (wa2(i-2)*cc(m,i-1,k, &
     567             :                &                           3)+wa2(i-1)*cc(m,i,k,3)))) - &
     568             :                (taui* ((wa1(i-2)*cc(m,i,k, &
     569             :                &                           2)-wa1(i-1)*cc(m,i-1,k, &
     570             :                &                           2))- (wa2(i-2)*cc(m,i,k, &
     571       34848 :                &                           3)-wa2(i-1)*cc(m,i-1,k,3))))
     572             :                ch(m,i,3,k) = (cc(m,i,k,1)+taur* &
     573             :                ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
     574             :                &                        2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
     575             :                i-1,k,3)))) + (taui* ((wa2(i-2)*cc(m,i-1,k, &
     576             :                &                        3)+wa2(i-1)*cc(m,i,k,3))- (wa1(i-2)*cc(m, &
     577       34848 :                i-1,k,2)+wa1(i-1)*cc(m,i,k,2))))
     578             :                ch(m,ic,2,k) = (taui* ((wa2(i-2)*cc(m,i-1,k, &
     579             :                &                         3)+wa2(i-1)*cc(m,i,k,3))- (wa1(i-2)*cc(m, &
     580             :                i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))) - &
     581             :                (cc(m,i,k,1)+taur* ((wa1(i-2)*cc(m,i,k, &
     582             :                &                         2)-wa1(i-1)*cc(m,i-1,k, &
     583             :                &                         2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
     584       34848 :                i-1,k,3))))
     585        6336 : 30          END DO
     586        6336 : 40       END DO
     587        3168 : 50    END DO
     588             :       RETURN
     589             :    END subroutine
     590      421086 :    SUBROUTINE vradf4(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3)
     591             : 
     592             : !     VRFFTPK, VERSION 1, AUGUST 1985
     593             : 
     594             :       DIMENSION cc(mdimc,ido,l1,4),ch(mdimc,ido,4,l1),wa1(ido),wa2(ido), &
     595             :          wa3(ido)
     596             : 
     597      421086 :       hsqt2 = sqrt(2.)/2.
     598     3078816 :       DO 20 k = 1,l1
     599    14459175 :          DO 10 m = 1,mp
     600             :             ch(m,1,1,k) = (cc(m,1,k,2)+cc(m,1,k,4)) + &
     601    11801445 :                           (cc(m,1,k,1)+cc(m,1,k,3))
     602             :             ch(m,ido,4,k) = (cc(m,1,k,1)+cc(m,1,k,3)) - &
     603    11801445 :                             (cc(m,1,k,2)+cc(m,1,k,4))
     604    11801445 :             ch(m,ido,2,k) = cc(m,1,k,1) - cc(m,1,k,3)
     605    11801445 :             ch(m,1,3,k) = cc(m,1,k,4) - cc(m,1,k,2)
     606     2657730 : 10       END DO
     607      421086 : 20    END DO
     608      421086 :       IF ( ido ==  2 ) THEN
     609             :          GOTO 70
     610      421086 :       ELSEIF ( ido < 2 ) THEN
     611             :          GOTO 100
     612             :       ENDIF
     613      270492 :       idp2 = ido + 2
     614      898734 :       DO 60 k = 1,l1
     615     1970040 :          DO 50 i = 3,ido,2
     616     1341798 :             ic = idp2 - i
     617     7348221 :             DO 40 m = 1,mp
     618             :                ch(m,i-1,1,k) = ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i, &
     619             :                k,2))+ (wa3(i-2)*cc(m,i-1,k, &
     620             :                &                          4)+wa3(i-1)*cc(m,i,k,4))) + &
     621             :                (cc(m,i-1,k,1)+ (wa2(i-2)*cc(m,i-1,k, &
     622     6006423 :                &                          3)+wa2(i-1)*cc(m,i,k,3)))
     623             :                ch(m,ic-1,4,k) = (cc(m,i-1,k,1)+ &
     624             :                (wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i, &
     625             :                k,3))) - ((wa1(i-2)*cc(m,i-1,k, &
     626             :                &                           2)+wa1(i-1)*cc(m,i,k,2))+ &
     627             :                (wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i, &
     628     6006423 :                k,4)))
     629             :                ch(m,i,1,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
     630             :                &                        2))+ (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
     631             :                i-1,k,4))) + (cc(m,i,k,1)+ &
     632             :                (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k, &
     633     6006423 :                &                        3)))
     634             :                ch(m,ic,4,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1, &
     635             :                k,2))+ (wa3(i-2)*cc(m,i,k, &
     636             :                &                         4)-wa3(i-1)*cc(m,i-1,k,4))) - &
     637             :                (cc(m,i,k,1)+ (wa2(i-2)*cc(m,i,k, &
     638     6006423 :                &                         3)-wa2(i-1)*cc(m,i-1,k,3)))
     639             :                ch(m,i-1,3,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1, &
     640             :                k,2))- (wa3(i-2)*cc(m,i,k, &
     641             :                &                          4)-wa3(i-1)*cc(m,i-1,k,4))) + &
     642             :                (cc(m,i-1,k,1)- (wa2(i-2)*cc(m,i-1,k, &
     643     6006423 :                &                          3)+wa2(i-1)*cc(m,i,k,3)))
     644             :                ch(m,ic-1,2,k) = (cc(m,i-1,k,1)- &
     645             :                (wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i, &
     646             :                k,3))) - ((wa1(i-2)*cc(m,i,k, &
     647             :                &                           2)-wa1(i-1)*cc(m,i-1,k,2))- &
     648             :                (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1, &
     649     6006423 :                k,4)))
     650             :                ch(m,i,3,k) = ((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,k, &
     651             :                &                        4))- (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m, &
     652             :                i,k,2))) + (cc(m,i,k,1)- &
     653             :                (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k, &
     654     6006423 :                &                        3)))
     655             :                ch(m,ic,2,k) = ((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i, &
     656             :                k,4))- (wa1(i-2)*cc(m,i-1,k, &
     657             :                &                         2)+wa1(i-1)*cc(m,i,k,2))) - &
     658             :                (cc(m,i,k,1)- (wa2(i-2)*cc(m,i,k, &
     659     6006423 :                &                         3)-wa2(i-1)*cc(m,i-1,k,3)))
     660     1341798 : 40          END DO
     661      628242 : 50       END DO
     662      270492 : 60    END DO
     663      270492 :       IF (mod(ido,2) == 1) RETURN
     664             : 70    CONTINUE
     665     1522116 :       DO 90 k = 1,l1
     666     3412971 :          DO 80 m = 1,mp
     667             :             ch(m,ido,1,k) = (hsqt2* (cc(m,ido,k,2)-cc(m,ido,k,4))) + &
     668     2786673 :                             cc(m,ido,k,1)
     669             :             ch(m,ido,3,k) = cc(m,ido,k,1) - &
     670     2786673 :                             (hsqt2* (cc(m,ido,k,2)-cc(m,ido,k,4)))
     671             :             ch(m,1,2,k) = (-hsqt2* (cc(m,ido,k,2)+cc(m,ido,k,4))) - &
     672     2786673 :                           cc(m,ido,k,3)
     673             :             ch(m,1,4,k) = (-hsqt2* (cc(m,ido,k,2)+cc(m,ido,k,4))) + &
     674     2786673 :                           cc(m,ido,k,3)
     675      626298 : 80       END DO
     676      420114 : 90    END DO
     677             : 100   RETURN
     678             :    END subroutine
     679        1804 :    SUBROUTINE vradf5(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3,wa4)
     680             : 
     681             : !     VRFFTPK, VERSION 1, AUGUST 1985
     682             : 
     683             :       USE m_constants, ONLY : pimach
     684             :       DIMENSION cc(mdimc,ido,l1,5),ch(mdimc,ido,5,l1),wa1(ido),wa2(ido), &
     685             :          wa3(ido),wa4(ido)
     686             : 
     687        1804 :       arg = 2.*pimach()/5.
     688        1804 :       tr11 = cos(arg)
     689        1804 :       ti11 = sin(arg)
     690        1804 :       tr12 = cos(2.*arg)
     691        1804 :       ti12 = sin(2.*arg)
     692       12628 :       DO 20 k = 1,l1
     693       70356 :          DO 10 m = 1,mp
     694             :             ch(m,1,1,k) = cc(m,1,k,1) + (cc(m,1,k,5)+cc(m,1,k,2)) + &
     695       59532 :                           (cc(m,1,k,4)+cc(m,1,k,3))
     696             :             ch(m,ido,2,k) = cc(m,1,k,1) + &
     697             :                             tr11* (cc(m,1,k,5)+cc(m,1,k,2)) + &
     698       59532 :                             tr12* (cc(m,1,k,4)+cc(m,1,k,3))
     699             :             ch(m,1,3,k) = ti11* (cc(m,1,k,5)-cc(m,1,k,2)) + &
     700       59532 :                           ti12* (cc(m,1,k,4)-cc(m,1,k,3))
     701             :             ch(m,ido,4,k) = cc(m,1,k,1) + &
     702             :                             tr12* (cc(m,1,k,5)+cc(m,1,k,2)) + &
     703       59532 :                             tr11* (cc(m,1,k,4)+cc(m,1,k,3))
     704             :             ch(m,1,5,k) = ti12* (cc(m,1,k,5)-cc(m,1,k,2)) - &
     705       59532 :                           ti11* (cc(m,1,k,4)-cc(m,1,k,3))
     706       10824 : 10       END DO
     707        1804 : 20    END DO
     708        1804 :       IF (ido == 1) RETURN
     709         902 :       idp2 = ido + 2
     710        2706 :       DO 50 k = 1,l1
     711        5412 :          DO 40 i = 3,ido,2
     712        3608 :             ic = idp2 - i
     713       23452 :             DO 30 m = 1,mp
     714             :                ch(m,i-1,1,k) = cc(m,i-1,k,1) + &
     715             :                ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i, &
     716             :                k,2))+ (wa4(i-2)*cc(m,i-1,k, &
     717             :                &                          5)+wa4(i-1)*cc(m,i,k,5))) + &
     718             :                ((wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i, &
     719             :                k,3))+ (wa3(i-2)*cc(m,i-1,k, &
     720       19844 :                &                          4)+wa3(i-1)*cc(m,i,k,4)))
     721             :                ch(m,i,1,k) = cc(m,i,k,1) + &
     722             :                ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
     723             :                &                        2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
     724             :                i-1,k,5))) + ((wa2(i-2)*cc(m,i,k, &
     725             :                &                        3)-wa2(i-1)*cc(m,i-1,k,3))+ &
     726             :                (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1,k, &
     727       19844 :                &                        4)))
     728             :                ch(m,i-1,3,k) = cc(m,i-1,k,1) + &
     729             :                                tr11* (wa1(i-2)*cc(m,i-1,k,2)+ &
     730             :                                       wa1(i-1)*cc(m,i,k,2)+ &
     731             :                                       wa4(i-2)*cc(m,i-1,k,5)+ &
     732             :                                       wa4(i-1)*cc(m,i,k,5)) + &
     733             :                                tr12* (wa2(i-2)*cc(m,i-1,k,3)+ &
     734             :                                       wa2(i-1)*cc(m,i,k,3)+ &
     735             :                                       wa3(i-2)*cc(m,i-1,k,4)+ &
     736             :                                       wa3(i-1)*cc(m,i,k,4)) + &
     737             :                                ti11* (wa1(i-2)*cc(m,i,k,2)- &
     738             :                                       wa1(i-1)*cc(m,i-1,k,2)- &
     739             :                                       (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m,i-1, &
     740             :                                                                         k,5))) + ti12* (wa2(i-2)*cc(m,i,k,3)- &
     741             :                                                                                         wa2(i-1)*cc(m,i-1,k,3)- &
     742             :                                                                                         (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1, &
     743       19844 :                                                                                                                           k,4)))
     744             :                ch(m,ic-1,2,k) = cc(m,i-1,k,1) + &
     745             :                tr11* (wa1(i-2)*cc(m,i-1,k,2)+ &
     746             :                wa1(i-1)*cc(m,i,k,2)+ &
     747             :                wa4(i-2)*cc(m,i-1,k,5)+ &
     748             :                wa4(i-1)*cc(m,i,k,5)) + &
     749             :                tr12* (wa2(i-2)*cc(m,i-1,k,3)+ &
     750             :                wa2(i-1)*cc(m,i,k,3)+ &
     751             :                wa3(i-2)*cc(m,i-1,k,4)+ &
     752             :                wa3(i-1)*cc(m,i,k,4)) - &
     753             :                (ti11* (wa1(i-2)*cc(m,i,k, &
     754             :                &                           2)-wa1(i-1)*cc(m,i-1,k, &
     755             :                &                           2)- (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
     756             :                i-1,k,5)))+ti12* (wa2(i-2)*cc(m,i,k, &
     757             :                &                           3)-wa2(i-1)*cc(m,i-1,k, &
     758             :                &                           3)- (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
     759       19844 :                i-1,k,4))))
     760             :                ch(m,i,3,k) = (cc(m,i,k,1)+tr11* &
     761             :                ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
     762             :                &                        2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
     763             :                i-1,k,5)))+tr12* ((wa2(i-2)*cc(m,i,k, &
     764             :                &                        3)-wa2(i-1)*cc(m,i-1,k,3))+ (wa3(i-2)*cc(m, &
     765             :                i,k,4)-wa3(i-1)*cc(m,i-1,k,4)))) + &
     766             :                (ti11* ((wa4(i-2)*cc(m,i-1,k, &
     767             :                &                        5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
     768             :                i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))+ &
     769             :                ti12* ((wa3(i-2)*cc(m,i-1,k, &
     770             :                &                        4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
     771       19844 :                i-1,k,3)+wa2(i-1)*cc(m,i,k,3))))
     772             :                ch(m,ic,2,k) = (ti11* ((wa4(i-2)*cc(m,i-1,k, &
     773             :                &                         5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
     774             :                i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))+ &
     775             :                ti12* ((wa3(i-2)*cc(m,i-1,k, &
     776             :                &                         4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
     777             :                i-1,k,3)+wa2(i-1)*cc(m,i,k,3)))) - &
     778             :                (cc(m,i,k,1)+tr11* ((wa1(i-2)*cc(m,i,k, &
     779             :                &                         2)-wa1(i-1)*cc(m,i-1,k, &
     780             :                &                         2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
     781             :                i-1,k,5)))+tr12* ((wa2(i-2)*cc(m,i,k, &
     782             :                &                         3)-wa2(i-1)*cc(m,i-1,k, &
     783             :                &                         3))+ (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
     784       19844 :                i-1,k,4))))
     785             :                ch(m,i-1,5,k) = (cc(m,i-1,k,1)+ &
     786             :                tr12* ((wa1(i-2)*cc(m,i-1,k, &
     787             :                &                          2)+wa1(i-1)*cc(m,i,k,2))+ (wa4(i-2)*cc(m, &
     788             :                i-1,k,5)+wa4(i-1)*cc(m,i,k,5)))+ &
     789             :                tr11* ((wa2(i-2)*cc(m,i-1,k, &
     790             :                &                          3)+wa2(i-1)*cc(m,i,k,3))+ (wa3(i-2)*cc(m, &
     791             :                i-1,k,4)+wa3(i-1)*cc(m,i,k,4)))) + &
     792             :                (ti12* ((wa1(i-2)*cc(m,i,k, &
     793             :                &                          2)-wa1(i-1)*cc(m,i-1,k, &
     794             :                &                          2))- (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
     795             :                i-1,k,5)))-ti11* ((wa2(i-2)*cc(m,i,k, &
     796             :                &                          3)-wa2(i-1)*cc(m,i-1,k, &
     797             :                &                          3))- (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
     798       19844 :                i-1,k,4))))
     799             :                ch(m,ic-1,4,k) = (cc(m,i-1,k,1)+ &
     800             :                tr12* ((wa1(i-2)*cc(m,i-1,k, &
     801             :                &                           2)+wa1(i-1)*cc(m,i,k, &
     802             :                &                           2))+ (wa4(i-2)*cc(m,i-1,k, &
     803             :                &                           5)+wa4(i-1)*cc(m,i,k,5)))+ &
     804             :                tr11* ((wa2(i-2)*cc(m,i-1,k, &
     805             :                &                           3)+wa2(i-1)*cc(m,i,k, &
     806             :                &                           3))+ (wa3(i-2)*cc(m,i-1,k, &
     807             :                &                           4)+wa3(i-1)*cc(m,i,k,4)))) - &
     808             :                (ti12* ((wa1(i-2)*cc(m,i,k, &
     809             :                &                           2)-wa1(i-1)*cc(m,i-1,k, &
     810             :                &                           2))- (wa4(i-2)*cc(m,i,k, &
     811             :                &                           5)-wa4(i-1)*cc(m,i-1,k,5)))- &
     812             :                ti11* ((wa2(i-2)*cc(m,i,k, &
     813             :                &                           3)-wa2(i-1)*cc(m,i-1,k, &
     814             :                &                           3))- (wa3(i-2)*cc(m,i,k, &
     815       19844 :                &                           4)-wa3(i-1)*cc(m,i-1,k,4))))
     816             :                ch(m,i,5,k) = (cc(m,i,k,1)+tr12* &
     817             :                ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
     818             :                &                        2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
     819             :                i-1,k,5)))+tr11* ((wa2(i-2)*cc(m,i,k, &
     820             :                &                        3)-wa2(i-1)*cc(m,i-1,k,3))+ (wa3(i-2)*cc(m, &
     821             :                i,k,4)-wa3(i-1)*cc(m,i-1,k,4)))) + &
     822             :                (ti12* ((wa4(i-2)*cc(m,i-1,k, &
     823             :                &                        5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
     824             :                i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))- &
     825             :                ti11* ((wa3(i-2)*cc(m,i-1,k, &
     826             :                &                        4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
     827       19844 :                i-1,k,3)+wa2(i-1)*cc(m,i,k,3))))
     828             :                ch(m,ic,4,k) = (ti12* ((wa4(i-2)*cc(m,i-1,k, &
     829             :                &                         5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
     830             :                i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))- &
     831             :                ti11* ((wa3(i-2)*cc(m,i-1,k, &
     832             :                &                         4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
     833             :                i-1,k,3)+wa2(i-1)*cc(m,i,k,3)))) - &
     834             :                (cc(m,i,k,1)+tr12* ((wa1(i-2)*cc(m,i,k, &
     835             :                &                         2)-wa1(i-1)*cc(m,i-1,k, &
     836             :                &                         2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
     837             :                i-1,k,5)))+tr11* ((wa2(i-2)*cc(m,i,k, &
     838             :                &                         3)-wa2(i-1)*cc(m,i-1,k, &
     839             :                &                         3))+ (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
     840       19844 :                i-1,k,4))))
     841        3608 : 30          END DO
     842        1804 : 40       END DO
     843         902 : 50    END DO
     844             :       RETURN
     845             :    END subroutine
     846           0 :    SUBROUTINE vradfg(mp,ido,ip,l1,idl1,cc,c1,c2,ch,ch2,mdimc,wa)
     847             : 
     848             : !     VRFFTPK, VERSION 1, AUGUST 1985
     849             : 
     850             :       USE m_constants, ONLY : pimach
     851             :       DIMENSION ch(mdimc,ido,l1,ip),cc(mdimc,ido,ip,l1), &
     852             :          c1(mdimc,ido,l1,ip),c2(mdimc,idl1,ip), &
     853             :          ch2(mdimc,idl1,ip),wa(ido)
     854             : 
     855           0 :       tpi = 2.*pimach()
     856           0 :       arg = tpi/real(ip)
     857           0 :       dcp = cos(arg)
     858           0 :       dsp = sin(arg)
     859           0 :       ipph = (ip+1)/2
     860           0 :       ipp2 = ip + 2
     861           0 :       idp2 = ido + 2
     862           0 :       nbd = (ido-1)/2
     863           0 :       IF (ido == 1) GO TO 250
     864           0 :       DO 20 ik = 1,idl1
     865           0 :          DO 10 m = 1,mp
     866           0 :             ch2(m,ik,1) = c2(m,ik,1)
     867           0 : 10       END DO
     868           0 : 20    END DO
     869           0 :       DO 50 j = 2,ip
     870           0 :          DO 40 k = 1,l1
     871           0 :             DO 30 m = 1,mp
     872           0 :                ch(m,1,k,j) = c1(m,1,k,j)
     873           0 : 30          END DO
     874           0 : 40       END DO
     875           0 : 50    END DO
     876           0 :       IF (nbd > l1) GO TO 100
     877           0 :       is = -ido
     878           0 :       DO 90 j = 2,ip
     879           0 :          is = is + ido
     880           0 :          idij = is
     881           0 :          DO 80 i = 3,ido,2
     882           0 :             idij = idij + 2
     883           0 :             DO 70 k = 1,l1
     884           0 :                DO 60 m = 1,mp
     885             :                   ch(m,i-1,k,j) = wa(idij-1)*c1(m,i-1,k,j) + &
     886           0 :                                   wa(idij)*c1(m,i,k,j)
     887             :                   ch(m,i,k,j) = wa(idij-1)*c1(m,i,k,j) - &
     888           0 :                                 wa(idij)*c1(m,i-1,k,j)
     889           0 : 60             END DO
     890           0 : 70          END DO
     891           0 : 80       END DO
     892           0 : 90    END DO
     893           0 :       GO TO 150
     894           0 : 100   is = -ido
     895           0 :       DO 140 j = 2,ip
     896           0 :          is = is + ido
     897           0 :          DO 130 k = 1,l1
     898           0 :             idij = is
     899           0 :             DO 120 i = 3,ido,2
     900           0 :                idij = idij + 2
     901           0 :                DO 110 m = 1,mp
     902             :                   ch(m,i-1,k,j) = wa(idij-1)*c1(m,i-1,k,j) + &
     903           0 :                                   wa(idij)*c1(m,i,k,j)
     904             :                   ch(m,i,k,j) = wa(idij-1)*c1(m,i,k,j) - &
     905           0 :                                 wa(idij)*c1(m,i-1,k,j)
     906           0 : 110            END DO
     907           0 : 120         END DO
     908           0 : 130      END DO
     909           0 : 140   END DO
     910           0 : 150   IF (nbd < l1) GO TO 200
     911           0 :       DO 190 j = 2,ipph
     912           0 :          jc = ipp2 - j
     913           0 :          DO 180 k = 1,l1
     914           0 :             DO 170 i = 3,ido,2
     915           0 :                DO 160 m = 1,mp
     916           0 :                   c1(m,i-1,k,j) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
     917           0 :                   c1(m,i-1,k,jc) = ch(m,i,k,j) - ch(m,i,k,jc)
     918           0 :                   c1(m,i,k,j) = ch(m,i,k,j) + ch(m,i,k,jc)
     919           0 :                   c1(m,i,k,jc) = ch(m,i-1,k,jc) - ch(m,i-1,k,j)
     920           0 : 160            END DO
     921           0 : 170         END DO
     922           0 : 180      END DO
     923           0 : 190   END DO
     924           0 :       GO TO 280
     925           0 : 200   DO 240 j = 2,ipph
     926           0 :          jc = ipp2 - j
     927           0 :          DO 230 i = 3,ido,2
     928           0 :             DO 220 k = 1,l1
     929           0 :                DO 210 m = 1,mp
     930           0 :                   c1(m,i-1,k,j) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
     931           0 :                   c1(m,i-1,k,jc) = ch(m,i,k,j) - ch(m,i,k,jc)
     932           0 :                   c1(m,i,k,j) = ch(m,i,k,j) + ch(m,i,k,jc)
     933           0 :                   c1(m,i,k,jc) = ch(m,i-1,k,jc) - ch(m,i-1,k,j)
     934           0 : 210            END DO
     935           0 : 220         END DO
     936           0 : 230      END DO
     937           0 : 240   END DO
     938           0 :       GO TO 280
     939           0 : 250   DO 270 ik = 1,idl1
     940           0 :          DO 260 m = 1,mp
     941           0 :             c2(m,ik,1) = ch2(m,ik,1)
     942           0 : 260      END DO
     943           0 : 270   END DO
     944           0 : 280   DO 310 j = 2,ipph
     945           0 :          jc = ipp2 - j
     946           0 :          DO 300 k = 1,l1
     947           0 :             DO 290 m = 1,mp
     948           0 :                c1(m,1,k,j) = ch(m,1,k,j) + ch(m,1,k,jc)
     949           0 :                c1(m,1,k,jc) = ch(m,1,k,jc) - ch(m,1,k,j)
     950           0 : 290         END DO
     951           0 : 300      END DO
     952           0 : 310   END DO
     953             : 
     954             :       ar1 = 1.
     955             :       ai1 = 0.
     956           0 :       DO 370 l = 2,ipph
     957           0 :          lc = ipp2 - l
     958           0 :          ar1h = dcp*ar1 - dsp*ai1
     959           0 :          ai1 = dcp*ai1 + dsp*ar1
     960           0 :          ar1 = ar1h
     961           0 :          DO 330 ik = 1,idl1
     962           0 :             DO 320 m = 1,mp
     963           0 :                ch2(m,ik,l) = c2(m,ik,1) + ar1*c2(m,ik,2)
     964           0 :                ch2(m,ik,lc) = ai1*c2(m,ik,ip)
     965           0 : 320         END DO
     966           0 : 330      END DO
     967             :          dc2 = ar1
     968             :          ds2 = ai1
     969             :          ar2 = ar1
     970             :          ai2 = ai1
     971           0 :          DO 360 j = 3,ipph
     972           0 :             jc = ipp2 - j
     973           0 :             ar2h = dc2*ar2 - ds2*ai2
     974           0 :             ai2 = dc2*ai2 + ds2*ar2
     975           0 :             ar2 = ar2h
     976           0 :             DO 350 ik = 1,idl1
     977           0 :                DO 340 m = 1,mp
     978           0 :                   ch2(m,ik,l) = ch2(m,ik,l) + ar2*c2(m,ik,j)
     979           0 :                   ch2(m,ik,lc) = ch2(m,ik,lc) + ai2*c2(m,ik,jc)
     980           0 : 340            END DO
     981           0 : 350         END DO
     982           0 : 360      END DO
     983           0 : 370   END DO
     984           0 :       DO 400 j = 2,ipph
     985           0 :          DO 390 ik = 1,idl1
     986           0 :             DO 380 m = 1,mp
     987           0 :                ch2(m,ik,1) = ch2(m,ik,1) + c2(m,ik,j)
     988           0 : 380         END DO
     989           0 : 390      END DO
     990           0 : 400   END DO
     991             : 
     992           0 :       IF (ido < l1) GO TO 440
     993           0 :       DO 430 k = 1,l1
     994           0 :          DO 420 i = 1,ido
     995           0 :             DO 410 m = 1,mp
     996           0 :                cc(m,i,1,k) = ch(m,i,k,1)
     997           0 : 410         END DO
     998           0 : 420      END DO
     999           0 : 430   END DO
    1000           0 :       GO TO 480
    1001           0 : 440   DO 470 i = 1,ido
    1002           0 :          DO 460 k = 1,l1
    1003           0 :             DO 450 m = 1,mp
    1004           0 :                cc(m,i,1,k) = ch(m,i,k,1)
    1005           0 : 450         END DO
    1006           0 : 460      END DO
    1007           0 : 470   END DO
    1008           0 : 480   DO 510 j = 2,ipph
    1009           0 :          jc = ipp2 - j
    1010           0 :          j2 = j + j
    1011           0 :          DO 500 k = 1,l1
    1012           0 :             DO 490 m = 1,mp
    1013           0 :                cc(m,ido,j2-2,k) = ch(m,1,k,j)
    1014           0 :                cc(m,1,j2-1,k) = ch(m,1,k,jc)
    1015           0 : 490         END DO
    1016           0 : 500      END DO
    1017           0 : 510   END DO
    1018           0 :       IF (ido == 1) RETURN
    1019           0 :       IF (nbd < l1) GO TO 560
    1020           0 :       DO 550 j = 2,ipph
    1021           0 :          jc = ipp2 - j
    1022           0 :          j2 = j + j
    1023           0 :          DO 540 k = 1,l1
    1024           0 :             DO 530 i = 3,ido,2
    1025           0 :                ic = idp2 - i
    1026           0 :                DO 520 m = 1,mp
    1027           0 :                   cc(m,i-1,j2-1,k) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
    1028           0 :                   cc(m,ic-1,j2-2,k) = ch(m,i-1,k,j) - ch(m,i-1,k,jc)
    1029           0 :                   cc(m,i,j2-1,k) = ch(m,i,k,j) + ch(m,i,k,jc)
    1030           0 :                   cc(m,ic,j2-2,k) = ch(m,i,k,jc) - ch(m,i,k,j)
    1031           0 : 520            END DO
    1032           0 : 530         END DO
    1033           0 : 540      END DO
    1034           0 : 550   END DO
    1035           0 :       RETURN
    1036           0 : 560   DO 600 j = 2,ipph
    1037           0 :          jc = ipp2 - j
    1038           0 :          j2 = j + j
    1039           0 :          DO 590 i = 3,ido,2
    1040           0 :             ic = idp2 - i
    1041           0 :             DO 580 k = 1,l1
    1042           0 :                DO 570 m = 1,mp
    1043           0 :                   cc(m,i-1,j2-1,k) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
    1044           0 :                   cc(m,ic-1,j2-2,k) = ch(m,i-1,k,j) - ch(m,i-1,k,jc)
    1045           0 :                   cc(m,i,j2-1,k) = ch(m,i,k,j) + ch(m,i,k,jc)
    1046           0 :                   cc(m,ic,j2-2,k) = ch(m,i,k,jc) - ch(m,i,k,j)
    1047           0 : 570            END DO
    1048           0 : 580         END DO
    1049           0 : 590      END DO
    1050           0 : 600   END DO
    1051             :       RETURN
    1052             :    END subroutine
    1053      155636 :    SUBROUTINE vrftf1(m,n,c,ch,mdimc,wa,fac)
    1054             : 
    1055             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1056             : 
    1057             :       DIMENSION ch(mdimc,n),c(mdimc,n),wa(n),fac(15)
    1058             : 
    1059      155636 :       nf = fac(2)
    1060      155636 :       na = 1
    1061      155636 :       l2 = n
    1062      155636 :       iw = n
    1063      590876 :       DO 110 k1 = 1,nf
    1064      435240 :          kh = nf - k1
    1065      435240 :          ip = fac(kh+3)
    1066      435240 :          l1 = l2/ip
    1067      435240 :          ido = n/l2
    1068      435240 :          idl1 = ido*l1
    1069      435240 :          iw = iw - (ip-1)*ido
    1070      435240 :          na = 1 - na
    1071      435240 :          IF (ip /= 4) GO TO 20
    1072      421086 :          ix2 = iw + ido
    1073      421086 :          ix3 = ix2 + ido
    1074      421086 :          IF (na /= 0) GO TO 10
    1075      269520 :          CALL vradf4(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3))
    1076      269520 :          GO TO 100
    1077      151566 : 10       CALL vradf4(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3))
    1078      151566 :          GO TO 100
    1079       14154 : 20       IF (ip /= 2) GO TO 40
    1080        5042 :          IF (na /= 0) GO TO 30
    1081        5042 :          CALL vradf2(m,ido,l1,c,ch,mdimc,wa(iw))
    1082        5042 :          GO TO 100
    1083           0 : 30       CALL vradf2(m,ido,l1,ch,c,mdimc,wa(iw))
    1084           0 :          GO TO 100
    1085        9112 : 40       IF (ip /= 3) GO TO 60
    1086        7308 :          ix2 = iw + ido
    1087        7308 :          IF (na /= 0) GO TO 50
    1088        4140 :          CALL vradf3(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2))
    1089        4140 :          GO TO 100
    1090        3168 : 50       CALL vradf3(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2))
    1091        3168 :          GO TO 100
    1092        1804 : 60       IF (ip /= 5) GO TO 80
    1093        1804 :          ix2 = iw + ido
    1094        1804 :          ix3 = ix2 + ido
    1095        1804 :          ix4 = ix3 + ido
    1096        1804 :          IF (na /= 0) GO TO 70
    1097         902 :          CALL vradf5(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    1098         902 :          GO TO 100
    1099         902 : 70       CALL vradf5(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    1100         902 :          GO TO 100
    1101           0 : 80       IF (ido == 1) na = 1 - na
    1102           0 :          IF (na /= 0) GO TO 90
    1103           0 :          CALL vradfg(m,ido,ip,l1,idl1,c,c,c,ch,ch,mdimc,wa(iw))
    1104           0 :          na = 1
    1105           0 :          GO TO 100
    1106           0 : 90       CALL vradfg(m,ido,ip,l1,idl1,ch,ch,ch,c,c,mdimc,wa(iw))
    1107           0 :          na = 0
    1108      435240 : 100      l2 = l1
    1109      155636 : 110   END DO
    1110      155636 :       scale = sqrt(1./n)
    1111      155636 :       IF (na == 1) GO TO 140
    1112    15597400 :       DO 130 j = 1,n
    1113    42654062 :          DO 120 i = 1,m
    1114    34917346 :             c(i,j) = scale*ch(i,j)
    1115     7736716 : 120      END DO
    1116      123968 : 130   END DO
    1117       31668 :       RETURN
    1118     1045044 : 140   DO 160 j = 1,n
    1119     2280096 :          DO 150 i = 1,m
    1120     1773408 :             c(i,j) = scale*c(i,j)
    1121      506688 : 150      END DO
    1122       31668 : 160   END DO
    1123             :       RETURN
    1124             :    END subroutine
    1125             : 
    1126         350 :    SUBROUTINE vrfftb(m,n,r,rt,mdimr,wsave)
    1127             : !***BEGIN PROLOGUE  VRFFTB
    1128             : !***DATE WRITTEN   850801   (YYMMDD)
    1129             : !***REVISION DATE  900509   (YYMMDD)
    1130             : !***CATEGORY NO.  J1A1
    1131             : !***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
    1132             : !             FOURIER SYNTHESIS, BACKWARD TRANSFORM, MULTIPLE SEQUENCES
    1133             : !***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
    1134             : !***PURPOSE  Backward real periodic transform, M sequences.
    1135             : !***DESCRIPTION
    1136             : 
    1137             : !  Subroutine VRFFTB computes the synthesis (backward transform) of a
    1138             : !  number of real periodic sequences from their Fourier coefficients.
    1139             : !  Specifically, for each set of independent Fourier coefficients
    1140             : !  F(K), the corresponding real periodic sequence is computed.
    1141             : 
    1142             : !  The array WSAVE which is used by subroutine VRFFTB must be
    1143             : !  initialized by calling subroutine VRFFTI(N,WSAVE).
    1144             : 
    1145             : !  Input Parameters
    1146             : 
    1147             : !  M       the number of sets of coefficients.
    1148             : 
    1149             : !  N       the length of the sequences of coefficients to be
    1150             : !          transformed.  The method is most efficient when N is a
    1151             : !          product of small primes, however n may be any positive
    1152             : !          integer.
    1153             : 
    1154             : !  R       areal two-dimensional array of size MDIMX x N containing the
    1155             : !          coefficients to be transformed.  Each set of coefficients
    1156             : !          F(K), K\0,1,..,N-1, is stored as a ROW of R.  Specifically,
    1157             : !          the I-th set of independent Fourier coefficients is stored
    1158             : 
    1159             : !                R(I,1) = REAL( F(I,0) ),
    1160             : 
    1161             : !                R(I,2*K) = REAL( F(I,K) )
    1162             : 
    1163             : !                R(I,2*K+1) = IMAG( F(I,K) )
    1164             : 
    1165             : !                   for K = 1, 2, . . . , M-1,
    1166             : 
    1167             : !                and, when N is even,
    1168             : 
    1169             : !                R(I,N) = REAL( F(I,N/2) ).
    1170             : 
    1171             : !  RT      a real two-dimensional work array of size MDIMX x N.
    1172             : 
    1173             : !  MDIMR   the row (or first) dimension of the arrays R and RT exactly
    1174             : !          as they appear in the calling program.  This parameter is
    1175             : !          used to specify the variable dimension of these arrays.
    1176             : 
    1177             : !  WSAVE   a real one-dimensional work array which must be dimensioned
    1178             : !          at least N+15.  The WSAVE array must be initialized by
    1179             : !          calling subroutine VRFFTI.  A different WSAVE array must be
    1180             : !          used for each different value of N.  This initialization does
    1181             : !          not have to be repeated so long as N remains unchanged.  The
    1182             : !          same WSAVE array may be used by VRFFTB and VRFFTB.
    1183             : 
    1184             : !  Output Parameters
    1185             : 
    1186             : !  R       contains M real periodic sequences corresponding to the given
    1187             : !          coefficients.  Specifically, the I-th row of R contains the
    1188             : !          real periodic sequence corresponding to the I-th set of
    1189             : !          independent Fourier coefficients F(I,K) stored as
    1190             : 
    1191             : !               R(I,J) = X(I,J-1) ,   J = 1, 2, . . . , N, where
    1192             : 
    1193             : !               X(I,J) = SQRT(1/N)* F(I,0) + (-1)**J*F(I,N/2)
    1194             : !                        + 2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N)
    1195             : !                        - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ]  ,
    1196             : 
    1197             : !                 when N is even, and
    1198             : 
    1199             : !               X(I,J) = SQRT(1/N)* F(I,0) +
    1200             : !                        2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N)
    1201             : !                        - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ]  ,
    1202             : 
    1203             : !                 when N is odd.
    1204             : 
    1205             : !  WSAVE   contains results which must not be destroyed between calls
    1206             : !          to VRFFTF or VRFFTB.
    1207             : 
    1208             : !  -----------------------------------------------------------------
    1209             : 
    1210             : !  NOTE  -  A call of VRFFTF followed immediately by a call of
    1211             : !           of VRFFTB will return the original sequences R.  Thus,
    1212             : !           VRFFTB is the correctly normalized inverse of VRFFTF.
    1213             : 
    1214             : !  -----------------------------------------------------------------
    1215             : 
    1216             : !  VRFFTB is a straightforward extension of the subprogram RFFTB to
    1217             : !  handle M simultaneous sequences.  RFFTB was originally developed
    1218             : !  by P. N. Swarztrauber of NCAR.
    1219             : 
    1220             : !              * * * * * * * * * * * * * * * * * * * * *
    1221             : !              *                                       *
    1222             : !              *         PROGRAM SPECIFICATIONS        *
    1223             : !              *                                       *
    1224             : !              * * * * * * * * * * * * * * * * * * * * *
    1225             : 
    1226             : !     Dimension of    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
    1227             : !     Arguments
    1228             : 
    1229             : !     Latest          AUGUST 1, 1985
    1230             : !     Revision
    1231             : 
    1232             : !     Subprograms     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
    1233             : !     required        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
    1234             : !                     VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
    1235             : 
    1236             : !     Special         NONE
    1237             : !     Conditions
    1238             : 
    1239             : !     Common          NONE
    1240             : !     blocks
    1241             : 
    1242             : !     I/O             NONE
    1243             : 
    1244             : !     Precision       SINGLE
    1245             : 
    1246             : !     Specialist      ROLAND SWEET
    1247             : 
    1248             : !     Language        FORTRAN
    1249             : 
    1250             : !     History         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
    1251             : !                     NATIONAL BUREAU OF STANDARDS (BOULDER).
    1252             : 
    1253             : !     Algorithm       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
    1254             : !                     OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
    1255             : 
    1256             : !     Portability     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
    1257             : !                     THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
    1258             : !                     THE FUNCTION PIMACH.
    1259             : 
    1260             : !     Required        COS,SIN
    1261             : !     resident
    1262             : !     Routines
    1263             : 
    1264             : !***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
    1265             : !               Computations, (G. Rodrigue, ed.), Academic Press, 1982,
    1266             : !               pp. 51-83.
    1267             : !***ROUTINES CALLED  VRFTB1
    1268             : !***END PROLOGUE  VRFFTB
    1269             : 
    1270             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1271             : 
    1272             :       IMPLICIT NONE
    1273             :       INTEGER, INTENT(IN) :: n, m, mdimr
    1274             :       REAL :: r(mdimr,n),rt(mdimr,n),wsave(n+15)
    1275             : 
    1276         350 :       IF (n == 1) RETURN
    1277         350 :       CALL vrftb1(m,n,r,rt,mdimr,wsave(1),wsave(n+1))
    1278             :       RETURN
    1279             : 
    1280             :    CONTAINS
    1281             : 
    1282          32 :       SUBROUTINE vradb2(mp,ido,l1,cc,ch,mdimc,wa1)
    1283             : 
    1284             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1285             : 
    1286             :          IMPLICIT NONE
    1287             :          INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
    1288             :          REAL, INTENT(IN)    :: cc(mdimc,ido,2,l1), wa1(ido)
    1289             :          REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,2)
    1290             : 
    1291             :          INTEGER :: i, m, k, ic, idp2
    1292             : 
    1293          64 :          DO 20 k = 1,l1
    1294        9584 :             DO 10 m = 1,mp
    1295        9552 :                ch(m,1,k,1) = cc(m,1,1,k) + cc(m,ido,2,k)
    1296        9552 :                ch(m,1,k,2) = cc(m,1,1,k) - cc(m,ido,2,k)
    1297          32 : 10          END DO
    1298          32 : 20       END DO
    1299          32 :          IF ( ido ==  2 ) THEN
    1300             :             GOTO 70
    1301          32 :          ELSEIF ( ido < 2 ) THEN
    1302             :             GOTO 100
    1303             :          ENDIF
    1304          32 :          idp2 = ido + 2
    1305          64 :          DO 60 k = 1,l1
    1306         236 :             DO 50 i = 3,ido,2
    1307         204 :                ic = idp2 - i
    1308       62220 :                DO 40 m = 1,mp
    1309       62016 :                   ch(m,i-1,k,1) = cc(m,i-1,1,k) + cc(m,ic-1,2,k)
    1310       62016 :                   ch(m,i,k,1) = cc(m,i,1,k) - cc(m,ic,2,k)
    1311             :                   ch(m,i-1,k,2) = wa1(i-2)* (cc(m,i-1,1,k)- &
    1312             :                                              cc(m,ic-1,2,k)) - wa1(i-1)* &
    1313       62016 :                                   (cc(m,i,1,k)+cc(m,ic,2,k))
    1314             :                   ch(m,i,k,2) = wa1(i-2)* (cc(m,i,1,k)+cc(m,ic,2,k)) + &
    1315       62016 :                                 wa1(i-1)* (cc(m,i-1,1,k)-cc(m,ic-1,2,k))
    1316         204 : 40             END DO
    1317          32 : 50          END DO
    1318          32 : 60       END DO
    1319          32 :          IF (mod(ido,2) == 1) RETURN
    1320          36 : 70       DO 90 k = 1,l1
    1321        3084 :             DO 80 m = 1,mp
    1322        3072 :                ch(m,ido,k,1) = cc(m,ido,1,k) + cc(m,ido,1,k)
    1323        3072 :                ch(m,ido,k,2) = - (cc(m,1,2,k)+cc(m,1,2,k))
    1324          12 : 80          END DO
    1325          12 : 90       END DO
    1326             : 100      RETURN
    1327             :       END SUBROUTINE vradb2
    1328             : 
    1329          36 :       SUBROUTINE vradb3(mp,ido,l1,cc,ch,mdimc,wa1,wa2)
    1330             : 
    1331             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1332             : 
    1333             :          USE m_constants, ONLY : pimach
    1334             :          IMPLICIT NONE
    1335             :          INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
    1336             :          REAL, INTENT(IN)    :: cc(mdimc,ido,3,l1), wa1(ido), wa2(ido)
    1337             :          REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,3)
    1338             : 
    1339             :          INTEGER :: i, k, idp2, ic, m
    1340             :          REAL :: arg, taur, taui
    1341             : 
    1342          36 :          arg = 2.*pimach()/3.
    1343          36 :          taur = cos(arg)
    1344          36 :          taui = sin(arg)
    1345         228 :          DO 20 k = 1,l1
    1346       55872 :             DO 10 m = 1,mp
    1347       55680 :                ch(m,1,k,1) = cc(m,1,1,k) + 2.*cc(m,ido,2,k)
    1348             :                ch(m,1,k,2) = cc(m,1,1,k) + (2.*taur)*cc(m,ido,2,k) - &
    1349       55680 :                              (2.*taui)*cc(m,1,3,k)
    1350             :                ch(m,1,k,3) = cc(m,1,1,k) + (2.*taur)*cc(m,ido,2,k) + &
    1351       55680 :                &                     2.*taui*cc(m,1,3,k)
    1352         192 : 10          END DO
    1353          36 : 20       END DO
    1354          36 :          IF (ido == 1) RETURN
    1355          12 :          idp2 = ido + 2
    1356          36 :          DO 50 k = 1,l1
    1357          48 :             DO 40 i = 3,ido,2
    1358          24 :                ic = idp2 - i
    1359        7800 :                DO 30 m = 1,mp
    1360             :                   ch(m,i-1,k,1) = cc(m,i-1,1,k) + &
    1361        7776 :                                   (cc(m,i-1,3,k)+cc(m,ic-1,2,k))
    1362        7776 :                   ch(m,i,k,1) = cc(m,i,1,k) + (cc(m,i,3,k)-cc(m,ic,2,k))
    1363             :                   ch(m,i-1,k,2) = wa1(i-2)* ((cc(m,i-1,1,k)+taur* (cc(m, &
    1364             :                                                                       i-1,3,k)+cc(m,ic-1,2,k)))- &
    1365             :                                              (taui* (cc(m,i,3,k)+cc(m,ic,2,k)))) - &
    1366             :                                   wa1(i-1)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
    1367             :                                                                     k)-cc(m,ic,2,k)))+ (taui* (cc(m,i-1,3, &
    1368        7776 :                                                                                                   k)-cc(m,ic-1,2,k))))
    1369             :                   ch(m,i,k,2) = wa1(i-2)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
    1370             :                   k)-cc(m,ic,2,k)))+ (taui* (cc(m,i-1,3, &
    1371             :                   k)-cc(m,ic-1,2,k)))) + &
    1372             :                   wa1(i-1)* ((cc(m,i-1,1,k)+taur* (cc(m,i-1, &
    1373             :                   &                        3,k)+cc(m,ic-1,2,k)))- &
    1374        7776 :                   (taui* (cc(m,i,3,k)+cc(m,ic,2,k))))
    1375             :                   ch(m,i-1,k,3) = wa2(i-2)* ((cc(m,i-1,1,k)+taur* (cc(m, &
    1376             :                                                                       i-1,3,k)+cc(m,ic-1,2,k)))+ &
    1377             :                                              (taui* (cc(m,i,3,k)+cc(m,ic,2,k)))) - &
    1378             :                                   wa2(i-1)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
    1379             :                                                                     k)-cc(m,ic,2,k)))- (taui* (cc(m,i-1,3, &
    1380        7776 :                                                                                                   k)-cc(m,ic-1,2,k))))
    1381             :                   ch(m,i,k,3) = wa2(i-2)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
    1382             :                   k)-cc(m,ic,2,k)))- (taui* (cc(m,i-1,3, &
    1383             :                   k)-cc(m,ic-1,2,k)))) + &
    1384             :                   wa2(i-1)* ((cc(m,i-1,1,k)+taur* (cc(m,i-1, &
    1385             :                   &                        3,k)+cc(m,ic-1,2,k)))+ &
    1386        7776 :                   (taui* (cc(m,i,3,k)+cc(m,ic,2,k))))
    1387          24 : 30             END DO
    1388          24 : 40          END DO
    1389          12 : 50       END DO
    1390             :          RETURN
    1391             :       END SUBROUTINE vradb3
    1392             : 
    1393         828 :       SUBROUTINE vradb4(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3)
    1394             : 
    1395             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1396             : 
    1397             :          IMPLICIT NONE
    1398             :          INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
    1399             :          REAL, INTENT(IN)    :: cc(mdimc,ido,4,l1), wa1(ido), wa2(ido), &
    1400             :                                 wa3(ido)
    1401             :          REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,4)
    1402             : 
    1403             :          INTEGER :: k, m, i, idp2, ic
    1404             :          REAL :: sqrt2
    1405             : 
    1406         828 :          sqrt2 = sqrt(2.)
    1407        5322 :          DO 20 k = 1,l1
    1408     1077678 :             DO 10 m = 1,mp
    1409             :                ch(m,1,k,3) = (cc(m,1,1,k)+cc(m,ido,4,k)) - &
    1410     1073184 :                              (cc(m,ido,2,k)+cc(m,ido,2,k))
    1411             :                ch(m,1,k,1) = (cc(m,1,1,k)+cc(m,ido,4,k)) + &
    1412     1073184 :                              (cc(m,ido,2,k)+cc(m,ido,2,k))
    1413             :                ch(m,1,k,4) = (cc(m,1,1,k)-cc(m,ido,4,k)) + &
    1414     1073184 :                              (cc(m,1,3,k)+cc(m,1,3,k))
    1415             :                ch(m,1,k,2) = (cc(m,1,1,k)-cc(m,ido,4,k)) - &
    1416     1073184 :                              (cc(m,1,3,k)+cc(m,1,3,k))
    1417        4494 : 10          END DO
    1418         828 : 20       END DO
    1419         828 :          IF ( ido ==  2 ) THEN
    1420             :             GOTO 70
    1421         828 :          ELSEIF ( ido < 2 ) THEN
    1422             :             GOTO 100
    1423             :          ENDIF
    1424         510 :          idp2 = ido + 2
    1425        1572 :          DO 60 k = 1,l1
    1426        3204 :             DO 50 i = 3,ido,2
    1427        2142 :                ic = idp2 - i
    1428      535038 :                DO 40 m = 1,mp
    1429             :                   ch(m,i-1,k,1) = (cc(m,i-1,1,k)+cc(m,ic-1,4,k)) + &
    1430      532896 :                                   (cc(m,i-1,3,k)+cc(m,ic-1,2,k))
    1431             :                   ch(m,i,k,1) = (cc(m,i,1,k)-cc(m,ic,4,k)) + &
    1432      532896 :                                 (cc(m,i,3,k)-cc(m,ic,2,k))
    1433             :                   ch(m,i-1,k,2) = wa1(i-2)* ((cc(m,i-1,1,k)-cc(m,ic-1,4, &
    1434             :                                                                k))- (cc(m,i,3,k)+cc(m,ic,2,k))) - &
    1435             :                                   wa1(i-1)* ((cc(m,i,1,k)+cc(m,ic,4,k))+ &
    1436      532896 :                                              (cc(m,i-1,3,k)-cc(m,ic-1,2,k)))
    1437             :                   ch(m,i,k,2) = wa1(i-2)* ((cc(m,i,1,k)+cc(m,ic,4,k))+ &
    1438             :                                            (cc(m,i-1,3,k)-cc(m,ic-1,2,k))) + &
    1439             :                                 wa1(i-1)* ((cc(m,i-1,1,k)-cc(m,ic-1,4,k))- &
    1440      532896 :                                            (cc(m,i,3,k)+cc(m,ic,2,k)))
    1441             :                   ch(m,i-1,k,3) = wa2(i-2)* ((cc(m,i-1,1,k)+cc(m,ic-1,4, &
    1442             :                                                                k))- (cc(m,i-1,3,k)+cc(m,ic-1,2,k))) - &
    1443             :                                   wa2(i-1)* ((cc(m,i,1,k)-cc(m,ic,4,k))- &
    1444      532896 :                                              (cc(m,i,3,k)-cc(m,ic,2,k)))
    1445             :                   ch(m,i,k,3) = wa2(i-2)* ((cc(m,i,1,k)-cc(m,ic,4,k))- &
    1446             :                                            (cc(m,i,3,k)-cc(m,ic,2,k))) + &
    1447             :                                 wa2(i-1)* ((cc(m,i-1,1,k)+cc(m,ic-1,4,k))- &
    1448      532896 :                                            (cc(m,i-1,3,k)+cc(m,ic-1,2,k)))
    1449             :                   ch(m,i-1,k,4) = wa3(i-2)* ((cc(m,i-1,1,k)-cc(m,ic-1,4, &
    1450             :                                                                k))+ (cc(m,i,3,k)+cc(m,ic,2,k))) - &
    1451             :                                   wa3(i-1)* ((cc(m,i,1,k)+cc(m,ic,4,k))- &
    1452      532896 :                                              (cc(m,i-1,3,k)-cc(m,ic-1,2,k)))
    1453             :                   ch(m,i,k,4) = wa3(i-2)* ((cc(m,i,1,k)+cc(m,ic,4,k))- &
    1454             :                                            (cc(m,i-1,3,k)-cc(m,ic-1,2,k))) + &
    1455             :                                 wa3(i-1)* ((cc(m,i-1,1,k)-cc(m,ic-1,4,k))+ &
    1456      532896 :                                            (cc(m,i,3,k)+cc(m,ic,2,k)))
    1457        2142 : 40             END DO
    1458        1062 : 50          END DO
    1459         510 : 60       END DO
    1460         510 :          IF (mod(ido,2) == 1) RETURN
    1461             : 70       CONTINUE
    1462        2574 :          DO 90 k = 1,l1
    1463      251310 :             DO 80 m = 1,mp
    1464             :                ch(m,ido,k,1) = (cc(m,ido,1,k)+cc(m,ido,3,k)) + &
    1465      250272 :                                (cc(m,ido,1,k)+cc(m,ido,3,k))
    1466             :                ch(m,ido,k,2) = sqrt2* ((cc(m,ido,1,k)-cc(m,ido,3,k))- &
    1467      250272 :                                        (cc(m,1,2,k)+cc(m,1,4,k)))
    1468             :                ch(m,ido,k,3) = (cc(m,1,4,k)-cc(m,1,2,k)) + &
    1469      250272 :                                (cc(m,1,4,k)-cc(m,1,2,k))
    1470             :                ch(m,ido,k,4) = -sqrt2* ((cc(m,ido,1,k)-cc(m,ido,3,k))+ &
    1471      250272 :                                         (cc(m,1,2,k)+cc(m,1,4,k)))
    1472        1038 : 80          END DO
    1473         816 : 90       END DO
    1474             : 100      RETURN
    1475             :       END SUBROUTINE vradb4
    1476             : 
    1477          16 :       SUBROUTINE vradb5(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3,wa4)
    1478             : 
    1479             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1480             : 
    1481             :          USE m_constants, ONLY : pimach
    1482             :          IMPLICIT NONE
    1483             :          INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
    1484             :          REAL, INTENT(IN)    :: cc(mdimc,ido,5,l1), wa1(ido), wa2(ido), &
    1485             :                                 wa3(ido),wa4(ido)
    1486             :          REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,5)
    1487             : 
    1488             :          INTEGER :: i, k, m, ic, idp2
    1489             :          REAL :: arg, tr11, ti11, tr12, ti12
    1490             : 
    1491          16 :          arg = 2.*pimach()/5.
    1492          16 :          tr11 = cos(arg)
    1493          16 :          ti11 = sin(arg)
    1494          16 :          tr12 = cos(2.*arg)
    1495          16 :          ti12 = sin(2.*arg)
    1496         112 :          DO 20 k = 1,l1
    1497       31200 :             DO 10 m = 1,mp
    1498             :                ch(m,1,k,1) = cc(m,1,1,k) + 2.*cc(m,ido,2,k) + &
    1499       31104 :                &                     2.*cc(m,ido,4,k)
    1500             :                ch(m,1,k,2) = (cc(m,1,1,k)+tr11*2.*cc(m,ido,2,k)+ &
    1501             :                               tr12*2.*cc(m,ido,4,k)) - &
    1502       31104 :                              (ti11*2.*cc(m,1,3,k)+ti12*2.*cc(m,1,5,k))
    1503             :                ch(m,1,k,3) = (cc(m,1,1,k)+tr12*2.*cc(m,ido,2,k)+ &
    1504             :                               tr11*2.*cc(m,ido,4,k)) - &
    1505       31104 :                              (ti12*2.*cc(m,1,3,k)-ti11*2.*cc(m,1,5,k))
    1506             :                ch(m,1,k,4) = (cc(m,1,1,k)+tr12*2.*cc(m,ido,2,k)+ &
    1507             :                               tr11*2.*cc(m,ido,4,k)) + &
    1508       31104 :                              (ti12*2.*cc(m,1,3,k)-ti11*2.*cc(m,1,5,k))
    1509             :                ch(m,1,k,5) = (cc(m,1,1,k)+tr11*2.*cc(m,ido,2,k)+ &
    1510             :                               tr12*2.*cc(m,ido,4,k)) + &
    1511       31104 :                              (ti11*2.*cc(m,1,3,k)+ti12*2.*cc(m,1,5,k))
    1512          96 : 10          END DO
    1513          16 : 20       END DO
    1514          16 :          IF (ido == 1) RETURN
    1515           8 :          idp2 = ido + 2
    1516          24 :          DO 50 k = 1,l1
    1517          48 :             DO 40 i = 3,ido,2
    1518          32 :                ic = idp2 - i
    1519       10400 :                DO 30 m = 1,mp
    1520             :                   ch(m,i-1,k,1) = cc(m,i-1,1,k) + &
    1521             :                                   (cc(m,i-1,3,k)+cc(m,ic-1,2,k)) + &
    1522       10368 :                                   (cc(m,i-1,5,k)+cc(m,ic-1,4,k))
    1523             :                   ch(m,i,k,1) = cc(m,i,1,k) + (cc(m,i,3,k)-cc(m,ic,2,k)) + &
    1524       10368 :                                 (cc(m,i,5,k)-cc(m,ic,4,k))
    1525             :                   ch(m,i-1,k,2) = wa1(i-2)* ((cc(m,i-1,1,k)+tr11* (cc(m, &
    1526             :                   i-1,3,k)+cc(m,ic-1,2,k))+tr12* (cc(m,i-1, &
    1527             :                   &                          5,k)+cc(m,ic-1,4,k)))- &
    1528             :                   (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
    1529             :                   k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
    1530             :                   wa1(i-1)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
    1531             :                   k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
    1532             :                   ic,4,k)))+ (ti11* (cc(m,i-1,3,k)-cc(m, &
    1533             :                   ic-1,2,k))+ti12* (cc(m,i-1,5,k)-cc(m, &
    1534       10368 :                   ic-1,4,k))))
    1535             :                   ch(m,i,k,2) = wa1(i-2)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
    1536             :                   k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
    1537             :                   ic,4,k)))+ (ti11* (cc(m,i-1,3,k)-cc(m,ic-1, &
    1538             :                   &                        2,k))+ti12* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
    1539             :                   k)))) + wa1(i-1)* ((cc(m,i-1,1, &
    1540             :                   k)+tr11* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
    1541             :                   k))+tr12* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))- &
    1542             :                   (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
    1543       10368 :                   k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k))))
    1544             :                   ch(m,i-1,k,3) = wa2(i-2)* ((cc(m,i-1,1,k)+tr12* (cc(m, &
    1545             :                   i-1,3,k)+cc(m,ic-1,2,k))+tr11* (cc(m,i-1, &
    1546             :                   &                          5,k)+cc(m,ic-1,4,k)))- &
    1547             :                   (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
    1548             :                   k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
    1549             :                   wa2(i-1)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
    1550             :                   k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
    1551             :                   ic,4,k)))+ (ti12* (cc(m,i-1,3,k)-cc(m, &
    1552             :                   ic-1,2,k))-ti11* (cc(m,i-1,5,k)-cc(m, &
    1553       10368 :                   ic-1,4,k))))
    1554             :                   ch(m,i,k,3) = wa2(i-2)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
    1555             :                   k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
    1556             :                   ic,4,k)))+ (ti12* (cc(m,i-1,3,k)-cc(m,ic-1, &
    1557             :                   &                        2,k))-ti11* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
    1558             :                   k)))) + wa2(i-1)* ((cc(m,i-1,1, &
    1559             :                   k)+tr12* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
    1560             :                   k))+tr11* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))- &
    1561             :                   (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
    1562       10368 :                   k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k))))
    1563             :                   ch(m,i-1,k,4) = wa3(i-2)* ((cc(m,i-1,1,k)+tr12* (cc(m, &
    1564             :                   i-1,3,k)+cc(m,ic-1,2,k))+tr11* (cc(m,i-1, &
    1565             :                   &                          5,k)+cc(m,ic-1,4,k)))+ &
    1566             :                   (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
    1567             :                   k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
    1568             :                   wa3(i-1)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
    1569             :                   k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
    1570             :                   ic,4,k)))- (ti12* (cc(m,i-1,3,k)-cc(m, &
    1571             :                   ic-1,2,k))-ti11* (cc(m,i-1,5,k)-cc(m, &
    1572       10368 :                   ic-1,4,k))))
    1573             :                   ch(m,i,k,4) = wa3(i-2)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
    1574             :                   k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
    1575             :                   ic,4,k)))- (ti12* (cc(m,i-1,3,k)-cc(m,ic-1, &
    1576             :                   &                        2,k))-ti11* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
    1577             :                   k)))) + wa3(i-1)* ((cc(m,i-1,1, &
    1578             :                   k)+tr12* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
    1579             :                   k))+tr11* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+ &
    1580             :                   (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
    1581       10368 :                   k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k))))
    1582             :                   ch(m,i-1,k,5) = wa4(i-2)* ((cc(m,i-1,1,k)+tr11* (cc(m, &
    1583             :                   i-1,3,k)+cc(m,ic-1,2,k))+tr12* (cc(m,i-1, &
    1584             :                   &                          5,k)+cc(m,ic-1,4,k)))+ &
    1585             :                   (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
    1586             :                   k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
    1587             :                   wa4(i-1)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
    1588             :                   k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
    1589             :                   ic,4,k)))- (ti11* (cc(m,i-1,3,k)-cc(m, &
    1590             :                   ic-1,2,k))+ti12* (cc(m,i-1,5,k)-cc(m, &
    1591       10368 :                   ic-1,4,k))))
    1592             :                   ch(m,i,k,5) = wa4(i-2)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
    1593             :                   k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
    1594             :                   ic,4,k)))- (ti11* (cc(m,i-1,3,k)-cc(m,ic-1, &
    1595             :                   &                        2,k))+ti12* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
    1596             :                   k)))) + wa4(i-1)* ((cc(m,i-1,1, &
    1597             :                   k)+tr11* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
    1598             :                   k))+tr12* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+ &
    1599             :                   (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
    1600       10368 :                   k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k))))
    1601          32 : 30             END DO
    1602          16 : 40          END DO
    1603           8 : 50       END DO
    1604             :          RETURN
    1605             :       END SUBROUTINE vradb5
    1606             : 
    1607           0 :       SUBROUTINE vradbg(mp,ido,ip,l1,idl1,cc,c1,c2,ch,ch2, &
    1608           0 :                         mdimc,wa)
    1609             : 
    1610             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1611             : 
    1612             :          USE m_constants, ONLY : pimach
    1613             :          IMPLICIT NONE
    1614             :          INTEGER, INTENT(IN) :: mp, ido, ip, l1, idl1, mdimc
    1615             :          REAL, INTENT(IN)    :: cc(mdimc,ido,ip,l1), wa(ido)
    1616             :          REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,ip), ch2(mdimc,idl1,ip), &
    1617             :                                 c1(mdimc,ido,l1,ip), c2(mdimc,idl1,ip)
    1618             : 
    1619             :          INTEGER :: i, j, j2, k, l, m, ik, is, idp2, nbd, ipp2, ipph, &
    1620             :                     ic, jc, lc, idij
    1621             :          REAL :: ai1, ar1, ar2, ai2, ar1h, ar2h, tpi, arg, dcp, dsp, &
    1622             :                  dc2, ds2
    1623             : 
    1624           0 :          tpi = 2.*pimach()
    1625           0 :          arg = tpi/real(ip)
    1626           0 :          dcp = cos(arg)
    1627           0 :          dsp = sin(arg)
    1628           0 :          idp2 = ido + 2
    1629           0 :          nbd = (ido-1)/2
    1630           0 :          ipp2 = ip + 2
    1631           0 :          ipph = (ip+1)/2
    1632           0 :          IF (ido < l1) GO TO 40
    1633           0 :          DO 30 k = 1,l1
    1634           0 :             DO 20 i = 1,ido
    1635           0 :                DO 10 m = 1,mp
    1636           0 :                   ch(m,i,k,1) = cc(m,i,1,k)
    1637           0 : 10             END DO
    1638           0 : 20          END DO
    1639           0 : 30       END DO
    1640           0 :          GO TO 80
    1641           0 : 40       DO 70 i = 1,ido
    1642           0 :             DO 60 k = 1,l1
    1643           0 :                DO 50 m = 1,mp
    1644           0 :                   ch(m,i,k,1) = cc(m,i,1,k)
    1645           0 : 50             END DO
    1646           0 : 60          END DO
    1647           0 : 70       END DO
    1648           0 : 80       DO 110 j = 2,ipph
    1649           0 :             jc = ipp2 - j
    1650           0 :             j2 = j + j
    1651           0 :             DO 100 k = 1,l1
    1652           0 :                DO 90 m = 1,mp
    1653           0 :                   ch(m,1,k,j) = cc(m,ido,j2-2,k) + cc(m,ido,j2-2,k)
    1654           0 :                   ch(m,1,k,jc) = cc(m,1,j2-1,k) + cc(m,1,j2-1,k)
    1655           0 : 90             END DO
    1656           0 : 100         END DO
    1657           0 : 110      END DO
    1658           0 :          IF (ido == 1) GO TO 210
    1659           0 :          IF (nbd < l1) GO TO 160
    1660           0 :          DO 150 j = 2,ipph
    1661           0 :             jc = ipp2 - j
    1662           0 :             DO 140 k = 1,l1
    1663           0 :                DO 130 i = 3,ido,2
    1664           0 :                   ic = idp2 - i
    1665           0 :                   DO 120 m = 1,mp
    1666           0 :                      ch(m,i-1,k,j) = cc(m,i-1,2*j-1,k) + cc(m,ic-1,2*j-2,k)
    1667             :                      ch(m,i-1,k,jc) = cc(m,i-1,2*j-1,k) - &
    1668           0 :                                       cc(m,ic-1,2*j-2,k)
    1669           0 :                      ch(m,i,k,j) = cc(m,i,2*j-1,k) - cc(m,ic,2*j-2,k)
    1670           0 :                      ch(m,i,k,jc) = cc(m,i,2*j-1,k) + cc(m,ic,2*j-2,k)
    1671           0 : 120               END DO
    1672           0 : 130            END DO
    1673           0 : 140         END DO
    1674           0 : 150      END DO
    1675           0 :          GO TO 210
    1676           0 : 160      DO 200 j = 2,ipph
    1677           0 :             jc = ipp2 - j
    1678           0 :             DO 190 i = 3,ido,2
    1679           0 :                ic = idp2 - i
    1680           0 :                DO 180 k = 1,l1
    1681           0 :                   DO 170 m = 1,mp
    1682           0 :                      ch(m,i-1,k,j) = cc(m,i-1,2*j-1,k) + cc(m,ic-1,2*j-2,k)
    1683             :                      ch(m,i-1,k,jc) = cc(m,i-1,2*j-1,k) - &
    1684           0 :                                       cc(m,ic-1,2*j-2,k)
    1685           0 :                      ch(m,i,k,j) = cc(m,i,2*j-1,k) - cc(m,ic,2*j-2,k)
    1686           0 :                      ch(m,i,k,jc) = cc(m,i,2*j-1,k) + cc(m,ic,2*j-2,k)
    1687           0 : 170               END DO
    1688           0 : 180            END DO
    1689           0 : 190         END DO
    1690           0 : 200      END DO
    1691             : 210      ar1 = 1.
    1692             :          ai1 = 0.
    1693           0 :          DO 270 l = 2,ipph
    1694           0 :             lc = ipp2 - l
    1695           0 :             ar1h = dcp*ar1 - dsp*ai1
    1696           0 :             ai1 = dcp*ai1 + dsp*ar1
    1697           0 :             ar1 = ar1h
    1698           0 :             DO 230 ik = 1,idl1
    1699           0 :                DO 220 m = 1,mp
    1700           0 :                   c2(m,ik,l) = ch2(m,ik,1) + ar1*ch2(m,ik,2)
    1701           0 :                   c2(m,ik,lc) = ai1*ch2(m,ik,ip)
    1702           0 : 220            END DO
    1703           0 : 230         END DO
    1704             :             dc2 = ar1
    1705             :             ds2 = ai1
    1706             :             ar2 = ar1
    1707             :             ai2 = ai1
    1708           0 :             DO 260 j = 3,ipph
    1709           0 :                jc = ipp2 - j
    1710           0 :                ar2h = dc2*ar2 - ds2*ai2
    1711           0 :                ai2 = dc2*ai2 + ds2*ar2
    1712           0 :                ar2 = ar2h
    1713           0 :                DO 250 ik = 1,idl1
    1714           0 :                   DO 240 m = 1,mp
    1715           0 :                      c2(m,ik,l) = c2(m,ik,l) + ar2*ch2(m,ik,j)
    1716           0 :                      c2(m,ik,lc) = c2(m,ik,lc) + ai2*ch2(m,ik,jc)
    1717           0 : 240               END DO
    1718           0 : 250            END DO
    1719           0 : 260         END DO
    1720           0 : 270      END DO
    1721           0 :          DO 300 j = 2,ipph
    1722           0 :             DO 290 ik = 1,idl1
    1723           0 :                DO 280 m = 1,mp
    1724           0 :                   ch2(m,ik,1) = ch2(m,ik,1) + ch2(m,ik,j)
    1725           0 : 280            END DO
    1726           0 : 290         END DO
    1727           0 : 300      END DO
    1728           0 :          DO 330 j = 2,ipph
    1729           0 :             jc = ipp2 - j
    1730           0 :             DO 320 k = 1,l1
    1731           0 :                DO 310 m = 1,mp
    1732           0 :                   ch(m,1,k,j) = c1(m,1,k,j) - c1(m,1,k,jc)
    1733           0 :                   ch(m,1,k,jc) = c1(m,1,k,j) + c1(m,1,k,jc)
    1734           0 : 310            END DO
    1735           0 : 320         END DO
    1736           0 : 330      END DO
    1737           0 :          IF (ido == 1) GO TO 430
    1738           0 :          IF (nbd < l1) GO TO 380
    1739           0 :          DO 370 j = 2,ipph
    1740           0 :             jc = ipp2 - j
    1741           0 :             DO 360 k = 1,l1
    1742           0 :                DO 350 i = 3,ido,2
    1743           0 :                   DO 340 m = 1,mp
    1744           0 :                      ch(m,i-1,k,j) = c1(m,i-1,k,j) - c1(m,i,k,jc)
    1745           0 :                      ch(m,i-1,k,jc) = c1(m,i-1,k,j) + c1(m,i,k,jc)
    1746           0 :                      ch(m,i,k,j) = c1(m,i,k,j) + c1(m,i-1,k,jc)
    1747           0 :                      ch(m,i,k,jc) = c1(m,i,k,j) - c1(m,i-1,k,jc)
    1748           0 : 340               END DO
    1749           0 : 350            END DO
    1750           0 : 360         END DO
    1751           0 : 370      END DO
    1752           0 :          GO TO 430
    1753           0 : 380      DO 420 j = 2,ipph
    1754           0 :             jc = ipp2 - j
    1755           0 :             DO 410 i = 3,ido,2
    1756           0 :                DO 400 k = 1,l1
    1757           0 :                   DO 390 m = 1,mp
    1758           0 :                      ch(m,i-1,k,j) = c1(m,i-1,k,j) - c1(m,i,k,jc)
    1759           0 :                      ch(m,i-1,k,jc) = c1(m,i-1,k,j) + c1(m,i,k,jc)
    1760           0 :                      ch(m,i,k,j) = c1(m,i,k,j) + c1(m,i-1,k,jc)
    1761           0 :                      ch(m,i,k,jc) = c1(m,i,k,j) - c1(m,i-1,k,jc)
    1762           0 : 390               END DO
    1763           0 : 400            END DO
    1764           0 : 410         END DO
    1765           0 : 420      END DO
    1766             : 430      CONTINUE
    1767           0 :          IF (ido == 1) RETURN
    1768           0 :          DO 450 ik = 1,idl1
    1769           0 :             DO 440 m = 1,mp
    1770           0 :                c2(m,ik,1) = ch2(m,ik,1)
    1771           0 : 440         END DO
    1772           0 : 450      END DO
    1773           0 :          DO 480 j = 2,ip
    1774           0 :             DO 470 k = 1,l1
    1775           0 :                DO 460 m = 1,mp
    1776           0 :                   c1(m,1,k,j) = ch(m,1,k,j)
    1777           0 : 460            END DO
    1778           0 : 470         END DO
    1779           0 : 480      END DO
    1780           0 :          IF (nbd > l1) GO TO 530
    1781           0 :          is = -ido
    1782           0 :          DO 520 j = 2,ip
    1783           0 :             is = is + ido
    1784           0 :             idij = is
    1785           0 :             DO 510 i = 3,ido,2
    1786           0 :                idij = idij + 2
    1787           0 :                DO 500 k = 1,l1
    1788           0 :                   DO 490 m = 1,mp
    1789             :                      c1(m,i-1,k,j) = wa(idij-1)*ch(m,i-1,k,j) - &
    1790           0 :                                      wa(idij)*ch(m,i,k,j)
    1791             :                      c1(m,i,k,j) = wa(idij-1)*ch(m,i,k,j) + &
    1792           0 :                                    wa(idij)*ch(m,i-1,k,j)
    1793           0 : 490               END DO
    1794           0 : 500            END DO
    1795           0 : 510         END DO
    1796           0 : 520      END DO
    1797           0 :          GO TO 580
    1798           0 : 530      is = -ido
    1799           0 :          DO 570 j = 2,ip
    1800           0 :             is = is + ido
    1801           0 :             DO 560 k = 1,l1
    1802           0 :                idij = is
    1803           0 :                DO 550 i = 3,ido,2
    1804           0 :                   idij = idij + 2
    1805           0 :                   DO 540 m = 1,mp
    1806             :                      c1(m,i-1,k,j) = wa(idij-1)*ch(m,i-1,k,j) - &
    1807           0 :                                      wa(idij)*ch(m,i,k,j)
    1808             :                      c1(m,i,k,j) = wa(idij-1)*ch(m,i,k,j) + &
    1809           0 :                                    wa(idij)*ch(m,i-1,k,j)
    1810           0 : 540               END DO
    1811           0 : 550            END DO
    1812           0 : 560         END DO
    1813           0 : 570      END DO
    1814             : 580      RETURN
    1815             :       END SUBROUTINE vradbg
    1816             : 
    1817         350 :       SUBROUTINE vrftb1(m,n,c,ch,mdimc,wa,fac)
    1818             : 
    1819             : !     VRFFTPK, VERSION 1, AUGUST 1985
    1820             : 
    1821             :          IMPLICIT NONE
    1822             :          INTEGER, INTENT(IN) :: m, n, mdimc
    1823             :          REAL, INTENT(IN)    :: wa(n), fac(15)
    1824             :          REAL, INTENT(INOUT) :: c(mdimc,n), ch(mdimc,n)
    1825             : 
    1826             :          INTEGER :: i, j, nf, na, l1, iw, k1, ip, ipo, idl1, ix2, ix3, ix4
    1827             :          INTEGER :: ido, l2
    1828             :          REAL    :: scale
    1829             : 
    1830         350 :          nf = fac(2)
    1831         350 :          na = 0
    1832         350 :          l1 = 1
    1833         350 :          iw = 1
    1834        1262 :          DO 160 k1 = 1,nf
    1835         912 :             ip = fac(k1+2)
    1836         912 :             l2 = ip*l1
    1837         912 :             ido = n/l2
    1838         912 :             idl1 = ido*l1
    1839         912 :             IF (ip /= 4) GO TO 30
    1840         828 :             ix2 = iw + ido
    1841         828 :             ix3 = ix2 + ido
    1842         828 :             IF (na /= 0) GO TO 10
    1843         498 :             CALL vradb4(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3))
    1844         498 :             GO TO 20
    1845         330 : 10          CALL vradb4(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3))
    1846         828 : 20          na = 1 - na
    1847         828 :             GO TO 150
    1848          84 : 30          IF (ip /= 2) GO TO 60
    1849          32 :             IF (na /= 0) GO TO 40
    1850          32 :             CALL vradb2(m,ido,l1,c,ch,mdimc,wa(iw))
    1851          32 :             GO TO 50
    1852           0 : 40          CALL vradb2(m,ido,l1,ch,c,mdimc,wa(iw))
    1853          32 : 50          na = 1 - na
    1854          32 :             GO TO 150
    1855          52 : 60          IF (ip /= 3) GO TO 90
    1856          36 :             ix2 = iw + ido
    1857          36 :             IF (na /= 0) GO TO 70
    1858          24 :             CALL vradb3(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2))
    1859          24 :             GO TO 80
    1860          12 : 70          CALL vradb3(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2))
    1861          36 : 80          na = 1 - na
    1862          36 :             GO TO 150
    1863          16 : 90          IF (ip /= 5) GO TO 120
    1864          16 :             ix2 = iw + ido
    1865          16 :             ix3 = ix2 + ido
    1866          16 :             ix4 = ix3 + ido
    1867          16 :             IF (na /= 0) GO TO 100
    1868           8 :             CALL vradb5(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    1869           8 :             GO TO 110
    1870           8 : 100         CALL vradb5(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    1871          16 : 110         na = 1 - na
    1872          16 :             GO TO 150
    1873           0 : 120         IF (na /= 0) GO TO 130
    1874           0 :             CALL vradbg(m,ido,ip,l1,idl1,c,c,c,ch,ch,mdimc,wa(iw))
    1875           0 :             GO TO 140
    1876           0 : 130         CALL vradbg(m,ido,ip,l1,idl1,ch,ch,ch,c,c,mdimc,wa(iw))
    1877           0 : 140         IF (ido == 1) na = 1 - na
    1878         912 : 150         l1 = l2
    1879         912 :             iw = iw + (ip-1)*ido
    1880         350 : 160      END DO
    1881         350 :          scale = sqrt(1./n)
    1882         350 :          IF (na == 0) GO TO 190
    1883       25060 :          DO 180 j = 1,n
    1884     3234856 :             DO 170 i = 1,m
    1885     3222432 :                c(i,j) = scale*ch(i,j)
    1886       12424 : 170         END DO
    1887         212 : 180      END DO
    1888         138 :          RETURN
    1889        4554 : 190      DO 210 j = 1,n
    1890      320160 :             DO 200 i = 1,m
    1891      317952 :                c(i,j) = scale*c(i,j)
    1892        2208 : 200         END DO
    1893         138 : 210      END DO
    1894             :          RETURN
    1895             :       END SUBROUTINE vrftb1
    1896             : 
    1897             :    END SUBROUTINE vrfftb
    1898             : 
    1899             : END MODULE m_rfft

Generated by: LCOV version 1.13