Date:    Thu, 23 Apr 1987 at 17:23 DNT (Copenhagen, GMT+2:00)
From: "Ole Holm Nielsen, NORDITA, Copenhagen" <OHNIELSE%DKNBI01.BITNET@wiscvm.wisc.edu>

Subject: 2,3-D FFT from Comput. Phys. Commun.
Ole H. Nielsen
Nordic Institute of Theoretical Physics, Copenhagen
E-mail: ohnielsen@nbivax.UUCP
        ohnielse @ dknbi01.BITNET

Cut here:
----------------------------------------------------------------------
C C.P.C. CODE: AALL
C MFFT: A PACKAGE FOR TWO- AND THREE-DIMENSIONAL VECTORIZED
C DISCRETE FOURIER TRANSFORMS.  A. NOBILE, V. ROBERTO.
C REF. IN COMP. PHYS. COMMUN. 42 (1986) 233
C
C***********************************************************************
C
C                           MFFT
C
C                  Version  1.1 - March 1986
C
C***********************************************************************
C
C**AUTHORS:
C                         A. NOBILE
C
C  INTERNATIONAL SCHOOL FOR ADVANCED STUDIES - TRIESTE, ITALY
C
C                             AND
C
C                         V. ROBERTO
C
C  DEPARTMENT OF THEORETICAL PHYSICS, UNIVERSITY OF TRIESTE, ITALY
C
C
C**PURPOSE:
C
C  MFFT IS A PACKAGE OF FORTRAN 77 SUBPROGRAMS TO PERFORM TWO-
C  AND THREE-DIMENSIONAL DISCRETE FOURIER TRANSFORMS (DFT).
C  THE DFT'S ARE CALCULATED IMPLEMENTING FAST FOURIER TRANSFORM (FFT)
C  ALGORITHMS IN VECTOR FORM.  THE PACKAGE IS SPECIALLY DESIGNED
C  FOR HIGH PERFORMANCES ON CRAY X-MP MACHINES.
C
C
C**USER INTERFACES:
C
C  THE PRESENT VERSION CONTAINS THE FOLLOWING USER ACCESSABLE ROUTINES:
C
C  1.  C2FFT  : 2-DIMENSIONAL TRANSFORM OF NX*NY COMPLEX DATA MATRICES
C  2.  C3FFT  : 3-DIMENSIONAL TRANSFORM OF NX*NY*NZ COMPL. DATA MATRICES
C  3.  R2FFT  : 2-DIMENSIONAL TRANSFORM OF NX*NY REAL DATA MATRICES
C  4.  R3FFT  : 3-DIMENSIONAL TRANSFORM OF NX*NY*NZ REAL DATA MATRICES
C
C
C**DATA SIZE:
C
C  EACH DIMENSION IN THE INPUT MATRICES HAS TO BE A PRODUCT OF INTEGER
C  POWERS OF 2,3,5 (SEE THE LONG WRITE-UP, REF.[2]).
C  IN REAL TRANSFORMS, THE FIRST DIMENSION OF DATA ARRAYS HAS TO BE AN
C  EVEN NUMBER GREATER OR EQUAL TO N1+2 (WHERE N1 IS THE ACTUAL SIZE OF
C  DATA, FIRST DIMENSION).
C
C
C
C**FURTHER DETAILS:
C
C  THE THEORETICAL GROUNDS ON WHICH MFFT IS BASED, AS WELL AS THE
C  DETAILS OF ITS IMPLEMENTATION ON A CRAY X-MP CAN BE FOUND IN REF.[1]
C  THE LONG WRITE-UP IS REPORTED IN REF.[2].
C  BOTH REFERENCES CONTAIN PERFORMANCE EVALUATIONS (CPU-TIME, OPERATION
C  COUNTS, MEGAFLOP RATES).
C
C
C**REFERENCES:
C
C  [1] A.NOBILE AND V.ROBERTO:'EFFICIENT IMPLEMENTATION OF MULTIDIMEN-
C      SIONAL FAST FOURIER TRANSFORMS ON A CRAY X-MP'-  PREPRINT
C      SISSA 82/85, TO APPEAR IN COMPUTER PHYSICS COMMUNICATIONS
C  [2] A.NOBILE AND V.ROBERTO: 'MFFT : A PACKAGE FOR TWO- AND
C      THREE-DIMENSIONAL VECTORIZED DISCRETE FOURIER TRANSFORMS'
C      SISSA TRIESTE 1986 PREPRINT, SUBMITTED FOR PUBLICATION TO
C      COMPUTER PHYSICS COMMUNICATIONS
C
C
C
C**REQUESTS:
C
C  FOR ANY PROBLEM CONCERNING MFFT, PLEASE CONTACT:
C
C  DR. V.ROBERTO, DEPT. OF THEORETICAL PHYSICS
C                 STRADA COSTIERA, 11
C                 I - 34014  TRIESTE   (ITALY)
C
C  TELEX NO.: 460932 ICTP I
C
C  MESSAGES THROUGH EARNET TO: 064 AT ITSSISSA
C
C
C***********************************************************************
C
C
      SUBROUTINE C2FFT(C,ID,NM,NN,WM,WN,ISIG,IORD,IWORK,IERR)
C***PURPOSE:
C       THIS ROUTINE PERFORMS A 2-DIMENSIONAL COMPLEX FOURIER TRANSFORM,
C       OF ORDER NM*NN.
C***USAGE:
C       THE USER IS EXPECTED TO PROVIDE THE DATA IN A 2-DIMENSIONAL
C       COMPLEX ARRAY C, DIMENSIONED IN THE CALLING PROGRAM C(ID,NN);
C       ID CAN BE DIFFERENT FROM NM, AND IT IS RECOMMENDED THAT IT IS
C       CHOSEN EQUAL TO NM+1 IF NM IS EVEN OR TO NM IF NM IS ODD.
C       THE ELEMENTS C(K,*),NM<K<=ID MUST BE ZEROED. THE ROUTINE IS
C       INTENDED FOR REPEATED USAGE, THUS SEPARATE SET-UP AND
C       OPERATING CALLS ARE AVAILABLE : THE USER SHOULD ALWAYS PERFORM
C       A SET-UP CALL ( ISIG=0 ) PASSING THE CHOSEN PARAMETERS, BEFORE
C       PERFORMING THE ACTUAL TRANSFORM ( ISIG= +1 OR -1 ); THE USER CAN
C       CHOOSE WHETHER TO OBTAIN THE RESULTS OF THE DIRECT TRANSFORM
C       IN NATURAL ORDER (ISIG=-1,IORD=1) OR LEAVE THEM IN THE
C      'BIT-REVERSED' ORDER( ISIG=-1,IORD=0); THIS CHOICE SAVES
C       SOME COMPUTER TIME, AND IT IS RECOMMENDED IN CASES DISCUSSED
C       IN THE LONG WRITE-UP. ANALOGOUSLY,THE INVERSE TRANSFORM ACCEPTS
C       INPUT ( PLEASE NOTE! ) DATA IN NATURAL ORDER ( ISIG=1,IORD=1),
C       OR DATA ALREADY SUBJECTED TO A BIT-REVERSAL PERMUTATION( ISIG=1
C       IORD=0 ).
C***PARAMETERS :
C       INPUT :
C       C : ARRAY TO BE TRANSFORMED; DECLARED COMPLEX C(ID,NN) IN THE
C           CALLING PROGRAM;
C       ID : FIRST DIMENSION OF C IN THE CALLING PROGRAM
C       ISIG : OPTION FLAG : ISIG=0 : SET-UP RUN, C NOT USED
C                            ISIG=-1: DIRECT TRANSFORM
C                            ISIG=+1: INVERSE TRANSFORM
C       WM,WN : INTEGER ARRAYS , USED TO HOST TABLES FOR THE TRANSFORM;
C               DIMENSIONED IN THE CALLING PROGRAM AT LEAST (4*NM+14)
C               AND (4*NN+14) RESPECTIVELY; IF ISIG.NE.0 ,THEY ARE
C               ASSUMED TO HAVE BEEN SET BY A PREVIOUS CALL WITH ISIG=0
C               (ALL OTHER ARGUMENTS UNCHANGED), AND NEVER HAVE BEEN
C               MODIFIED SINCE THEN;
C               IF NM=NN, WM AND WN DO NOT NEED TO BE DISTINCT;
C       NM : ORDER OF THE TRANSFORM ALONG THE COLUMNS OF C
C       NN : ORDER OF THE TRANSFORM ALONG THE ROWS OF C
C       IORD : OPTION FLAG : =1 : OUTPUT IN NATURAL ORDER (ISIG=-1)
C                                 INPUT IN NATURAL ORDER  (ISIG=+1)
C                            =0 : OUTPUT IN BIT-REVERSED ORDER(ISIG=-1)
C                                 INPUT IN BIT-REVERSED ORDER(ISIG=+1)
C       IWORK : INTEGER ARRAY, USED AS WORK AREA FOR REORDERING IF IORD=
C                1, IT MUST BE AT LEAST MAX(NM,NN) WORDS LONG.
C                IT IS UNUSED IF IORD=0.
C
C***OUTPUT :
C       C : TRANSFORMED ARRAY
C       WM, WN : ONLY IF ISIG=0, WM AND WN ARE FILLED WITH THE
C               APPROPRIATE  TABLES (TWIDDLE FACTORS).
C       IWORK : UNDEFINED
C       IERR  : ERROR CODE : =0 : SUCCESSFUL
C                            =1 : DATA DIMENSIONS ARE NOT CORRECT
C                            =2 : PRIME FACTORS DIFFERENT FROM 2,3,5
C                                 ARE PRESENT IN DATA DIMENSIONS
C                            =3 : TABLES NOT CORRECTLY INITIALIZED
C
      COMPLEX C(*)
      INTEGER WM(-14:*),WN(-14:*)
      INTEGER IWORK(*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      IF(ID.LT.NM)THEN
        IERR=IDERR
        RETURN
      ENDIF
      IERR=0
      IF(ISIG.EQ.0)THEN
        CALL MFFTP(NM,WM,0,IERR)
        IF(IERR.NE.0)RETURN
        IF(NN.NE.NM)THEN
          CALL MFFTP(NN,WN,0,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTZ0(WM,1,2*NM+7,WN,1)
        ENDIF
        RETURN

      ELSE IF(ISIG.GT.0)THEN

        IF(IORD.NE.0)THEN
         CALL MFFTOV(C,1,ID,NM,NN,WM(3*NM),IWORK)
         CALL MFFTOV(C,ID,1,NN,NM,WN(3*NN),IWORK)
        ENDIF

        CALL MFFTIV(C,1,ID,NM,NN,WM,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTIV(C,ID,1,NN,NM,WN,IERR)
        IF(IERR.NE.0)RETURN

      ELSE
        CALL MFFTDV(C,1,ID,NM,NN,WM,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTDV(C,ID,1,NN,NM,WN,IERR)
        IF(IERR.NE.0)RETURN

        IF(IORD.NE.0)THEN
         CALL MFFTOV(C,1,ID,NM,NN,WM(2*NM),IWORK)
         CALL MFFTOV(C,ID,1,NN,NM,WN(2*NN),IWORK)
        ENDIF
      ENDIF

      END
      SUBROUTINE C3FFT(C,ID,NL,NM,NN,WL,WM,WN,IOPT,ISIG,IORD,IWORK,IERR)

C***PURPOSE:
C       THIS ROUTINE PERFORMS A 3-DIMENSIONAL COMPLEX FOURIER TRANSFORM,
C       OF ORDER NL*NM*NN  .
C       IT IS A USER INTERFACE FOR A PACKAGE (MFFT), SPECIALLY DESIGNED
C       FOR HIGH PERFORMANCE ON CRAY X-MP MACHINES.
C***USAGE:
C       THE USER IS EXPECTED TO PROVIDE THE DATA IN A 3-DIMENSIONAL
C       COMPLEX ARRAY C, DIMENSIONED IN THE CALLING PROGRAM C(ID,NM,NN);
C       ID CAN BE DIFFERENT FROM NL, AND IT IS RECOMMENDED THAT IT IS
C       CHOSEN ODD FOR MAXIMUM PERFORMANCE. THE ROUTINE IS
C       INTENDED FOR REPEATED USAGE, THUS SEPARATE SET-UP AND
C       OPERATING CALLS ARE AVAILABLE : THE USER SHOULD ALWAYS PERFORM
C       A SET-UP CALL ( ISIG=0 ) PASSING THE CHOSEN PARAMETERS, BEFORE
C       PERFORMING THE ACTUAL TRANSFORM ( ISIG= +1 OR -1 ); THE USER CAN
C       CHOOSE WHETHER TO OBTAIN THE RESULTS OF THE DIRECT TRANSFORM
C       IN NATURAL ORDER (ISIG=-1,IORD=1) OR LEAVE THEM IN THE
C       BIT-REVERSED  ORDER( ISIG=-1,IORD=0); THIS CHOICE SAVES
C       SOME COMPUTER TIME, AND IT IS RECOMMENDED IN CASES DISCUSSED
C       IN THE LONG WRITE-UP. ANALOGOUSLY, THE INVERSE TRANSFORM ACCEPTS
C       INPUT ( PLEASE NOTE! ) DATA IN NATURAL ORDER ( ISIG=1,IORD=1),
C       OR DATA ALREADY SUBJECTED TO A BIT-REVERSAL PERMUTATION( ISIG=1
C       IORD=0).
C       A SPECIAL TREATMENT IS AVAILABLE TO SPEED UP THE TRANSFORM OF
C       SMALL MATRICES. THIS TREATMENT IS ACTIVATED BY THE FLAG IOPT. IN
C       THIS CASE THE TABLES FOR THE SECOND DIMENSION ( WM ) ARE LARGER,
C       BUT THE INCREASE IN PERFORMANCE IS SUBSTANTIAL WHEN NM<32.
C***ARGUMENTS :
C       INPUT :
C       C : ARRAY TO BE TRANSFORMED; DECLARED COMPLEX C(ID,NM,NN) IN THE
C           CALLING PROGRAM;
C       ID : FIRST DIMENSION OF C IN THE CALLING PROGRAM
C       ISIG : OPTION FLAG : ISIG=0 : SET-UP RUN, C NOT USED
C                            ISIG=-1: DIRECT TRANSFORM
C                            ISIG=+1: INVERSE TRANSFORM
C       WL,WM,WN : INTEGER ARRAYS,USED TO HOST TABLES FOR THE TRANSFORMS
C               DIMENSIONED IN THE CALLING PROGRAM AT LEAST (4*NL+14)
C               (4*NM+14) AND (4*NN+14) RESPECTIVELY; IF IOPT=1
C               WM MUST BE DIMENSIONED AT LEAST (4*NM*(ID+1)+14)
C               IF ISIG.NE.0, THEY ARE ASSUMED TO HAVE BEEN SET BY A
C               PREVIOUS CALL WITH ISIG=0 AND OTHER ARGUMENTS EQUAL, AND
C               NEVER HAVE BEEN MODIFIED ;
C               WHEN THE CORRESPONDING ORDERS ARE EQUAL, THEY DO NOT
C               NEED TO BE DISTINCT.
C       NL : ORDER OF THE TRANSFORM ALONG THE COLUMNS OF C
C       NM : ORDER OF THE TRANSFORM ALONG THE ROWS OF C
C       NN : ORDER OF THE TRANSFORM ALONG THE THIRD DIMENSION OF C
C       IOPT : OPTION FLAG : =0 : NORMAL TREATMENT
C                            =1 : SPECIAL TREATMENT FOR IMPROVING
C                                 VECTORIZATION ON MATRICES WITH
C                                 SMALL NL; REQUIRES LONG WM(SEE);IF
C                                 REQUESTED, MUST BE PRESENT IN BOTH
C                                 THE SET-UP AND TRANSFORM CALLS;
C       IORD : OPTION FLAG : =1 : OUTPUT IN NATURAL ORDER (ISIG=-1)
C                                 INPUT IN NATURAL ORDER  (ISIG=+1)
C                            =0 : OUTPUT IN BIT-REVERSED ORDER(ISIG=-1)
C                                 INPUT IN BIT-REVERSED ORDER(ISIG=+1)
C       IWORK : INTEGER ARRAY, USED AS WORK AREA FOR REORDERING IF
C               IORD=1; IT MUST BE AT LEAST MAX(NL,NM,NN) WORDS LONG.
C
C        OUTPUT :
C       C : TRANSFORMED ARRAY
C       WL, WM, WN : ONLY IF ISIG=0, WL,WM AND WN ARE FILLED WITH THE
C                     APPROPRIATE TABLES
C       IWORK : UNDEFINED
C       IERR  : ERROR CODE : =0 : SUCCESSFUL
C                          : =1 : DATA DIMENSIONS ARE NOT CORRECT
C                            =2 : PRIME FACTORS DIFFERENT FROM 2,3,5
C                                 ARE PRESENT IN DATA DIMENSIONS
C                            =3 : TABLES NOT CORRECTLY INITIALIZED

      COMPLEX C(0:*)
      INTEGER WM(-14:*),WN(-14:*),WL(-14:*)
      INTEGER IWORK(*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)

      IF(ID.LT.NL)THEN
        IERR=IDERR
        RETURN
      ENDIF
      IERR=0
      IF(ISIG.EQ.0)THEN
        CALL MFFTP(NM,WM,ID*IOPT,IERR)
        IF(IERR.NE.0)RETURN
        IF(NM.NE.NN)THEN
          CALL MFFTP(NN,WN,0,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTZ0(WM,1,2*NM+7,WN,1)
        ENDIF
        IF(NL.EQ.NM)THEN
          CALL MFFTZ0(WM,1,2*NM+7,WL,1)
        ELSE IF(NL.EQ.NN)THEN
          CALL MFFTZ0(WN,1,2*NN+7,WL,1)
        ELSE
          CALL MFFTP(NL,WL,0,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
        RETURN

      ELSE   IF(ISIG.GT.0)THEN

        IF(IORD.NE.0)THEN
          CALL MFFTOV(C,1,ID,NL,NM*NN,WL(3*NL),IWORK)
          CALL MFFTOV(C,ID*NM,1,NN,ID*NM,WN(3*NN),IWORK)
          CALL MFFTOM(C,ID,ID*NM,1,NM,NN,NL,WM(3*NM),
     $                IWORK)
        ENDIF
        CALL MFFTIV(C,1,ID,NL,NM*NN,WL,IERR)
        IF(IERR.NE.0)RETURN
        IF(IOPT.EQ.0)THEN
          CALL MFFTIM(C,ID,ID*NM,1,NM,NN,NL,WM,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTIS(C,ID,ID*NM,1,NM,NN,NL,WM,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
        CALL MFFTIV(C,ID*NM,1,NN,ID*NM,WN,IERR)
        IF(IERR.NE.0)RETURN

      ELSE

        CALL MFFTDV(C,1,ID,NL,NM*NN,WL,IERR)
        IF(IERR.NE.0)RETURN
        IF(IOPT.EQ.0)THEN
          CALL MFFTDM(C,ID,ID*NM,1,NM,NN,NL,WM,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTDS(C,ID,ID*NM,1,NM,NN,NL,WM,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
        CALL MFFTDV(C,ID*NM,1,NN,ID*NM,WN,IERR)
        IF(IERR.NE.0)RETURN
        IF(IORD.NE.0)THEN
          CALL MFFTOV(C,1,ID,NL,NM*NN,WL(2*NL),IWORK)
          CALL MFFTOV(C,ID*NM,1,NN,ID*NM,WN(2*NN),IWORK)
          CALL MFFTOM(C,ID,ID*NM,1,NM,NN,NL,WM(2*NM),
     $                IWORK)
        ENDIF

      ENDIF

      END
       SUBROUTINE R2FFT(C,ID,NM,NN,WM,WN,ISIG,IORD,IWORK,IERR)
C
C***PURPOSE:
C       THIS ROUTINE PERFORMS A 2-DIMENSIONAL REAL FOURIER TRANSFORM,
C       OF ORDER NM*NN.
C***USAGE:
C       THE USER IS EXPECTED TO PROVIDE THE DATA IN A 2-DIMENSIONAL
C       REAL ARRAY C, DIMENSIONED IN THE CALLING PROGRAM C(ID,NN)
C       AND EQUIVALENCED TO A COMPLEX ARRAY CCOM(ID/2,NN).
C       WHEN A DIRECT TRASNFORM IS PERFORMED, HALF OF THE OUTPUT
C       COEFFICIENTS ARE CONTAINED IN CCOM(I,J) AS FOLLOWS
C        X(I,J)= CCOM(I,J) WITH I=1,(NM/2+1) AND J=1,NN
C       THE REMAINING ONES ARE OBTAINED EXPLOITING THE CONJUGATED
C       SYMMETRY RELATION
C        X(I,J) = CONJG(CCOM(NM-I+2,NN-J+2))
C       WITH I=NM/2+2,NM  AND  J=1,NN
C       NM (ACTUAL FIRST DIMENSION OF DATA) MUST BE AN EVEN NUMBER;
C       ID (DECLARED FIRST DIMENSION) MUST BE GR. OR EQ. TO NM+2.
C       THE ELEMENTS C(K,*), NM<K<=ID MUST BE ZEROED. THE ROUTINE IS
C       INTENDED FOR REPEATED USAGE, THUS SEPARATE SET-UP AND
C       OPERATING CALLS ARE AVAILABLE : THE USER SHOULD ALWAYS PERFORM
C       A SET-UP CALL ( ISIG=0 ) PASSING THE CHOSEN PARAMETERS, BEFORE
C       PERFORMING THE ACTUAL TRANSFORM ( ISIG= +1 OR -1 ); THE USER CAN
C       CHOOSE WHETHER TO OBTAIN THE RESULTS OF THE DIRECT TRANSFORM
C       IN NATURAL ORDER (ISIG=-1,IORD=1) OR LEAVE THEM IN THE
C       BIT-REVERSED  ORDER( ISIG=-1,IORD=0); THIS CHOICE SAVES
C       SOME COMPUTER TIME, AND IT IS RECOMMENDED IN CASES DISCUSSED
C       IN THE LONG WRITE-UP. ANALOGOUSLY, THE INVERSE TRANSFORM ACCEPTS
C       INPUT ( PLEASE NOTE! ) DATA IN NATURAL ORDER ( ISIG=1,IORD=1),
C       OR DATA ALREADY SUBJECTED TO A BIT-REVERSAL PERMUTATION( ISIG=1
C       IORD=0 ).
C***PARAMETERS :
C       INPUT :
C       C : ARRAY TO BE TRANSFORMED; DECLARED REAL C(ID,NN) IN THE
C           CALLING PROGRAM;
C       ID : FIRST DIMENSION OF C IN THE CALLING PROGRAM
C            IT MUST BE AN EVEN NUMBER EQUAL TO NM+2.
C       ISIG : OPTION FLAG : ISIG=0 : SET-UP RUN, C NOT USED
C                            ISIG=-1: DIRECT TRANSFORM
C                            ISIG=+1: INVERSE TRANSFORM
C       WM,WN : INTEGER ARRAYS , USED TO HOST TABLES FOR THE TRANSFORM;
C               DIMENSIONED IN THE CALLING PROGRAM IN THIS WAY :
C               WM : AT LEAST (6*NM+14)
C               WN : AT LEAST (4*NN+14)
C               IF ISIG.NE.0 ,THEY ARE ASSUMED TO HAVE BEEN SET
C                BY A PREVIOUS CALL WITH ISIG=0 AND
C               ALL THE OTHER ARGUMENTS UNCHANGED, AND NEVER HAVE BEEN
C               MODIFIED SINCE THEN;
C               IF NM=NN, WM AND WN DO NOT NEED TO BE DISTINCT;
C       NM : ORDER OF THE TRANSFORM ALONG THE COLUMNS OF C
C       NN : ORDER OF THE TRANSFORM ALONG THE ROWS OF C
C       IORD : OPTION FLAG : =1 : OUTPUT IN NATURAL ORDER (ISIG=-1)
C                                 INPUT IN NATURAL ORDER  (ISIG=+1)
C                            =0 : OUTPUT IN BIT-REVERSED ORDER(ISIG=-1)
C                                 INPUT IN BIT-REVERSED ORDER(ISIG=+1)
C                                 WARNING:THE FIRST DIMENSION IS ORDERED
C                                 IN ANY CASE BECAUSE OF POST-PROCESSING
C                                 (FOR DIRECT) AND PRE-PROCESSING FOR
C                                 INVERSE REAL TRANSFORMS
C       IWORK : INTEGER ARRAY, USED AS WORK AREA FOR REORDERING IF IORD=
C               1, MUST BE AT LEAST MAX(NM,NN) WORDS LONG.
C
C***OUTPUT :
C       C : TRANSFORMED ARRAY
C       WM, WN : ONLY IF ISIG=0, WM AND WN ARE FILLED WITH THE
C               APPROPRIATE  TABLES
C       IWORK : UNDEFINED
C       IERR  : ERROR CODE : =0 : SUCCESSFUL
C                            =1 : WRONG ID PARAMETER
C                            =2 : PRIME FACTORS DIFFERENT FROM 2,3,5
C                                 ARE PRESENT IN DATA DIMENSIONS
C                            =3 : TABLES NOT CORRECTLY INITIALIZED
C                            =4 : FIRST DIMENSION IS AN ODD NUMBER
      COMPLEX C(*)
      INTEGER WM(-14:*),WN(-14:*)
      INTEGER IWORK(*)
C
      INTEGER IDERR,FACERR,TBERR,ODDERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3,ODDERR=4)
C
      IF(ID.LT.NM+2)THEN
        IERR=IDERR
        RETURN
      ENDIF
      NM1=NM/2
      IF(NM1*2.NE.NM)THEN
        IERR=ODDERR
        RETURN
      ENDIF
      IERR=0
C
      IF(ISIG.EQ.0) THEN
        CALL MFFTRP(NM,WM(4*NM))
        CALL MFFTP(NM1,WM,0,IERR)
        IF(IERR.NE.0)RETURN
        IF(NN.NE.NM1) THEN
          CALL MFFTP(NN,WN,0,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTZ0(WM,1,2*NM1+7,WN,1)
        ENDIF
        RETURN
      ELSE  IF(ISIG.GT.0) THEN

        IF(IORD.NE.0) THEN
          CALL MFFTOV(C,ID/2,1,NN,NM1+1,WN(3*NN),IWORK)
        ENDIF
        CALL MFFTIV(C,ID/2,1,NN,NM1+1,WN,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTRI(C,1,ID/2,NM1,NN,WM(4*NM))
        CALL MFFTOV(C,1,ID/2,NM1,NN,WM(3*NM1),IWORK)
        CALL MFFTIV(C,1,ID/2,NM1,NN,WM,IERR)
        IF(IERR.NE.0)RETURN

      ELSE
        CALL MFFTDV(C,1,ID/2,NM1,NN,WM,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTOV(C,1,ID/2,NM1,NN,WM(NM1*2),IWORK)
        CALL MFFTRD(C,1,ID/2,NM1,NN,WM(4*NM))
        CALL MFFTDV(C,ID/2,1,NN,NM1+1,WN,IERR)
        IF(IERR.NE.0)RETURN

        IF(IORD.NE.0) THEN
          CALL MFFTOV(C,ID/2,1,NN,NM1+1,WN(NN*2),IWORK)
        ENDIF
C
      ENDIF
C
      END
      SUBROUTINE R3FFT(C,ID,NL,NM,NN,WL,WM,WN,IOPT,ISIG,IORD,IWORK,IERR)
C
C***PURPOSE:
C       THIS ROUTINE PERFORMS A 3-DIMENSIONAL REAL FOURIER TRANSFORM,
C       OF ORDER NL*NM*NN  .
C***USAGE:
C       THE USER IS EXPECTED TO PROVIDE THE DATA IN A 3-DIMENSIONAL
C       REAL ARRAY C, DIMENSIONED IN THE CALLING PROGRAM C(ID,NM,NN);
C       ID HAS TO BE AN EVEN INTEGER, EQUAL TO NL+2.
C       FOR OUTPUT DATA ARRENGEMENT SEE NOTES TO R2FFT HERE ABOVE.
C       THIS ROUTINE IS
C       INTENDED FOR REPEATED USAGE, THUS SEPARATE SET-UP AND
C       AND OPERATING CALLS ARE AVAILABLE: THE USER SHOULD IN ANY CASE
C       PERFORM A SET-UP CALL (ISIG=0) PASSING THE PARAMETERS BEFORE
C       PERFORMING AN ACTUAL TRANSFORM ( ISIG= +1 OR -1 ); THE USER CAN
C       CHOOSE WHETHER TO OBTAIN THE RESULTS OF THE DIRECT TRANSFORM
C       IN NATURAL ORDER (ISIG=-1,IORD=1) OR LEAVE THEM IN THE
C       BIT-REVERSED  ORDER( ISIG=-1,IORD=0); THIS CHOICE SAVES
C       SOME COMPUTER TIME, AND IT IS RECOMMENDED IN CASES DISCUSSED
C       IN THE LONG WRITE-UP. ANALOGOUSLY, THE INVERSE TRANSFORM ACCEPTS
C       INPUT ( PLEASE NOTE! ) DATA IN NATURAL ORDER ( ISIG=1,IORD=1),
C       OR DATA ALREADY SUBJECTED TO A BIT-REVERSAL PERMUTATION( ISIG=1
C       IORD=0).
C       A SPECIAL TREATMENT IS AVAILABLE TO SPEED UP THE TRANSFORM OF
C       SMALL MATRICES. THIS TREATMENT IS ACTIVATED BY THE FLAG IOPT. IN
C       THIS CASE THE TABLES FOR THE SECOND DIMENSION ( WM ) ARE LARGER,
C       BUT THE INCREASE IN PERFORMANCE IS SUBSTANTIAL WHEN NM<32.
C***ARGUMENTS :
C       INPUT :
C       C : ARRAY TO BE TRANSFORMED; DECLARED COMPLEX C(ID,NM,NN) IN THE
C           CALLING PROGRAM;
C       ID : FIRST DIMENSION OF C IN THE CALLING PROGRAM
C            IT HAS TO BE AN EVEN INTEGER .GE. NL+2.
C       ISIG : OPTION FLAG : ISIG=0 : SET-UP RUN, C NOT USED
C                            ISIG=-1: DIRECT TRANSFORM
C                            ISIG=+1: INVERSE TRANSFORM
C       WL,WM,WN : INTEGER ARRAYS,USED TO HOST TABLES FOR THE TRANSFORMS
C               DIMENSIONED IN THE CALLING PROGRAM AT LEAST (6*NL+14)
C               (4*NM+14) AND (4*NN+14) RESPECTIVELY; IF
C               IOPT=1, WM MUST BE DIMENSIONED AT LEAST 4*NM*(ID/2+1)+14
C               IF ISIG.NE.0, THEY ARE ASSUMED TO HAVE BEEN SET BY A
C               PREVIOUS CALL WITH ISIG=0 AND OTHER ARGUMENTS EQUAL, AND
C               NEVER HAVE BEEN MODIFIED ;
C               WHEN THE CORRESPONDING ORDERS ARE EQUAL, THEY DO NOT
C               NEED TO BE DISTINCT
C       NL : ORDER OF THE TRANSFORM ALONG THE COLUMNS OF C
C            IT HAS TO BE AN EVEN INTEGER.
C       NM : ORDER OF THE TRANSFORM ALONG THE ROWS OF C
C       NN : ORDER OF THE TRANSFORM ALONG THE THIRD DIMENSION OF C
C       IOPT : OPTION FLAG : =0 : NORMAL TREATMENT
C                            =1 : SPECIAL TREATMENT FOR IMPROVING
C                                 VECTORIZATION ON MATRICES WITH
C                                 SMALL NL; REQUIRES LONG WM(SEE);IF
C                                 REQUESTED, MUST BE PRESENT IN BOTH
C                                 THE SET-UP AND TRANSFORM CALLS;
C       IORD : OPTION FLAG : =1 : OUTPUT IN NATURAL ORDER (ISIG=-1)
C                                 INPUT IN NATURAL ORDER  (ISIG=+1)
C                            =0 : OUTPUT IN BIT-REVERSED ORDER(ISIG=-1)
C                                 INPUT IN BIT-REVERSED ORDER(ISIG=+1)
C       IWORK : INTEGER ARRAY, USED AS WORK AREA FOR REORDERING IF
C               IORD=1; IT MUST BE AT LEAST MAX(NL,NM,NN) WORDS LONG.
C
C        OUTPUT :
C       C : TRANSFORMED ARRAY
C       WL, WM, WN : ONLY IF ISIG=0, WL,WM AND WN ARE FILLED WITH THE
C                     APPROPRIATE TABLES
C       IWORK : UNDEFINED
C       IERR  : ERROR CODE : =0 : SUCCESSFUL
C                            =1 : WRONG ID PARAMETER
C                            =2 : PRIME FACTORS DIFFERENT FROM 2,3,5
C                                 ARE PRESENT IN DATA DIMENSIONS
C                            =3 : TABLES NOT CORRECTLY INITIALIZED
C                            =4 : FIRST DIMENSION IS AN ODD NUMBER
C
      COMPLEX C(*)
      INTEGER WL(-14:*),WM(-14:*),WN(-14:*)
      INTEGER  IWORK(*)
C
C
      INTEGER IDERR,FACERR,TBERR,ODDERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3,ODDERR=4)
C
      IF(ID.LT.NL+2)THEN
        IERR=IDERR
        RETURN
      ENDIF
      NL1=NL/2
      IF(NL1*2.NE.NL)THEN
        IERR=ODDERR
        RETURN
      ENDIF
      IERR=0

      NMPN=NM*NN
C
C
      IF(ISIG.EQ.0) THEN
        CALL MFFTP(NM,WM,ID/2*IOPT,IERR)
        IF(IERR.NE.0)RETURN

        CALL MFFTRP(NL,WL(4*NL))
        IF(NL1.NE.NM) THEN
          CALL MFFTP(NL1,WL,0,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTZ0(WM,1,2*NM+7,WL,1)
        ENDIF
C
        IF(NM.EQ.NN) THEN
          CALL MFFTZ0(WM,1,2*NM+7,WN,1)
        ELSE IF(NN.EQ.NL1) THEN
          CALL MFFTZ0(WL,1,2*NL1+7,WN,1)
        ELSE
          CALL MFFTP(NN,WN,0,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
        RETURN
C
      ELSE   IF(ISIG.GT.0) THEN
C
        IF(IORD.NE.0) THEN
          CALL MFFTOM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM(NM*3),IWORK)
          CALL MFFTOV(C,ID/2*NM,1,NN,ID/2*NM,WN(NN*3),IWORK)
        ENDIF
C
        CALL MFFTIV(C,ID/2*NM,1,NN,ID/2*NM,WN,IERR)
        IF(IERR.NE.0)RETURN
C
        IF(IOPT.EQ.0) THEN
          CALL MFFTIM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTIS(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
C
        CALL MFFTRI(C,1,ID/2,NL1,NMPN,WL(4*NL))
        CALL MFFTOV(C,1,ID/2,NL1,NMPN,WL(NL1*3),IWORK)
        CALL MFFTIV(C,1,ID/2,NL1,NMPN,WL,IERR)
        IF(IERR.NE.0)RETURN
C
C
      ELSE
C
C
        CALL MFFTDV(C,1,ID/2,NL1,NMPN,WL,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTOV(C,1,ID/2,NL1,NMPN,WL(NL1*2),IWORK)
        CALL MFFTRD(C,1,ID/2,NL1,NMPN,WL(4*NL))
C
C
        IF(IOPT.EQ.0) THEN
          CALL MFFTDM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTDS(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM,IERR)
          IF(IERR.NE.0)RETURN
        ENDIF
C
        CALL MFFTDV(C,ID/2*NM,1,NN,ID/2*NM,WN,IERR)
        IF(IERR.NE.0)RETURN
C
        IF(IORD.NE.0) THEN
          CALL MFFTOV(C,ID/2*NM,1,NN,ID/2*NM,WN(NN*2),IWORK)
          CALL MFFTOM(C,ID/2,ID/2*NM,1,NM,NN,NL1+1,WM(NM*2),IWORK)
        ENDIF
C
      ENDIF
C
      END
      SUBROUTINE MFFTA4(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 2 STEP APPLIED TO A VECTOR-OF
C       VECTORS-OF-COMPLEX C[IVS,NV [IES,NE]]. SEE REF.[1] FOR NOTATIONS
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDV, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:1),FAC(0:*)
      COMPLEX T0,F


        IF(MX.GT.2*LX)THEN
          DO 100 LAM=0,LLIM,LSTEP
            MUF=LX
            DO 90 MU=LAM+MSTEP,LAM+MD2LIM,MSTEP
              F=CONJG(FAC(MUF))
              DO 80 I=MU,MU+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=(T0-C(I,1))*F
80            CONTINUE
              MUF=MUF+LX
90          CONTINUE
            MUF=MUF-LX
            DO 91 MU=LAM+MD2LIM+2*MSTEP,LAM+MLIM,MSTEP
              F=-FAC(MUF)
              DO 81 I=MU,MU+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=(T0-C(I,1))*F
81            CONTINUE
              MUF=MUF-LX
91          CONTINUE
            DO 82 I=MD2LIM+MSTEP+LAM,MD2LIM+MSTEP+LAM+ILIM,IES
               T0=C(I,0)
              C(I,0)=T0+C(I,1)
              C(I,1)=CMPLX(AIMAG(T0-C(I,1)),-REAL(T0-C(I,1)))
82          CONTINUE
            DO 83 I=LAM,LAM+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=T0-C(I,1)
83          CONTINUE
100       CONTINUE
        ELSE
          IF(MX.EQ.1)GOTO 1000
C IF MX > 1 COME HERE
            MUF=LX
            DO 200 MU=MSTEP,MD2LIM,MSTEP
              F=CONJG(FAC(MUF))
              DO 190 LAM=MU,MU+LLIM,LSTEP
                DO 180 I=LAM,LAM+ILIM,IES
                  T0=C(I,0)
                  C(I,0)=T0+C(I,1)
                  C(I,1)=(T0-C(I,1))*F
  180           CONTINUE
  190         CONTINUE
              MUF=MUF+LX
200         CONTINUE
            MUF=MUF-LX
            DO 201 MU=MD2LIM+2*MSTEP,MLIM,MSTEP
              F=-FAC(MUF)
              DO 191 LAM=MU,MU+LLIM,LSTEP
                DO 181 I=LAM,LAM+ILIM,IES
                  T0=C(I,0)
                  C(I,0)=T0+C(I,1)
                  C(I,1)=(T0-C(I,1))*F
  181           CONTINUE
  191         CONTINUE
              MUF=MUF-LX
201         CONTINUE
             DO 192 LAM=MD2LIM+MSTEP,MD2LIM+MSTEP+LLIM,LSTEP
             DO 182 I=LAM,LAM+ILIM,IES
               T0=C(I,0)
               C(I,0)=T0+C(I,1)
               C(I,1)=CMPLX(AIMAG(T0-C(I,1)),-REAL(T0-C(I,1)))
182           CONTINUE
192         CONTINUE
1000        DO 193 LAM=0,LLIM,LSTEP
              DO 183 I=LAM,LAM+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=T0-C(I,1)
183           CONTINUE
193         CONTINUE
        ENDIF
      END
      SUBROUTINE MFFTA5(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 2 STEP APPLIED TO A VECTOR-OF
C       VECTORS-OF-COMPLEX C[IVS,NV [IES,NE]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIV, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:*),T0,F
      COMPLEX FAC(0:*)

        IF(2*MX.GE.LX)THEN
          IF(LX.EQ.1)GOTO 1000
          LAMF=MX
          DO 100 LAM=LSTEP,LD2LIM,LSTEP
            F=FAC(LAMF)
            DO 90 MU=LAM,LAM+MLIM,MSTEP
              DO 80 I=MU,MU+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
80            CONTINUE
90          CONTINUE
            LAMF=LAMF+MX
100       CONTINUE
          LAMF=LAMF-MX
          DO 101 LAM=LD2LIM+2*LSTEP,LLIM,LSTEP
            F=-CONJG(FAC(LAMF))
            DO 91 MU=LAM,LAM+MLIM,MSTEP
              DO 81 I=MU,MU+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
81            CONTINUE
91          CONTINUE
            LAMF=LAMF-MX
101       CONTINUE
          DO 93 MU=LD2LIM+LSTEP,LD2LIM+LSTEP+MLIM,MSTEP

            DO 83 I=MU,MU+ILIM,IES
              T0=CMPLX(-AIMAG(C(I,1)),REAL(C(I,1)))
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
83          CONTINUE
93        CONTINUE
1000      DO 92 MU=0,MLIM,MSTEP

            DO 82 I=MU,MU+ILIM,IES
              T0=C(I,1)
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
82          CONTINUE
92        CONTINUE
        ELSE
          DO 200 MU=0,MLIM,MSTEP
            LAMF=MX
            DO 190 LAM=MU+LSTEP,MU+LD2LIM,LSTEP
              F=FAC(LAMF)
              DO 180 I=LAM,LAM+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
180           CONTINUE
              LAMF=LAMF+MX
190         CONTINUE
            LAMF=LAMF-MX
            DO 191 LAM=MU+LD2LIM+2*LSTEP,MU+LLIM,LSTEP
              F=-CONJG(FAC(LAMF))
              DO 181 I=LAM,LAM+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
181           CONTINUE
              LAMF=LAMF-MX
191         CONTINUE
            DO 182 I=MU,MU+ILIM,IES
              T0=C(I,1)
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
182         CONTINUE
            DO 183 I=MU+LD2LIM+LSTEP,MU+LD2LIM+LSTEP+ILIM,IES
              T0=CMPLX(-AIMAG(C(I,1)),REAL(C(I,1)))
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
183         CONTINUE
200       CONTINUE
        ENDIF
      END
      SUBROUTINE MFFTA6(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 2 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:1),FAC(0:*)
      COMPLEX T0,F


        IF(MX.GT.2*LX)THEN
          DO 100 LAM=0,LLIM,LSTEP
            MUF=LX
            DO 90 MU=LAM+MSTEP,LAM+MD2LIM,MSTEP
              F=CONJG(FAC(MUF))
              DO 80  IV=MU,MU+IVLIM,IVS
              DO 80 I=IV,IV+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=(T0-C(I,1))*F
80            CONTINUE
              MUF=MUF+LX
90          CONTINUE
            MUF=MUF-LX
            DO 91 MU=LAM+MD2LIM+2*MSTEP,LAM+MLIM,MSTEP
              F=-FAC(MUF)
              DO 81 IV=MU,MU+IVLIM,IVS
              DO 81 I=IV,IV+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=(T0-C(I,1))*F
81            CONTINUE
              MUF=MUF-LX
91          CONTINUE
            DO 82 IV=MD2LIM+MSTEP+LAM,MD2LIM+MSTEP+LAM+IVLIM,IVS
            DO 82 I=IV,IV+ILIM,IES
               T0=C(I,0)
              C(I,0)=T0+C(I,1)
              C(I,1)=CMPLX(AIMAG(T0-C(I,1)),-REAL(T0-C(I,1)))
82          CONTINUE
            DO 83 IV=LAM,LAM+IVLIM,IVS
            DO 83 I=IV,IV+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=T0-C(I,1)
83          CONTINUE
100       CONTINUE
        ELSE
          IF(MX.EQ.1)GOTO 1000
C IF MX > 1 COME HERE
            MUF=LX
            DO 200 MU=MSTEP,MD2LIM,MSTEP
              F=CONJG(FAC(MUF))
              DO 190 LAM=MU,MU+LLIM,LSTEP
                DO 180 IV=LAM,LAM+IVLIM,IVS
                DO 180 I=IV,IV+ILIM,IES
                  T0=C(I,0)
                  C(I,0)=T0+C(I,1)
                  C(I,1)=(T0-C(I,1))*F
  180           CONTINUE
  190         CONTINUE
              MUF=MUF+LX
200         CONTINUE
            MUF=MUF-LX
            DO 201 MU=MD2LIM+2*MSTEP,MLIM,MSTEP
              F=-FAC(MUF)
              DO 191 LAM=MU,MU+LLIM,LSTEP
                DO 181 IV=LAM,LAM+IVLIM,IVS
                DO 181 I=IV,IV+ILIM,IES
                  T0=C(I,0)
                  C(I,0)=T0+C(I,1)
                  C(I,1)=(T0-C(I,1))*F
  181           CONTINUE
  191         CONTINUE
              MUF=MUF-LX
201         CONTINUE
             DO 192 LAM=MD2LIM+MSTEP,MD2LIM+MSTEP+LLIM,LSTEP
             DO 182 IV=LAM,LAM+IVLIM,IVS
             DO 182 I=IV,IV+ILIM,IES
               T0=C(I,0)
               C(I,0)=T0+C(I,1)
               C(I,1)=CMPLX(AIMAG(T0-C(I,1)),-REAL(T0-C(I,1)))
182           CONTINUE
192         CONTINUE
1000        DO 193 LAM=0,LLIM,LSTEP
              DO 183 IV=LAM,LAM+IVLIM,IVS
              DO 183 I=IV,IV+ILIM,IES
                T0=C(I,0)
                C(I,0)=T0+C(I,1)
                C(I,1)=T0-C(I,1)
183           CONTINUE
193         CONTINUE
        ENDIF
      END
      SUBROUTINE MFFTA7(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 2 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:1),FAC(0:*)
      COMPLEX T0,F

        IF(2*MX.GT.LX)THEN
          IF(LX.EQ.1)GOTO 1000
C   ELSE HERE
          LAMF=MX
          DO 100 LAM=LSTEP,LD2LIM,LSTEP
            F=FAC(LAMF)
            DO 90 MU=LAM,LAM+MLIM,MSTEP
              DO 80 IV=MU,MU+IVLIM,IVS
              DO 80 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
80            CONTINUE
90          CONTINUE
            LAMF=LAMF+MX
100       CONTINUE
          LAMF=LAMF-MX
          DO 101 LAM=LD2LIM+2*LSTEP,LLIM,LSTEP
            F=-CONJG(FAC(LAMF))
            DO 91 MU=LAM,LAM+MLIM,MSTEP
              DO 81 IV=MU,MU+IVLIM,IVS
              DO 81 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
81            CONTINUE
91          CONTINUE
            LAMF=LAMF-MX
101       CONTINUE
          DO 93 MU=LD2LIM+LSTEP,LD2LIM+LSTEP+MLIM,MSTEP

            DO 83 IV=MU,MU+IVLIM,IVS
            DO 83 I=IV,IV+ILIM,IES
              T0=CMPLX(-AIMAG(C(I,1)),REAL(C(I,1)))
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
83          CONTINUE
93        CONTINUE
1000      DO 92 MU=0,MLIM,MSTEP

            DO 82 IV=MU,MU+IVLIM,IVS
            DO 82 I=IV,IV+ILIM,IES
              T0=C(I,1)
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
82          CONTINUE
92        CONTINUE
        ELSE
          DO 200 MU=0,MLIM,MSTEP
            LAMF=MX
            DO 190 LAM=MU+LSTEP,MU+LD2LIM,LSTEP
              F=FAC(LAMF)
              DO 180 IV=LAM,LAM+IVLIM,IVS
              DO 180 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
180           CONTINUE
              LAMF=LAMF+MX
190         CONTINUE
            LAMF=LAMF-MX
            DO 191 LAM=MU+LD2LIM+2*LSTEP,MU+LLIM,LSTEP
              F=-CONJG(FAC(LAMF))
              DO 181 IV=LAM,LAM+IVLIM,IVS
              DO 181 I=IV,IV+ILIM,IES
                T0=C(I,1)*F
                C(I,1)=C(I,0)-T0
                C(I,0)=C(I,0)+T0
181           CONTINUE
              LAMF=LAMF-MX
191         CONTINUE
            DO 182 IV=MU,MU+IVLIM,IVS
            DO 182 I=IV,IV+ILIM,IES
              T0=C(I,1)
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
182         CONTINUE
            DO 183 IV=MU+LD2LIM+LSTEP,MU+LD2LIM+LSTEP+IVLIM,IVS
            DO 183 I=IV,IV+ILIM,IES
              T0=CMPLX(-AIMAG(C(I,1)),REAL(C(I,1)))
              C(I,1)=C(I,0)-T0
              C(I,0)=C(I,0)+T0
183         CONTINUE
200       CONTINUE
        ENDIF
      END
      SUBROUTINE MFFTA8(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX-2 STEP APPLIED TO A VECTOR-OF-
C       2-VECTORS-OF-COMPLEX [IMS,NM [IVS,NV [IES,NE]]], OPTIMIZED FOR
C       SMALL NE MATRICES.
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDM, WHICH CONTROLS
C       ITS OPERATION THROUGH COMMON MFFTPA.
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:1),FAC(0:*),T0
C
      IF (MX.NE.1) THEN
C
      DO 200 LAM=0,LLIM,LSTEP
        DO 150 IV=LAM,LAM+IVLIM,IVS
          IMUF=0
          DO 100 IMU=IV,IV+ILIM
            T0=C(IMU,0)
            C(IMU,0)=T0+C(IMU,1)
            C(IMU,1)=(T0-C(IMU,1))*FAC(IMUF)
            IMUF=IMUF+1
100       CONTINUE
150      CONTINUE
200   CONTINUE
C
      ELSE
        DO 400 LAM=0,LLIM,LSTEP
          DO 350 IV=LAM,LAM+IVLIM,IVS
            DO 300 IMU=IV,IV+ILIM
            T0=C(IMU,0)
            C(IMU,0)=T0+C(IMU,1)
            C(IMU,1)=T0-C(IMU,1)
 300        CONTINUE
 350      CONTINUE
 400    CONTINUE
      ENDIF
C
      END
      SUBROUTINE MFFTA9(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX-2 STEP APPLIED TO A VECTOR-OF-
C       2-VECTORS-OF-COMPLEX [IMS,NM [IVS,NV [IES,NE]]], OPTIMIZED FOR
C       SMALL NE MATRICES.
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIM, WHICH CONTROLS
C       ITS OPERATION THROUGH COMMON MFFTPA.
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:1),FAC(0:*),T0
C
      IF(LX.NE.1) THEN
C
      DO 200 MU=0,MLIM,MSTEP
        DO 150 IV=MU,MU+IVLIM,IVS
          ILAMF=0
          DO 100 ILAM=IV,IV+ILIM
            T0=C(ILAM,1)*FAC(ILAMF)
            C(ILAM,1)=C(ILAM,0)-T0
            C(ILAM,0)=C(ILAM,0)+T0
            ILAMF=ILAMF+1
100       CONTINUE
150     CONTINUE
200   CONTINUE
C
      ELSE
      DO 400 MU=0,MLIM,MSTEP
        DO 350 IV=MU,MU+IVLIM,IVS
          DO 300 ILAM=IV,IV+ILIM
          T0=C(ILAM,1)
          C(ILAM,1)=C(ILAM,0)-T0
          C(ILAM,0)=C(ILAM,0)+T0
 300      CONTINUE
 350    CONTINUE
 400  CONTINUE
      ENDIF
C
      END
      SUBROUTINE MFFTB4(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 3 STEP APPLIED TO A VECTOR-OF
C       VECTORS-OF-COMPLEX C[IVS,NV [IES,NE]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDV, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:2),FAC(0:*)
      COMPLEX T0,T1,T2,F1,F2

      REAL SIN60
      PARAMETER ( SIN60 =  8.6602540378443864E-1)

C..  MU > 1
            MUF=LX
            DO 200 MU=MSTEP,MLIM,MSTEP
              F1=CONJG(FAC(MUF))
              F2=CONJG(FAC(2*MUF))
              DO 190 LAM=MU,MU+LLIM,LSTEP
                DO 180 I=LAM,LAM+ILIM,IES
                T0=C(I,1)+C(I,2)
                T1=C(I,0)-0.5*T0
                T2=(C(I,1)-C(I,2))*SIN60
                C(I,0)=C(I,0)+T0
                C(I,1)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))*F1
                C(I,2)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))*F2
  180           CONTINUE
  190         CONTINUE
              MUF=MUF+LX
200         CONTINUE

C..  MU=0
1000        DO 193 LAM=0,LLIM,LSTEP
              DO 183 I=LAM,LAM+ILIM,IES
                T0=C(I,1)+C(I,2)
                T1=C(I,0)-0.5*T0
                T2=(C(I,1)-C(I,2))*SIN60
                C(I,0)=C(I,0)+T0
                C(I,1)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
                C(I,2)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
183           CONTINUE
193         CONTINUE
      END
      SUBROUTINE MFFTB5(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 3 STEP APPLIED TO A VECTOR-OF
C       VECTORS-OF-COMPLEX C[IVS,NV [IES,NE]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIV, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:2),FAC(0:*)
      COMPLEX T0,T1,T2,F1,F2
      REAL SIN60
      PARAMETER ( SIN60 =  8.6602540378443864E-1)

C..  LAM > 0
          LAMF=MX
          DO 100 LAM=LSTEP,LLIM,LSTEP
            F1=FAC(LAMF)
            F2=FAC(2*LAMF)
            DO 90 MU=LAM,LAM+MLIM,MSTEP
              DO 80 I=MU,MU+ILIM,IES
                T0=C(I,1)*F1+C(I,2)*F2
                T2=(C(I,1)*F1-C(I,2)*F2)*SIN60
                T1=C(I,0)-0.5*T0
                C(I,0)=C(I,0)+T0
                C(I,1)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
                C(I,2)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
80            CONTINUE
90          CONTINUE
            LAMF=LAMF+MX
100       CONTINUE

C.. LAM=0
          DO 92 MU=0,MLIM,MSTEP
            DO 82 I=MU,MU+ILIM,IES
                T0=C(I,1)+C(I,2)
                T2=(C(I,1)-C(I,2))*SIN60
                T1=C(I,0)-0.5*T0
                C(I,0)=C(I,0)+T0
                C(I,1)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
                C(I,2)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
82          CONTINUE
92        CONTINUE
      END
      SUBROUTINE MFFTB6(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 3 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:2),FAC(0:*)
      COMPLEX T0,T1,T2,F1,F2
      REAL SIN60
      PARAMETER ( SIN60 =  8.6602540378443864E-1)


C..  MU > 0
            MUF=LX
            DO 200 MU=MSTEP,MLIM,MSTEP
              F1=CONJG(FAC(MUF))
              F2=CONJG(FAC(MUF*2))
              DO 190 LAM=MU,MU+LLIM,LSTEP
                DO 180 IV=LAM,LAM+IVLIM,IVS
                DO 180 I=IV,IV+ILIM,IES
                  T0=C(I,1)+C(I,2)
                  T2=(C(I,1)-C(I,2))*SIN60
                  T1=C(I,0)-0.5*T0
                  C(I,0)=C(I,0)+T0
                  C(I,1)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))*F1
                  C(I,2)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))*F2
  180           CONTINUE
  190         CONTINUE
              MUF=MUF+LX
200         CONTINUE
            DO 193 LAM=0,LLIM,LSTEP
              DO 183 IV=LAM,LAM+IVLIM,IVS
              DO 183 I=IV,IV+ILIM,IES
                T0=C(I,1)+C(I,2)
                T2=(C(I,1)-C(I,2))*SIN60
                T1=C(I,0)-0.5*T0
                C(I,0)=C(I,0)+T0
                C(I,1)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
                C(I,2)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
183           CONTINUE
193         CONTINUE
      END
      SUBROUTINE MFFTB7(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 3 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:2),FAC(0:*)
      COMPLEX T0,T1,T2,F1,F2
      REAL SIN60
      PARAMETER ( SIN60 =  8.6602540378443864E-1)
C..  LAM>0
          LAMF=MX
          DO 100 LAM=LSTEP,LLIM,LSTEP
            F1=FAC(LAMF)
            F2=FAC(LAMF*2)
            DO 90 MU=LAM,LAM+MLIM,MSTEP
              DO 80 IV=MU,MU+IVLIM,IVS
              DO 80 I=IV,IV+ILIM,IES
                T0=C(I,1)*F1+C(I,2)*F2
                T2=(C(I,1)*F1-C(I,2)*F2)*SIN60
                T1=C(I,0)-0.5*T0
                C(I,0)=C(I,0)+T0
                C(I,1)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
                C(I,2)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
80            CONTINUE
90          CONTINUE
            LAMF=LAMF+MX
100       CONTINUE
          DO 92 MU=0,MLIM,MSTEP
            DO 82 IV=MU,MU+IVLIM,IVS
            DO 82 I=IV,IV+ILIM,IES
                T0=C(I,1)+C(I,2)
                T2=(C(I,1)-C(I,2))*SIN60
                T1=C(I,0)-0.5*T0
                C(I,0)=C(I,0)+T0
                C(I,1)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
                C(I,2)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
82          CONTINUE
92        CONTINUE
      END
      SUBROUTINE MFFTB8(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 3 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:1),FAC(0:*)
      COMPLEX T0,T1,T2,F1,F2
      REAL SIN60
      PARAMETER ( SIN60 =  8.6602540378443864E-1)


      IF(MX.NE.1)THEN
            DO 200 LAM=0,LLIM,LSTEP
              DO 190 IV=LAM,LAM+IVLIM,IVS
                IMUF=0
                DO 180 IMU=IV,IV+ILIM

                  T0=C(IMU,1)+C(IMU,2)
                  T2=(C(IMU,1)-C(IMU,2))*SIN60
                  T1=C(IMU,0)-0.5*T0
                  C(IMU,0)=C(IMU,0)+T0
                  C(IMU,1)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))*FAC(IMUF)
                  C(IMU,2)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))*FAC(IMUF+
     $                     NUSTEP)
                  IMUF=IMUF+1
  180           CONTINUE
  190         CONTINUE
200         CONTINUE
        ELSE
            DO 400 LAM=0,LLIM,LSTEP
              DO 390 IV=LAM,LAM+IVLIM,IVS
                DO 380 IMU=IV,IV+ILIM
                  T0=C(IMU,1)+C(IMU,2)
                  T2=(C(IMU,1)-C(IMU,2))*SIN60
                  T1=C(IMU,0)-0.5*T0
                  C(IMU,0)=C(IMU,0)+T0
                  C(IMU,1)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
                  C(IMU,2)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
380             CONTINUE
390           CONTINUE
400         CONTINUE
        ENDIF
      END
      SUBROUTINE MFFTB9(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 3 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:1),FAC(0:*)
      COMPLEX T0,T1,T2,F1,F2
      REAL SIN60
      PARAMETER ( SIN60 =  8.6602540378443864E-1)

      IF(LX.NE.1)THEN
          DO 200 MU=0,MLIM,MSTEP
            DO 150 IV=MU,MU+IVLIM,IVS
              ILAMF=0
              DO 100 ILAM=IV,IV+ILIM
                T0=C(ILAM,1)*FAC(ILAMF)+C(ILAM,2)*FAC(ILAMF+NUSTEP)
                T2=(C(ILAM,1)*FAC(ILAMF)-C(ILAM,2)*FAC(ILAMF+NUSTEP))*
     $             SIN60
                T1=C(ILAM,0)-0.5*T0
                C(ILAM,0)=C(ILAM,0)+T0
                C(ILAM,1)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
                C(ILAM,2)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
                ILAMF=ILAMF+1
100           CONTINUE
150         CONTINUE
200       CONTINUE
      ELSE
          DO 400 MU=0,MLIM,MSTEP
            DO 350 IV=MU,MU+IVLIM,IVS
              DO 300 ILAM=IV,IV+ILIM
                T0=C(ILAM,1)+C(ILAM,2)
                T2=(C(ILAM,1)-C(ILAM,2))*SIN60
                T1=C(ILAM,0)-0.5*T0
                C(ILAM,0)=C(ILAM,0)+T0
                C(ILAM,1)=(T1+CMPLX(-AIMAG(T2),REAL(T2)))
                C(ILAM,2)=(T1-CMPLX(-AIMAG(T2),REAL(T2)))
300           CONTINUE
350         CONTINUE
400       CONTINUE
      ENDIF
      END
      SUBROUTINE MFFTC4(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 5 STEP APPLIED TO A VECTOR-OF
C       VECTORS-OF-COMPLEX C[IVS,NVE[IES,NE]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDV, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:4),FAC(0:*)
      COMPLEX T1,T2,T3,T4,T5,F1,F2,F3,F4
      REAL SIN72,RAD5D4,S36D72
      PARAMETER (
     $ SIN72 =  9.51056516295153572116439333E-1,
     $ RAD5D4 =  5.59016994374947424102293417E-1,
     $ S36D72 =  6.18033988749894848204586834E-1 )
C..  MU > 1
            MUF=LX
            DO 200 MU=MSTEP,MLIM,MSTEP
              F1=CONJG(FAC(MUF))
              F2=CONJG(FAC(2*MUF))
              F3=CONJG(FAC(3*MUF))
              F4=CONJG(FAC(4*MUF))
              DO 190 LAM=MU,MU+LLIM,LSTEP
                DO 180 I=LAM,LAM+ILIM,IES
                  T1=C(I,1)+C(I,4)
                  T2=C(I,2)+C(I,3)
                  T3=(C(I,1)-C(I,4))*SIN72
                  T4=(C(I,2)-C(I,3))*SIN72
                  T5=T1+T2
                  T1=RAD5D4*(T1-T2)
                  T2=C(I,0)-0.25*T5
                  C(I,0)=C(I,0)+T5
                  T5=T2+T1
                  T2=T2-T1
                  T1=T3+S36D72*T4
                  T3=S36D72*T3-T4
                  C(I,1)=(T5-CMPLX(-AIMAG(T1),REAL(T1)))*F1
                  C(I,4)=(T5+CMPLX(-AIMAG(T1),REAL(T1)))*F4
                  C(I,2)=(T2-CMPLX(-AIMAG(T3),REAL(T3)))*F2
                  C(I,3)=(T2+CMPLX(-AIMAG(T3),REAL(T3)))*F3
  180           CONTINUE
  190         CONTINUE
              MUF=MUF+LX
200         CONTINUE
C..  MU=0
1000        DO 193 LAM=0,LLIM,LSTEP
              DO 183 I=LAM,LAM+ILIM,IES
                T1=C(I,1)+C(I,4)
                T2=C(I,2)+C(I,3)
                T3=(C(I,1)-C(I,4))*SIN72
                T4=(C(I,2)-C(I,3))*SIN72
                T5=T1+T2
                T1=RAD5D4*(T1-T2)
                T2=C(I,0)-0.25*T5
                C(I,0)=C(I,0)+T5
                T5=T2+T1
                T2=T2-T1
                T1=T3+S36D72*T4
                T3=S36D72*T3-T4
                C(I,1)=T5-CMPLX(-AIMAG(T1),REAL(T1))
                C(I,4)=T5+CMPLX(-AIMAG(T1),REAL(T1))
                C(I,2)=T2-CMPLX(-AIMAG(T3),REAL(T3))
                C(I,3)=T2+CMPLX(-AIMAG(T3),REAL(T3))
183           CONTINUE
193         CONTINUE
      END
      SUBROUTINE MFFTC5(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 5 STEP APPLIED TO A VECTOR-OF
C       VECTORS-OF-COMPLEX C[IVS,NV [IES,NE]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIV, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:4),FAC(0:*)
      COMPLEX T1,T2,T3,T4,T5,F1,F2,F3,F4
      REAL SINT2,RAD5D4,S36D72
      PARAMETER (
     $ SIN72 =  9.51056516295153572116439333E-1,
     $ RAD5D4 =  5.59016994374947424102293417E-1,
     $ S36D72 =  6.18033988749894848204586834E-1 )
C..  LAM > 0
          LAMF=MX
          DO 100 LAM=LSTEP,LLIM,LSTEP
            F1=FAC(LAMF)
            F2=FAC(2*LAMF)
            F3=FAC(3*LAMF)
            F4=FAC(4*LAMF)
            DO 90 MU=LAM,LAM+MLIM,MSTEP
              DO 80 I=MU,MU+ILIM,IES
                T1=C(I,1)*F1+C(I,4)*F4
                T2=C(I,2)*F2+C(I,3)*F3
                T3=(C(I,1)*F1-C(I,4)*F4)*SIN72
                T4=(C(I,2)*F2-C(I,3)*F3)*SIN72
                T5=T1+T2
                T1=RAD5D4*(T1-T2)
                T2=C(I,0)-0.25*T5
                C(I,0)=C(I,0)+T5
                T5=T2+T1
                T2=T2-T1
                T1=T3+S36D72*T4
                T3=S36D72*T3-T4
                C(I,1)=T5+CMPLX(-AIMAG(T1),REAL(T1))
                C(I,4)=T5-CMPLX(-AIMAG(T1),REAL(T1))
                C(I,2)=T2+CMPLX(-AIMAG(T3),REAL(T3))
                C(I,3)=T2-CMPLX(-AIMAG(T3),REAL(T3))
80            CONTINUE
90          CONTINUE
            LAMF=LAMF+MX
100       CONTINUE

C.. LAM=0
          DO 92 MU=0,MLIM,MSTEP
            DO 82 I=MU,MU+ILIM,IES
              T1=C(I,1)+C(I,4)
              T2=C(I,2)+C(I,3)
              T3=(C(I,1)-C(I,4))*SIN72
              T4=(C(I,2)-C(I,3))*SIN72
              T5=T1+T2
              T1=RAD5D4*(T1-T2)
              T2=C(I,0)-0.25*T5
              C(I,0)=C(I,0)+T5
              T5=T2+T1
              T2=T2-T1
              T1=T3+S36D72*T4
              T3=S36D72*T3-T4
              C(I,1)=T5+CMPLX(-AIMAG(T1),REAL(T1))
              C(I,4)=T5-CMPLX(-AIMAG(T1),REAL(T1))
              C(I,2)=T2+CMPLX(-AIMAG(T3),REAL(T3))
              C(I,3)=T2-CMPLX(-AIMAG(T3),REAL(T3))
82          CONTINUE
92        CONTINUE
      END
      SUBROUTINE MFFTC6(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 5 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:4),FAC(0:*)
      COMPLEX T1,T2,T3,T4,T5,F1,F2,F3,F4
      REAL SIN72,RAD5D4,S36D72
      PARAMETER (
     $ SIN72 =  9.51056516295153572116439333E-1,
     $ RAD5D4 =  5.59016994374947424102293417E-1,
     $ S36D72 =  6.18033988749894848204586834E-1 )
C..  MU > 0
            MUF=LX
            DO 200 MU=MSTEP,MLIM,MSTEP
              F1=CONJG(FAC(MUF))
              F2=CONJG(FAC(MUF*2))
              F3=CONJG(FAC(MUF*3))
              F4=CONJG(FAC(MUF*4))
              DO 190 LAM=MU,MU+LLIM,LSTEP
                DO 180 IV=LAM,LAM+IVLIM,IVS
                DO 180 I=IV,IV+ILIM,IES
                  T1=C(I,1)+C(I,4)
                  T2=C(I,2)+C(I,3)
                  T3=(C(I,1)-C(I,4))*SIN72
                  T4=(C(I,2)-C(I,3))*SIN72
                  T5=T1+T2
                  T1=RAD5D4*(T1-T2)
                  T2=C(I,0)-0.25*T5
                  C(I,0)=C(I,0)+T5
                  T5=T2+T1
                  T2=T2-T1
                  T1=T3+S36D72*T4
                  T3=S36D72*T3-T4
                  C(I,1)=(T5-CMPLX(-AIMAG(T1),REAL(T1)))*F1
                  C(I,4)=(T5+CMPLX(-AIMAG(T1),REAL(T1)))*F4
                  C(I,2)=(T2-CMPLX(-AIMAG(T3),REAL(T3)))*F2
                  C(I,3)=(T2+CMPLX(-AIMAG(T3),REAL(T3)))*F3
  180           CONTINUE
  190         CONTINUE
              MUF=MUF+LX
200         CONTINUE
            DO 193 LAM=0,LLIM,LSTEP
              DO 183 IV=LAM,LAM+IVLIM,IVS
              DO 183 I=IV,IV+ILIM,IES
                T1=C(I,1)+C(I,4)
                T2=C(I,2)+C(I,3)
                T3=(C(I,1)-C(I,4))*SIN72
                T4=(C(I,2)-C(I,3))*SIN72
                T5=T1+T2
                T1=RAD5D4*(T1-T2)
                T2=C(I,0)-0.25*T5
                C(I,0)=C(I,0)+T5
                T5=T2+T1
                T2=T2-T1
                T1=T3+S36D72*T4
                T3=S36D72*T3-T4
                C(I,1)=T5-CMPLX(-AIMAG(T1),REAL(T1))
                C(I,4)=T5+CMPLX(-AIMAG(T1),REAL(T1))
                C(I,2)=T2-CMPLX(-AIMAG(T3),REAL(T3))
                C(I,3)=T2+CMPLX(-AIMAG(T3),REAL(T3))
183           CONTINUE
193         CONTINUE
      END
      SUBROUTINE MFFTC7(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 5 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIM, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:4),FAC(0:*)
      COMPLEX T1,T2,T3,T4,T5,F1,F2,F3,F4
      REAL SIN72,RAD5D4,S36D72
      PARAMETER (
     $ SIN72 =  9.51056516295153572116439333E-1,
     $ RAD5D4 =  5.59016994374947424102293417E-1,
     $ S36D72 =  6.18033988749894848204586834E-1 )
C..  LAM>0
          LAMF=MX
          DO 100 LAM=LSTEP,LLIM,LSTEP
            F1=FAC(LAMF)
            F2=FAC(LAMF*2)
            F3=FAC(LAMF*3)
            F4=FAC(LAMF*4)
            DO 90 MU=LAM,LAM+MLIM,MSTEP
              DO 80 IV=MU,MU+IVLIM,IVS
              DO 80 I=IV,IV+ILIM,IES
                T1=C(I,1)*F1+C(I,4)*F4
                T2=C(I,2)*F2+C(I,3)*F3
                T3=(C(I,1)*F1-C(I,4)*F4)*SIN72
                T4=(C(I,2)*F2-C(I,3)*F3)*SIN72
                T5=T1+T2
                T1=RAD5D4*(T1-T2)
                T2=C(I,0)-0.25*T5
                C(I,0)=C(I,0)+T5
                T5=T2+T1
                T2=T2-T1
                T1=T3+S36D72*T4
                T3=S36D72*T3-T4
                C(I,1)=T5+CMPLX(-AIMAG(T1),REAL(T1))
                C(I,4)=T5-CMPLX(-AIMAG(T1),REAL(T1))
                C(I,2)=T2+CMPLX(-AIMAG(T3),REAL(T3))
                C(I,3)=T2-CMPLX(-AIMAG(T3),REAL(T3))
80            CONTINUE
90          CONTINUE
            LAMF=LAMF+MX
100       CONTINUE
          DO 92 MU=0,MLIM,MSTEP
            DO 82 IV=MU,MU+IVLIM,IVS
            DO 82 I=IV,IV+ILIM,IES
              T1=C(I,1)+C(I,4)
              T2=C(I,2)+C(I,3)
              T3=(C(I,1)-C(I,4))*SIN72
              T4=(C(I,2)-C(I,3))*SIN72
              T5=T1+T2
              T1=RAD5D4*(T1-T2)
              T2=C(I,0)-0.25*T5
              C(I,0)=C(I,0)+T5
              T5=T2+T1
              T2=T2-T1
              T1=T3+S36D72*T4
              T3=S36D72*T3-T4
              C(I,1)=T5+CMPLX(-AIMAG(T1),REAL(T1))
              C(I,4)=T5-CMPLX(-AIMAG(T1),REAL(T1))
              C(I,2)=T2+CMPLX(-AIMAG(T3),REAL(T3))
              C(I,3)=T2-CMPLX(-AIMAG(T3),REAL(T3))
82          CONTINUE
92        CONTINUE
      END
      SUBROUTINE MFFTC8(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY GENTLEMAN-SANDE RADIX 5 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTDS, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:4),FAC(0:*)
      COMPLEX T1,T2,T3,T4,T5
      REAL SIN72,RAD5D4,S32D72
      PARAMETER (
     $ SIN72 =  9.51056516295153572116439333E-1,
     $ RAD5D4 =  5.59016994374947424102293417E-1,
     $ S36D72 =  6.18033988749894848204586834E-1 )
      IF(MX.NE.1)THEN
            DO 200 LAM=0,LLIM,LSTEP
              DO 190 IV=LAM,LAM+IVLIM,IVS
                IMUF=0
                IMUF2=NUSTEP
                IMUF3=2*NUSTEP
                IMUF4=3*NUSTEP
                DO 180 IMU=IV,IV+ILIM
                  T1=C(IMU,1)+C(IMU,4)
                  T2=C(IMU,2)+C(IMU,3)
                  T3=(C(IMU,1)-C(IMU,4))*SIN72
                  T4=(C(IMU,2)-C(IMU,3))*SIN72
                  T5=T1+T2
                  T1=RAD5D4*(T1-T2)
                  T2=C(IMU,0)-0.25*T5
                  C(IMU,0)=C(IMU,0)+T5
                  T5=T2+T1
                  T2=T2-T1
                  T1=T3+S36D72*T4
                  T3=S36D72*T3-T4
                  C(IMU,1)=(T5-CMPLX(-AIMAG(T1),REAL(T1)))*FAC(IMUF )
                  C(IMU,4)=(T5+CMPLX(-AIMAG(T1),REAL(T1)))*FAC(IMUF4)
                  C(IMU,2)=(T2-CMPLX(-AIMAG(T3),REAL(T3)))*FAC(IMUF2)
                  C(IMU,3)=(T2+CMPLX(-AIMAG(T3),REAL(T3)))*FAC(IMUF3)
                  IMUF=IMUF+1
                  IMUF2=IMUF2+1
                  IMUF3=IMUF3+1
                  IMUF4=IMUF4+1
  180           CONTINUE
  190         CONTINUE
200         CONTINUE
        ELSE
            DO 400 LAM=0,LLIM,LSTEP
              DO 390 IV=LAM,LAM+IVLIM,IVS
                DO 380 IMU=IV,IV+ILIM
                  T1=C(IMU,1)+C(IMU,4)
                  T2=C(IMU,2)+C(IMU,3)
                  T3=(C(IMU,1)-C(IMU,4))*SIN72
                  T4=(C(IMU,2)-C(IMU,3))*SIN72
                  T5=T1+T2
                  T1=RAD5D4*(T1-T2)
                  T2=C(IMU,0)-0.25*T5
                  C(IMU,0)=C(IMU,0)+T5
                  T5=T2+T1
                  T2=T2-T1
                  T1=T3+S36D72*T4
                  T3=S36D72*T3-T4
                  C(IMU,1)=T5-CMPLX(-AIMAG(T1),REAL(T1))
                  C(IMU,4)=T5+CMPLX(-AIMAG(T1),REAL(T1))
                  C(IMU,2)=T2-CMPLX(-AIMAG(T3),REAL(T3))
                  C(IMU,3)=T2+CMPLX(-AIMAG(T3),REAL(T3))
380             CONTINUE
390           CONTINUE
400         CONTINUE
        ENDIF
      END
      SUBROUTINE MFFTC9(C,FAC)
C
C   PURPOSE:
C       ELEMENTARY COOLEY-TUKEY RADIX 5 STEP APPLIED TO A VECTOR-OF
C       2-VECTORS-OF-COMPLEX C[IMS,NM [IVS,NV [IES,NE]]].
C       SEE REF.[1] FOR NOTATIONS.
C       THIS ROUTINE CAN BE USED ONLY BY ROUTINE MFFTIS, WHICH CONTROLS
C       ITS OPERATION THROUGH THE MFFTPA COMMON
C
C   DUMMY ARGUMENTS :
C
C   C   ARRAY BEING FOURIER  TRANSFORMED
C   FAC PHASE FACTORS, PREPARED BY MFFTP; NOT MODIFIED IN OUTPUT
C
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      COMPLEX C(0:NUSTEP-1,0:4),FAC(0:*)
      COMPLEX T1,T2,T3,T4,T5
      REAL SIN72,RAD5D4,S36D72
      PARAMETER (
     $ SIN72 =  9.51056516295153572116439333E-1,
     $ RAD5D4 =  5.59016994374947424102293417E-1,
     $ S36D72 =  6.18033988749894848204586834E-1 )

      IF(LX.NE.1)THEN
          DO 200 MU=0,MLIM,MSTEP
            DO 150 IV=MU,MU+IVLIM,IVS
              ILAMF=0
              ILAMF2=NUSTEP
              ILAMF3=2*NUSTEP
              ILAMF4=3*NUSTEP
              DO 100 ILAM=IV,IV+ILIM
                T1=C(ILAM,1)*FAC(ILAMF)+C(ILAM,4)*FAC(ILAMF4)
                T2=C(ILAM,2)*FAC(ILAMF2)+C(ILAM,3)*FAC(ILAMF3)
                T3=(C(ILAM,1)*FAC(ILAMF)-C(ILAM,4)*FAC(ILAMF4))*SIN72
                T4=(C(ILAM,2)*FAC(ILAMF2)-C(ILAM,3)*FAC(ILAMF3))*SIN72
                T5=T1+T2
                T1=RAD5D4*(T1-T2)
                T2=C(ILAM,0)-0.25*T5
                C(ILAM,0)=C(ILAM,0)+T5
                T5=T2+T1
                T2=T2-T1
                T1=T3+S36D72*T4
                T3=S36D72*T3-T4
                C(ILAM,1)=T5+CMPLX(-AIMAG(T1),REAL(T1))
                C(ILAM,4)=T5-CMPLX(-AIMAG(T1),REAL(T1))
                C(ILAM,2)=T2+CMPLX(-AIMAG(T3),REAL(T3))
                C(ILAM,3)=T2-CMPLX(-AIMAG(T3),REAL(T3))
                ILAMF=ILAMF+1
                ILAMF2=ILAMF2+1
                ILAMF3=ILAMF3+1
                ILAMF4=ILAMF4+1
100           CONTINUE
150         CONTINUE
200       CONTINUE
      ELSE
          DO 400 MU=0,MLIM,MSTEP
            DO 350 IV=MU,MU+IVLIM,IVS
              DO 300 ILAM=IV,IV+ILIM
                T1=C(ILAM,1)+C(ILAM,4)
                T2=C(ILAM,2)+C(ILAM,3)
                T3=(C(ILAM,1)-C(ILAM,4))*SIN72
                T4=(C(ILAM,2)-C(ILAM,3))*SIN72
                T5=T1+T2
                T1=RAD5D4*(T1-T2)
                T2=C(ILAM,0)-0.25*T5
                C(ILAM,0)=C(ILAM,0)+T5
                T5=T2+T1
                T2=T2-T1
                T1=T3+S36D72*T4
                T3=S36D72*T3-T4
                C(ILAM,1)=T5+CMPLX(-AIMAG(T1),REAL(T1))
                C(ILAM,4)=T5-CMPLX(-AIMAG(T1),REAL(T1))
                C(ILAM,2)=T2+CMPLX(-AIMAG(T3),REAL(T3))
                C(ILAM,3)=T2-CMPLX(-AIMAG(T3),REAL(T3))
300           CONTINUE
350         CONTINUE
400       CONTINUE
      ENDIF
      END
      SUBROUTINE MFFTDM(C,IMSX,IVSX,IESX,NMX,NVX,NEX,TABLES,IERR)
C   PURPOSE:
C       THIS SUBROUTINE PERFORMS A DIRECT FOURIER TRANSFORM ALONG
C       THE SECOND DIMENSION OF A 3-DIMENSIONAL MATRIX, USING THE
C       GENTLEMAN-SANDE ALGORITHM.
C       THE SEQUENCE TO BE TRANSFORMED IS C[IMSX,NMX], WHOSE COMPONENTS
C       ARE THE 2-VECTORS C(M)[IVSX,NVX [IESX,NEX]].
C       SEE REF.[1] FOR NOTATIONS.
C  EXAMPLE:
C       LET C BE A 3-D MATRIX C(N1,N2,N3), DECLARED VIA
C              DIMENSION C(NN1,N2,N3)
C       WITH NN1.GE.N1
C       THEN ITS DFT ALONG THE SECOND DIMENSION IS OBTAINED BY
C              CALL MFFTDM(C,NN1,NN1*N2,1,N2,N3,N1,TABLES,IERR)
C  IMPLEMENTATION:
C       THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
C       "BUTTERFLY" TRANSFORMATION MFFT?6; PARAMETERS OF THE "BUTTERFLY"
C       ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
C  ARGUMENTS:
C    INPUT  :
C       C : ARRAY TO BE TRANSFORMED.
C       IMSX,IVSX,IESX,NMX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE
C           OF C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED
C           ON OUTPUT
C       TABLES : ARRAY PREPARED BY MFFTP. IT  IS NOT CHANGED ON OUTPUT.
C                IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
C                IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
C    OUTPUT:
C       C : TRANSFORM OF THE ORIGINAL ARRAY; "BIT REVERSED" ORDER
C       IERR : ERROR CODE : =0 : SUCCESSFUL
C                         : =3 :  'TABLES' NOT CORRECTLY INITIALIZED
      COMPLEX C(*)
      INTEGER TABLES(-14:*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
C  LOADING THE COMMON BLOCK : CONSTANTS
      IMS=IMSX
      IVS=IVSX
      IES=IESX
      NM=NMX
      NV=NVX
      NE=NEX

      IVLIM=(NV-1)*IVS
      ILIM=(NE-1)*IES
      MSTEP=IMS
C  LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
      LX=1
      MX=NM
C  SELECT THE HIGHEST FACTOR OF NM
      IFAC=TABLES(-1)
      GOTO(200,300,500)IFAC
      IERR=TBERR
      RETURN
C...  RADIX 5 LOOP
500   CONTINUE
      DO 510 IM=1,TABLES(-12)
        MX=MX/5
        NUSTEP=MX*MSTEP
        MLIM=NUSTEP-MSTEP
        LSTEP=NUSTEP*5
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTC6(C,TABLES(0))
        LX=LX*5
510   CONTINUE
C..  RADIX 3 LOOP
300   CONTINUE
      DO 310 IM=1,TABLES(-13)
        MX=MX/3
        NUSTEP=MX*MSTEP
        MLIM=NUSTEP-MSTEP
        LSTEP=NUSTEP*3
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTB6(C,TABLES(0))
        LX=LX*3
310   CONTINUE
C..  RADIX 2 LOOP
200   CONTINUE
      DO 210 IM=1,TABLES(-14)
        MX=MX/2
        NUSTEP=MX*MSTEP
        MLIM=NUSTEP-MSTEP
        MD2LIM=NUSTEP/2-MSTEP
        LSTEP=NUSTEP*2
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTA6(C,TABLES(0))
        LX=LX+LX
210   CONTINUE
      END
      SUBROUTINE MFFTDS(C,IMSX,IVSX,IESX,NMX,NVX,NEX,TABLES,IERR)

C   PURPOSE:
C       THE SAME AS MFFTDM.    IT IS A VARIANT OF MFFTDM
C       AIMED AT OPTIMAL PERFORMANCE ON MATRICES HAVING THE
C       FIRST DIMENSION SMALLER THAN 64. IT REQUIRES THAT MFFTP HAS
C       BEEN CALLED WITH ID .NE. 0
C  IMPLEMENTATION:
C       THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
C       "BUTTERFLY" TRANSFORMATION MFFT?8; PARAMETERS OF THE "BUTTERFLY"
C       ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
C  ARGUMENTS:
C    INPUT  :
C       C : ARRAY TO BE TRANSFORMED.
C       IMSX,IVSX,IESX,NMX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE
C           OF C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED
C           ON OUTPUT
C       TABLES : ARRAY PREPARED BY MFFTP. IT  IS NOT CHANGED ON OUTPUT.
C                IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
C                IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
C    OUTPUT:
C       C : TRANSFORM OF THE ORIGINAL ARRAY; "BIT REVERSED" ORDER
C       IERR : ERROR CODE : =0 : SUCCESSFUL
C                         : =3 :  'TABLES' NOT CORRECTLY INITIALIZED
      COMPLEX C(*)
      INTEGER TABLES(-14:*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      SAVE IBASE
C  LOADING THE COMMON BLOCK : CONSTANTS
      IMS=IMSX
      IVS=IVSX
      NM=NMX
      NV=NVX
      NE=NEX
      IVLIM=(NV-1)*IVS
      MSTEP=IMS
C  LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
      LX=1
      MX=NM

      IBASE=2*NM
C  SELECT THE HIGHEST FACTOR OF NM
      IFAC=TABLES(-1)
      GOTO(200,300,500)IFAC
      IERR=TBERR
      RETURN
C...  RADIX 5 LOOP
500   CONTINUE
      DO 510 IM=1,TABLES(-12)
        MX=MX/5
        NUSTEP=MX*MSTEP
        ILIM=NUSTEP-1
        LSTEP=NUSTEP*5
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTC8(C,TABLES(2*IBASE))
        LX=LX*5
        IBASE=IBASE+NUSTEP*4
510   CONTINUE
C..  RADIX 3 LOOP
300   CONTINUE
      DO 310 IM=1,TABLES(-13)
        MX=MX/3
        NUSTEP=MX*MSTEP
        ILIM=NUSTEP-1
        LSTEP=NUSTEP*3
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTB8(C,TABLES(2*IBASE))
        LX=LX*3
        IBASE=IBASE+NUSTEP*2
310   CONTINUE
C..  RADIX 2 LOOP
200   CONTINUE
      DO 210 IM=1,TABLES(-14)
        MX=MX/2
        NUSTEP=MX*MSTEP
        ILIM=NUSTEP-1
        MD2LIM=NUSTEP/2-MSTEP
        LSTEP=NUSTEP*2
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTA8(C,TABLES(2*IBASE))
        LX=LX+LX
        IBASE=IBASE+NUSTEP
210   CONTINUE

      END
      SUBROUTINE MFFTDV(C,IVSX,IESX,NVX,NEX,TABLES,IERR)
C   PURPOSE:
C       THIS SUBROUTINE PERFORMS A DIRECT FOURIER TRANSFORM ALONG
C       ONE DIMENSION OF A 2-DIMENSIONAL MATRIX, USING THE
C       GENTLEMAN-SANDE ALGORITHM; NO REORDERING IS PERFORMED.
C       THE SEQUENCE TO BE TRANSFORMED IS C[IVSX,NVX], WHOSE COMPONENTS
C       ARE THE VECTORS C(M)[IESX,NEX].
C       SEE REF.[1] FOR NOTATIONS.
C  EXAMPLE:
C       LET C BE A 2-D MATRIX C(N1,N2) DECLARED VIA
C                 DIMENSION C(ID,N2)
C       WITH ID.GE.N1.
C       THEN THE DFT ALONG THE FIRST DIMENSION IS OBTAINED BY
C                 CALL MFFTDV(C,1,ID,N1,N2,TABLES,IERR)
C       THE DFT ALONG THE SECOND DIMENSION IS OBTAINED BY
C                 CALL MFFTDV(C,ID,1,N2,N1,TABLES,IERR)
C  IMPLEMENTATION:
C       THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
C       "BUTTERFLY" TRANSFORMATION MFFT?4; PARAMETERS OF THE "BUTTERFLY"
C       ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
C  ARGUMENTS:
C       C : ARRAY TO BE TRANSFORMED.
C       IVSX,IESX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE OF
C           C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED ON
C           OUTPUT;
C       TABLES : ARRAY PREPARED BY MFFTP. IT  IS NOT CHANGED ON OUTPUT.
C                IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
C                IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
C    OUTPUT:
C       C : TRANSFORM OF THE ORIGINAL ARRAY; "BIT REVERSED" ORDER
C       IERR : ERROR CODE : =0 : SUCCESSFUL
C                         : =3 :  'TABLES' NOT CORRECTLY INITIALIZED
      COMPLEX C(*)
      INTEGER TABLES(-14:*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
C  LOADING THE COMMON BLOCK : CONSTANTS
      IVS=IVSX
      IES=IESX
      NV=NVX
      NE=NEX
      ILIM=(NE-1)*IES
      MSTEP=IVS
C  LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
      LX=1
      MX=NV
C  SELECT THE HIGHEST FACTOR OF NV
      IFAC=TABLES(-1)
      GOTO(200,300,500)IFAC
      IERR=TBERR
      RETURN
C...  RADIX 5 LOOP
500   CONTINUE
      DO 510 IM=1,TABLES(-12)
        MX=MX/5
        NUSTEP=MX*MSTEP
        MLIM=NUSTEP-MSTEP
        LSTEP=NUSTEP*5
        LLIM=NV*MSTEP-LSTEP
        CALL MFFTC4(C,TABLES(0))
        LX=LX*5
510   CONTINUE
C..  RADIX 3 LOOP
300   CONTINUE
      DO 310 IM=1,TABLES(-13)
        MX=MX/3
        NUSTEP=MX*MSTEP
        MLIM=NUSTEP-MSTEP
        LSTEP=NUSTEP*3
        LLIM=NV*MSTEP-LSTEP
        CALL MFFTB4(C,TABLES(0))
        LX=LX*3
310   CONTINUE
C..  RADIX 2 LOOP
200   CONTINUE
      DO 210 IM=1,TABLES(-14)
        MX=MX/2
        NUSTEP=MX*MSTEP
        MLIM=NUSTEP-MSTEP
        MD2LIM=NUSTEP/2-MSTEP
        LSTEP=NUSTEP*2
        LLIM=NV*MSTEP-LSTEP
        CALL MFFTA4(C,TABLES(0))
        LX=LX+LX
210   CONTINUE
      END
      SUBROUTINE MFFTIM(C,IMSX,IVSX,IESX,NMX,NVX,NEX,TABLES,IERR)

C   PURPOSE:
C       THIS SUBROUTINE PERFORMS AN INVERSE FOURIER TRANSFORM ALONG
C       THE SECOND DIMENSION OF A 3-DIMENSIONAL MATRIX, USING THE
C       COOLEY-TUKEY ALGORITHM.
C       THE INPUT SEQUENCE IS ASSUMED TO HAVE BEEN SUBJECTED TO A
C       "BIT REVERSAL" PERMUTATION, THROUGH A CALL TO MFFTOM, OR
C       BECAUSE IT IS THE OUTPUT OF MFFTDM.
C       THE SEQUENCE TO BE TRANSFORMED IS C[IMSX,NMX], WHOSE COMPONENTS
C       ARE THE 2-VECTORS C(M)[IVSX,NVX [IESX,NEX]].
C       SEE REF.[1] FOR NOTATIONS.
C  EXAMPLE:
C       LET C BE A 3-D MATRIX C(N1,N2,N3), DECLARED VIA
C              DIMENSION C(NN1,N2,N3)
C       WITH NN1.GE.N1
C       THEN ITS IDFT ALONG THE SECOND DIMENSION IS OBTAINED BY
C              CALL MFFTIM(C,NN1,NN1*N2,1,N2,N3,N1,TABLES,IERR)
C  IMPLEMENTATION:
C       THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
C       "BUTTERFLY" TRANSFORMATION MFFT?7; PARAMETERS OF THE "BUTTERFLY"
C       ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
C  ARGUMENTS:
C    INPUT  :
C       C : ARRAY TO BE TRANSFORMED.
C       IMSX,IVSX,IESX,NMX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE
C           OF C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED
C           ON OUTPUT
C       TABLES : ARRAY PREPARED BY MFFTP. IT  IS NOT CHANGED ON OUTPUT.
C                IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
C                IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
C    OUTPUT:
C       C : TRANSFORM OF THE ORIGINAL ARRAY;
C       IERR : ERROR CODE : =0 : SUCCESSFUL
C                         : =3 :  'TABLES' NOT CORRECTLY INITIALIZED
      COMPLEX C(*)
      INTEGER TABLES(-14:*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
C  LOADING THE COMMON BLOCK : CONSTANTS
      IMS=IMSX
      IVS=IVSX
      IES=IESX
      NM=NMX
      NV=NVX
      NE=NEX
      IVLIM=(NV-1)*IVS
      ILIM=(NE-1)*IES
      LSTEP=IMS
C  LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
      LX=1
      MX=NM
      IFAC=TABLES(-1)
      IF(IFAC.GT.3)THEN
        IERR=TBERR
        RETURN
      ENDIF
C..   RADIX 2 LOOP
200   CONTINUE
      DO 210 IM=1,TABLES(-14)
        MX=MX/2
        NUSTEP=LX*LSTEP
        LLIM=NUSTEP-LSTEP
        LD2LIM=NUSTEP/2-LSTEP
        MSTEP=NUSTEP*2
        MLIM=NM*LSTEP-MSTEP
        CALL MFFTA7(C,TABLES(0))
        LX=LX+LX
210   CONTINUE
      IF(IFAC.EQ.1)RETURN
C..   RADIX 3 LOOP
300   CONTINUE
      DO 310 IM=1,TABLES(-13)
        MX=MX/3
        NUSTEP=LX*LSTEP
        LLIM=NUSTEP-LSTEP
        MSTEP=NUSTEP*3
        MLIM=NM*LSTEP-MSTEP
        CALL MFFTB7(C,TABLES(0))
        LX=LX*3
310   CONTINUE
      IF(IFAC.EQ.2)RETURN
C..   RADIX 5 LOOP
500   CONTINUE
      DO 510 IM=1,TABLES(-12)
        MX=MX/5
        NUSTEP=LX*LSTEP
        LLIM=NUSTEP-LSTEP
        MSTEP=NUSTEP*5
        MLIM=NM*LSTEP-MSTEP
        CALL MFFTC7(C,TABLES(0))
        LX=LX*5
510   CONTINUE
      END
      SUBROUTINE MFFTIS(C,IMSX,IVSX,IESX,NMX,NVX,NEX,TABLES,IERR)
C   PURPOSE:
C       THE SAME AS MFFTIM.  IT IS A VARIANT OF MFFTIM, OPTIMIZED FOR
C       MAXIMUM PERFORMANCE ON MATRICES WHOSE FIRST DIMENSION IS
C       LESS THAN 64.
C       IT REQUIRES THAT TABLES HAS BEEN PREPARED BY A CALL TO MFFTP
C       WITH ID .NE. 0 .
C  IMPLEMENTATION:
C       THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
C       "BUTTERFLY" TRANSFORMATION MFFT?9; PARAMETERS OF THE "BUTTERFLY"
C       ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
C  ARGUMENTS:
C    INPUT:
C       C : ARRAY TO BE TRANSFORMED.
C       IMSX,IVSX,IESX,NMX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE
C           OF C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED
C           ON OUTPUT
C       TABLES : ARRAY PREPARED BY MFFTP. IT  IS NOT CHANGED ON OUTPUT.
C                IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
C                IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
C    OUTPUT:
C       C : TRANSFORM OF THE ORIGINAL ARRAY; "BIT REVERSED" ORDER
C       IERR : ERROR CODE : =0 : SUCCESSFUL
C                         : =3 :  'TABLES' NOT CORRECTLY INITIALIZED
      COMPLEX C(*)
      INTEGER TABLES(-14:*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      SAVE IBASE
C  LOADING THE COMMON BLOCK : CONSTANTS
      IMS=IMSX
      IVS=IVSX
      IES=1
      NM=NMX
      NV=NVX
      NE=NEX

      IVLIM=(NV-1)*IVS
      LSTEP=IMS
C  LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
      LX=1
      MX=NM

      IBASE=(2+IMS)*NM
      IFAC=TABLES(-1)
      IF(IFAC.GT.3)THEN
        IERR=TBERR
        RETURN
      ENDIF
C..   RADIX 2 LOOP
200   CONTINUE
      DO 210 IM=1,TABLES(-14)
        MX=MX/2
        NUSTEP=LX*LSTEP
        ILIM=NUSTEP-1
        MSTEP=NUSTEP*2
        MLIM=NM*LSTEP-MSTEP
        CALL MFFTA9(C,TABLES(2*IBASE))
        LX=LX+LX
        IBASE=IBASE+NUSTEP
210   CONTINUE
      IF(IFAC.EQ.1)RETURN
C..   RADIX 3 LOOP
300   CONTINUE
      DO 310 IM=1,TABLES(-13)
        MX=MX/3
        NUSTEP=LX*LSTEP
        ILIM=NUSTEP-1
        MSTEP=NUSTEP*3
        MLIM=NM*LSTEP-MSTEP
        CALL MFFTB9(C,TABLES(2*IBASE))
        LX=LX*3
        IBASE=IBASE+NUSTEP*2
310   CONTINUE
      IF(IFAC.EQ.2)RETURN
C..   RADIX 5 LOOP
500   CONTINUE
      DO 510 IM=1,TABLES(-12)
        MX=MX/5
        NUSTEP=LX*LSTEP
        ILIM=NUSTEP-1
        MSTEP=NUSTEP*5
        MLIM=NM*LSTEP-MSTEP
        CALL MFFTC9(C,TABLES(2*IBASE))
        LX=LX*5
        IBASE=IBASE+NUSTEP*4
510   CONTINUE
      END
      SUBROUTINE MFFTIV(C,IVSX,IESX,NVX,NEX,TABLES,IERR)
C
C   PURPOSE:
C       THIS SUBROUTINE PERFORMS AN INVERSE FOURIER TRANSFORM ALONG
C       ONE DIMENSION OF A 2-DIMENSIONAL MATRIX, USING THE
C       COOLEY-TUKEY ALGORITHM.
C       THE INPUT MATRIX IS ASSUMED TO HAVE BEEN SUBJECTED
C        TO A "BIT-REVERSAL" REORDERING ( THROUGH MFFTOV OR
C       BECAUSE IT IS THE OUTPUT OF MFFDV AND HAS NOT BEEN REORDERED)
C       THE SEQUENCE TO BE TRANSFORMED IS C[IVSX,NVX], WHOSE COMPONENTS
C       ARE THE VECTORS C(M)[IESX,NEX].
C       SEE REF.[1] FOR NOTATIONS.
C  EXAMPLE:
C       LET C BE A 2-D MATRIX C(N1,N2) DECLARED VIA
C                 DIMENSION C(ID,N2)
C       WITH ID.GE.N1.
C       THEN THE IDFT ALONG THE FIRST DIMENSION IS OBTAINED BY
C                 CALL MFFTIV(C,1,ID,N1,N2,TABLES,IERR)
C       THE IDFT ALONG THE SECOND DIMENSION IS OBTAINED BY
C                 CALL MFFTIV(C,ID,1,N2,N1,TABLES,IERR)
C  IMPLEMENTATION:
C       THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
C       "BUTTERFLY" TRANSFORMATION MFFT?5; PARAMETERS OF THE "BUTTERFLY"
C       ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
C  ARGUMENTS:
C       C : ARRAY TO BE TRANSFORMED.
C       IVSX,IESX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE OF
C           C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED ON
C           OUTPUT;
C       TABLES : ARRAY PREPARED BY MFFTP. IT  IS NOT CHANGED ON OUTPUT.
C                IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
C                IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
C  OUTPUT:
C       C : TRANSFORM OF THE ORIGINAL ARRAY; "BIT REVERSED" ORDER
C       IERR : ERROR CODE : =0 : SUCCESSFUL
C                         : =3 :  'TABLES' NOT CORRECTLY INITIALIZED
      COMPLEX C(*)
      INTEGER TABLES(-14:*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
C  LOADING THE COMMON BLOCK : CONSTANTS
      IVS=IVSX
      IES=IESX
      NV=NVX
      NE=NEX
      ILIM=(NE-1)*IES
      LSTEP=IVS
C  LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
      LX=1
      MX=NV
C  SELECT THE HIGHEST FACTOR OF NV
      IFAC=TABLES(-1)
      IF(IFAC.GT.3)THEN
        IERR=TBERR
        RETURN
      ENDIF
C..  RADIX 2 LOOP
200   CONTINUE
      DO 210 IM=1,TABLES(-14)
        MX=MX/2
        NUSTEP=LX*LSTEP
        LLIM=NUSTEP-LSTEP
        LD2LIM=NUSTEP/2-LSTEP
        MSTEP=NUSTEP*2
        MLIM=NV*LSTEP-MSTEP
        CALL MFFTA5(C,TABLES(0))
        LX=LX+LX
210   CONTINUE
      IF(IFAC.EQ.1)RETURN
C..  RADIX 3 LOOP
300   CONTINUE
      DO 310 IM=1,TABLES(-13)
        MX=MX/3
        NUSTEP=LX*LSTEP
        LLIM=NUSTEP-LSTEP
        MSTEP=NUSTEP*3
        MLIM=NV*LSTEP-MSTEP
        CALL MFFTB5(C,TABLES(0))
        LX=LX*3
310   CONTINUE
      IF(IFAC.EQ.2)RETURN
C..   RADIX 5 LOOP
500   CONTINUE
      DO 510 IM=1,TABLES(-12)
        MX=MX/5
        NUSTEP=LX*LSTEP
        LLIM=NUSTEP-LSTEP
        MSTEP=NUSTEP*5
        MLIM=NV*LSTEP-MSTEP
        CALL MFFTC5(C,TABLES(0))
        LX=LX*5
510   CONTINUE
      END
      SUBROUTINE MFFTOM(C,IMS,IVS,IES,NM,NV,NE,INDEX,ITEMP)
C
C   PURPOSE :
C     THIS ROUTINE PERFORMS A REORDERING OF A VECTOR-OF-2 VECTORS
C     OF COMPLEX C[IMS,NM [IVS,NV [IES,NE]]], ACCORDING TO A
C     PERMUTATION INDEX "INDEX".
C     SEE REF.[1] FOR NOTATIONS, AND COMMENTS TO MFFTDM ABOVE.
C  ARGUMENTS
C     C : VECTOR-OF-2 VECTORS TO BE REORDERED
C     IMS,IVS,IES,NM,NV,NE: THESE ARGUMENTS DESCRIBE THE STRUCTURE OF
C     C, ACCORDING TO THE ABOVE DEFINITION;
C     INDEX: INTEGER ARRAY , CONTAINING THE PERMUTATION INDEX  ;
C            IT IS NM ELEMENTS LONG ; PREPARED BY MFFTP.
C     IWORK: INTEGER ARRAY, OF LENGTH AT LEAST NM, USED AS WORKSPACE;
C
      INTEGER INDEX(0:*),ITEMP(0:*)
      COMPLEX C(0:IMS-1,0:*),T

      I3LIM=(NV-1)*IVS
      JLIM=(NE-1)*IES
      DO 1 I=0,NM-1
1     ITEMP(I)=INDEX(I)

      DO 4 I=1,NM-3
2       IF(ITEMP(I).NE.I)THEN
          IDEST=ITEMP(I)
          DO 3 I3=0,I3LIM,IVS
          DO 3 J=I3,I3+JLIM,IES
            T=C(J,I)
            C(J,I)=C(J,IDEST)
            C(J,IDEST)=T
3         CONTINUE
          ITEMP(I)=ITEMP(IDEST)
          ITEMP(IDEST)=IDEST
        GOTO 2
        ENDIF
4     CONTINUE
      END
      SUBROUTINE MFFTOV(C,IVS,IES,NV,NE,INDEX,ITEMP)
C
C   PURPOSE :
C     THIS ROUTINE PERFORMS A REORDERING OF A VECTOR-OF-VECTORS
C     OF COMPLEX C[IVS,NV [IES,NE]], ACCORDING TO A
C     PERMUTATION INDEX "INDEX".
C     SEE REF.[1] FOR NOTATIONS, AND COMMENTS TO MFFTDV.
C  ARGUMENTS
C     C : VECTOR-OF-VECTORS TO BE REORDERED
C     IVS,IES,NV,NE: THESE ARGUMENTS DESCRIBE THE STRUCTURE OF
C     C, ACCORDING TO THE ABOVE DEFINITION;
C     INDEX: INTEGER ARRAY , CONTAINING THE PERMUTATION INDEX ;
C            IT IS NV ELEMENTS LONG; PREPARED BY MFFTP.
C     IWORK: INTEGER ARRAY, OF LENGTH AT LEAST NV, USED AS WORKSPACE;
C
      INTEGER INDEX(0:NV-1),ITEMP(0:NV-1)
      COMPLEX C(IVS,0:*),T

      NEIES=NE*IES
      DO 1 I=0,NV-1
1     ITEMP(I)=INDEX(I)

      DO 4 I=1,NV-3
2       IF(ITEMP(I).NE.I)THEN
          IDEST=ITEMP(I)
          DO 3 J=1,NEIES,IES
            T=C(J,I)
            C(J,I)=C(J,IDEST)
            C(J,IDEST)=T
3         CONTINUE
          ITEMP(I)=ITEMP(IDEST)
          ITEMP(IDEST)=IDEST
        GOTO 2
        ENDIF
4     CONTINUE
      END
      SUBROUTINE MFFTP(N,W,ID,IERR)
C
C     THIS SUBROUTINE PREPARES TABLES FOR USE BY THE MFFT
C     ROUTINES. THE PARAMETERS ARE
C
C     N : IS THE ORDER OF THE TRANSFORM;
C
C     W : ARRAY OF LENGTH AT LEAST 4*N+14  WORDS IF ID=0,
C      AND 4*N*(ID+1)+14 IF ID > 0 ; THIS ARRAY IS FILLED
C      BY THE PRESENT ROUTINE, AND SHOULD NOT BE MODIFIED BY THE
C      USER; IT IS REQUIRED BY ALL THE OPERATING ROUTINES;
C      THE ROUTINES MFFTIS AND MFFTDS
C      REQUIRE THAT W HAS BEEN FILLED BY A CALL TO MFFTP WITH
C      ID .GT. 0 ; ALL THE ROUTINES DO NOT MODIFY THE CONTENTS
C      OF W ;
C      WARNING: DIFFERENT PORTIONS OF W ARE HANDLED AS COMPLEX OR
C               INTEGER VARIABLES BY DIFFERENT ROUTINES.
C
C     ID : IF ID .EQ. 0 THE TABLES ARE SET FOR A NORMAL
C      TRANSFORM; IF ID .GT. 0 THE TABLES ARE SET FOR BOTH
C      A NORMAL AND A "SPECIAL" TRANSFORM (I.E. OPTIMIZED FOR SMALL
C      FIRST DIMENSION DATA ARRAYS, SEE REF.[2]); IN THIS CASE IT
C      SHOULD BE EQUAL TO THE FIRST DIMENSION OF THE ARRAY TO BE
C      TRANSFORMED, AS DECLARED IN THE CALLING PROGRAM.
C
C     IERR : ERROR CODE : =0 : SUCCESSFUL
C                       : =2 : FACTORIZATION ERROR
C
C
C************************************************************
C
C      REFERENCE INFORMATION : LAYOUT OF W
C
C  IN ALL CASES
C
C WORD ADDRESS      TYPE    N. OF EL. WORD LENGTH   PURPOSE
C
C  0                INTEGER      14      14     FACTORIZATION OF N
C
C  14               COMPLEX     N       2*N     EXP(I*P/(2*N)*K),K=
C                                               0,N-1
C  3*N+14           INTEGER     N        N      PERMUTATION INDEX:KOFI
C  2*N+14           INTEGER     N        N      PERMUTATION INDEX:IOFK
C
C  ONLY IF ID .GT. 0
C
C  4*N+14           COMPLEX    N*ID     2*N*ID  TABLES FOR MFFTDS
C  4*N+2*N*ID+14    COMPLEX    N*ID     2*N*ID  TABLES FOR MFFTIS
C
C
C
C******************************************************************
      COMPLEX W(-7:*)

C...  FACTORIZATION OF N IN W(-7)..W(-1)
      CALL MFFTP1(W(-7),N,IERR)
      IF(IERR.NE.0)RETURN
C
C     PREPARATION OF PERMUTATION INDEXES IN W(N)..W(2*N-1)
C     WARNING : W(0)..W(N-1) USED AS A WORK SPACE
      CALL MFFTP2(W(N),W(0),W(-7),N)
C
C     PREPARATION OF PHASE FACTOR TABLE IN W(0)..W(N-1)
      W(0)=1.0
      PI2DN=ATAN(1.)*8./N
      DO 1  I=1,N-1
        W(I)=CMPLX(COS(PI2DN*I),SIN(PI2DN*I))
1     CONTINUE
C
C     IF TABLES FOR SPECIAL TRANSFORM ARE REQUESTED
      IF(ID.GT.0)THEN
        CALL MFFTP4(W(0),W(2*N),W(-7),N,ID)
      ENDIF
C
      END
      SUBROUTINE MFFTP1(W,NM,IERR)
C
C   PURPOSE :
C     FACTORIZATION OF NM STORING THE POWERS OF EACH FACTOR IN W
C     ( NM = 2**W(1)*3**W(2)*... )
C     THE MAXIMUM FACTOR FOUND  IS STORED IN W(14)
C
      PARAMETER (MAXFAC=3)
      INTEGER W(14),FACTOR(MAXFAC)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      DATA (FACTOR(I),I=1,MAXFAC)/2,3,5/
      N=NM
      DO 100 I=1,MAXFAC
      W(I)=0
 10   IF(MOD(N,FACTOR(I)).EQ.0)THEN
      W(I)=W(I)+1
      N=N/FACTOR(I)
      GOTO 10
      ENDIF
      IF(N.EQ.1)GOTO 200
100   CONTINUE
      IERR=FACERR
      RETURN
200   W(14)=I
      END
      SUBROUTINE MFFTP2(INDX,I2,IW,NM)
C
C     THIS ROUTINE COMPUTES TWO INDEX TABLES FOR THE
C     PERMUTATION DUE TO REPRESENTATION INVERSION ("BIT REVERSAL")
C     THE INDEX TABLES ARE STORED ONE AFTER THE OTHER, AND ARE
C     RECIPROCAL.
C
C     WARNING: BY "BIT REVERSAL" WE MEAN THE SHUFFLING OF INDEXES
C              AS REQUIRED BY IN-PLACE FFT ALGORITHMS, REGARDLESS
C              OF WHAT IS THEIR RADIX (I.E. 2,3,5).
C              THE SHUFFLING IS EFFECTIVELY AN 'INTEGER BIT REVERSAL'
C              ONLY IN CASE OF RADIX-2 ALGORITHMS.
C
      INTEGER INDX(0:*),I2(0:*),IW(14)
      PARAMETER (MAXFAC=3)
      INTEGER FACTOR(MAXFAC)
      DATA (FACTOR(I),I=1,MAXFAC)/2,3,5/
      DO 1 I=0,NM-1
        INDX(I)=I
1     CONTINUE
      IFAC=IW(14)
      LX=1
      MX=NM
      DO 30 IFAC=1,IFAC
        DO 10 I=1,IW(IFAC)-1,2
          MX=MX/FACTOR(IFAC)
          CALL MFFTP3(INDX,I2,MX,FACTOR(IFAC),LX)
          LX=LX*FACTOR(IFAC)
          MX=MX/FACTOR(IFAC)
          CALL MFFTP3(I2,INDX,MX,FACTOR(IFAC),LX)
          LX=LX*FACTOR(IFAC)
10      CONTINUE

C...  IF W(IFAC)ODD,THEN
        IF (I.EQ.IW(IFAC))THEN
          MX=MX/FACTOR(IFAC)
          CALL MFFTP3(INDX,I2,MX,FACTOR(IFAC),LX)
          LX=LX*FACTOR(IFAC)
          DO 20 I=0,NM-1
            INDX(I)=I2(I)
20        CONTINUE
        ENDIF
30    CONTINUE
C...     INVERSE PERMUTATION
         DO 40 I=0,NM-1
           I2(I)=INDX(I)
           INDX(I+NM)=I
40       CONTINUE
         DO 59 I=1,NM-3
51         IF(I2(I).NE.I)THEN
             IDEST=I2(I)
             IT=INDX(I+NM)
             INDX(I+NM)=INDX(IDEST+NM)
             INDX(IDEST+NM)=IT
             I2(I)=I2(IDEST)
             I2(IDEST)=IDEST
           GOTO 51
           ENDIF
59       CONTINUE

      END
      SUBROUTINE MFFTP3(INDX,I1,MX,NX,LX)
C
C     THIS SUBROUTINE PERFORMS A "BIT REVERSAL" PERMUTATION
C
      INTEGER INDX(MX,NX,LX),I1(MX,LX,NX)


      DO 1 NU=1,NX
      DO 1 MU=1,MX
      DO 1 LAM=1,LX
        I1(MU,LAM,NU)=INDX(MU,NU,LAM)
1     CONTINUE
      END
      SUBROUTINE MFFTP4(EXPTAB,SPETAB,FACTAB,N,N1)
C
C     THIS SUBROUTINE BUILDS THE TWIDDLE FACTOR TABLES FOR USE
C     OF MFFT?S ROUTINES (SPECIAL OPTIMIZATION FOR SMALL
C     DATA MATRICES); IT MUST BE USED BY MFFTP ONLY.
C
C     PARAMETERS:
C     EXPTAB: TWIDDLE FACTOR TABLE
C     SPETAB: SPECIAL TWIDDLE FACTOR TABLE
C     FACTAB: FACTORIZATION OF N
C     N: ORDER OF THE TRANSFORM
C     N1: FIRST DIMENSION OF THE ARRAY TO BE TRANSFORMED (.GE.N)
C
      INTEGER MAXFAC
      PARAMETER(MAXFAC=3)
      INTEGER FACTAB(14)
      COMPLEX EXPTAB(0:*),SPETAB(0:*)
      INTEGER FACTOR(MAXFAC)
      DATA FACTOR/2,3,5/

      MX=N
      LX=1
      J=0
      DO 50 IFACT=FACTAB(14),1,-1
        NX=FACTOR(IFACT)
        DO 40 IPOW=1,FACTAB(IFACT)
          MX=MX/NX
          DO 30 NU=1,NX-1
            DO 20 MU=0,MX-1
              DO 10 I1=0,N1-1
                SPETAB(J)=CONJG(EXPTAB(MU*LX*NU))
                J=J+1
10            CONTINUE
20          CONTINUE
30        CONTINUE
          LX=LX*NX
40      CONTINUE
50    CONTINUE
      MX=N
      LX=1
      J=N*N1
      DO 100 IFACT=1,FACTAB(14)
        NX=FACTOR(IFACT)
        DO 90 IPOW=1,FACTAB(IFACT)
          MX=MX/NX
          DO 80 NU=1,NX-1
            DO 70 LAMBDA=0,LX-1
              DO 60 I1=0,N1-1
                SPETAB(J)=EXPTAB(MX*LAMBDA*NU)
                J=J+1
60            CONTINUE
70          CONTINUE
80        CONTINUE
          LX=LX*NX
90      CONTINUE
100   CONTINUE
      END
      SUBROUTINE MFFTRD(C,ISV,ISE,NV,NE,RW)
C
C   PURPOSE:
C       THIS ROUTINE PERFORMS THE POST-PROCESSING PHASE FOR
C       REAL 2-DIMENSIONAL DFT'S, ACCORDING TO FORMULA (2.7)
C       IN REF.[1].
C       POST-PROCESSING ACTS AFTER COMPUTING THE COMPLEX DFT
C       AND EVENTUAL REORDERING (CALLS TO MFFTDV AND MFFTOV).
C       IT APPLIES TO A VECTOR-OF-VECTORS-OF-COMPLEX
C               C[IVS,NV [IES,NE]].
C       SEE REF.[1] FOR NOTATIONS.
C
C   ARGUMENTS:
C      INPUT:
C       C : DATA ARRAY, OUTPUT FROM MFFTDV, MFFTOV; TO BE DECLARED
C                      REAL C(ISE*2,NE)
C           IN THE CALLING PROGRAM.
C     ISV : SEPARATION OF ELEMENTS IN A COLUMN OF C (USUALLY 1)
C     ISE : SEPARATION OF ELEMENTS IN A ROW OF C, DIVIDED BY 2
C      NV : NO. OF ELEMENTS TO BE PROCESSED IN A COLUMN OF C
C      NE : NO. OF ELEMENTS IN A ROW OF C, DIVIDED BY 2.
C      RW : COMPLEX ARRAY OF LENGHT AT LEAST NV; IT MUST
C           BE INITIALIZED BY A CALL TO MFFTRP; IT REMAINS
C           UNCHANGED IN OUTPUT.
C
C  OUTPUT : POST-PROCESSED ARRAY C
C
C
      COMPLEX C(0:ISV-1,0:*),RW(0:*),T1,T2
C
      IF (NV.GT.1) THEN
      DO 200 IV=1,(NV-1)/2
       DO 190 IE=0,(NE-1)*ISE,ISE
C
      T1=C(IE,IV)
      T2=C(IE,NV-IV)
      C(IE,IV)=((T1+CONJG(T2))+(RW(IV)*(T1-CONJG(T2))))*0.5
      C(IE,NV-IV)=(CONJG(T1+CONJG(T2))-CONJG(RW(IV)*(T1-CONJG(T2))))*0.5
C
 190   CONTINUE
 200  CONTINUE
C
C
      IF(2*IV.EQ.NV) THEN
        DO 210 IE=0,(NE-1)*ISE,ISE
        C(IE,IV)=CONJG(C(IE,IV))
 210  CONTINUE
C
      ENDIF
C
       ENDIF

      DO 300 IE=0,(NE-1)*ISE,ISE
        T1=C(IE,0)
        C(IE,0)=(REAL(T1)+AIMAG(T1))
        C(IE,NV)=(REAL(T1)-AIMAG(T1))
 300  CONTINUE
C
C
      END
      SUBROUTINE MFFTRI(C,ISV,ISE,NV,NE,RW)
C
C   PURPOSE:
C       THIS ROUTINE PERFORMS THE PRE-PROCESSING PHASE FOR
C       REAL 2-DIMENSIONAL INVERSE DFT'S (SEE FORMULA (2.7)
C       IN REF.[1]).
C       PRE-PROCESSING ACTS ON SEQUENTAL DATA BEFORE A CALL TO
C       MFFTOV AN BEFORE COMPUTING THE IDFT (A CALL TO MFFTIV)
C       AND EVENTUAL REORDERING (CALLS TO MFFTIV AND MFFTOV).
C       IT APPLIES TO A VECTOR-OF-VECTORS-OF-COMPLEX
C               C[IVS,NV [IES,NE]].
C       SEE REF.[1] FOR NOTATIONS.
C
C   ARGUMENTS:
C      INPUT:
C       C : DATA ARRAY, OUTPUT FROM MFFTDV, MFFTOV; TO BE DECLARED
C                      REAL C(ISE*2,NE)
C           IN THE CALLING PROGRAM.
C     ISV : SEPARATION OF ELEMENTS IN A COLUMN OF C (USUALLY 1)
C     ISE : SEPARATION OF ELEMENTS IN A ROW OF C, DIVIDED BY 2
C      NV : NO. OF ELEMENTS TO BE PROCESSED IN A COLUMN OF C
C      NE : NO. OF ELEMENTS IN A ROW OF C, DIVIDED BY 2.
C      RW : COMPLEX ARRAY OF LENGHT AT LEAST NV; IT MUST
C           BE INITIALIZED BY A CALL TO MFFTRP; IT REMAINS
C           UNCHANGED IN OUTPUT.
C
C  OUTPUT : POST-PROCESSED ARRAY C
C
C
      COMPLEX C(0:ISV-1,0:*),RW(0:*),T1,T2
C
C
      IF(NV.GT.1) THEN
C
      DO 200 IV=1,(NV-1)/2
        DO 190 IE=0,(NE-1)*ISE,ISE
        T1=C(IE,IV)
        T2=C(IE,NV-IV)
C
        C(IE,IV)=(T1+CONJG(T2))+(CONJG(RW(IV))*(T1-CONJG(T2)))
        C(IE,NV-IV)=CONJG(T1+CONJG(T2))-CONJG(CONJG(RW(IV))*(T1-
     $              CONJG(T2)))


 190    CONTINUE
 200   CONTINUE
C
C
        IF(2*IV.EQ.NV) THEN
C
        DO 210 IE=0,(NE-1)*ISE,ISE
        C(IE,IV)=2*CONJG(C(IE,IV))
 210   CONTINUE
C
C
      ENDIF
C
      ENDIF
C
C
      DO 220 IE=0,(NE-1)*ISE,ISE
         RP=REAL(C(IE,0))+REAL(C(IE,NV))
         RM=REAL(C(IE,0))-REAL(C(IE,NV))
         C(IE,0)=CMPLX(RP,RM)
 220  CONTINUE
C
C
      END
      SUBROUTINE MFFTRP(NN,RW)
C
C  THIS ROUTINE PREPARES THE PHASE FACTOR TABLES TO BE USED IN THE
C  POST-PROCESSING PHASE IN REAL TRANSFORMS, I.E. R2FFT AND R3FFT.
C
C  THE DEFINITION OF PHASE FACTORS WE USE IS THE FOLLOWING
C
C             RW(IV) = -I*EXP(-I*2*PI/NN*IV)
C
C  WHERE THE SYMBOLS ARE:
C
C  I: IMAGINARY UNIT
C  NN : NUMBER OF REAL ELEMENTS
C  IV: RUNNING INDEX FROM 0 TO NN-1
C
C
C
C
      COMPLEX RW(0:*)
C
C
      PI2=8*ATAN(1.)
      PHASE=PI2/NN
      RW(0)=(0.,-1.)
C
C
      DO 10 I=1,NN-1
      RW(I)=CMPLX(-SIN(I*PHASE),-COS(I*PHASE))
 10   CONTINUE
C
C
      END
         SUBROUTINE MFFTZ0(W1,ISE1,NE,W2,ISE2)
C*** PURPOSE : VECTOR COPY ROUTINE
C    PARAMETERS :
C        W1 : VECTOR TO BE COPIED
C        ISE1: STRIDE OF W1
C        NE : NUMBER OF ELEMENTS
C        W2 : DESTINATION VECTOR
C        ISE2: STRIDE OF W2
         COMPLEX W1(0:*),W2(0:*)
         J=0
         DO 1 I=0,(NE-1)*ISE1,ISE1
         W2(J)=W1(I)
         J=J+ISE2
1        CONTINUE
         END
ID       STATION
JOB,JN=TK3FT,MFL,T=10.
ACCOUNT,AC=TK3TSZ13,US=TK3TSZ13,APW=TTT5NHPR,UPW=TTT5NHPR.
ACCESS,DN=NOTE,OWN=TK3TSZ13.
NOTE,TEXT='USR=(vito,cray)'.
CFT,OPT=FASTMD:BTREG,ON=FZ,OFF=TC.
ACCESS,DN=FFTE.
LDR,LIB=FFTE.
EXIT.
DUMPJOB..
DEBUG,BLOCKS.
/EOF
         PROGRAM TEST
C
C***************** TESTING PROGRAM OF THE COMPLEX 2-DIM FFT ROUTINE
C                                     C2FFT
      INTEGER NPRNT,YPRNT
      PARAMETER(NPRNT=0,YPRNT=1)
      PARAMETER(MAXTEST=1)
      PARAMETER(NSIZE=513)
C
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      INTEGER SIZERR,DFFTERR,IFFTERR
      PARAMETER (SIZERR=10,DFFTERR=100,IFFTERR=200)
C
      PRINT '(////,3X,''TEST OF COMPLEX TWO DIMENSIONAL MIXED-RADIX FFT
     $ROUTINES'',/)'
      PRINT'(T4,''INPUT PARAMETERS ARE :'',//,T6,''ID'',T14,
     $       ''NL'',T23,''NM'',T31,''II'',T39,''JJ'',T46,
     $       ''IORD'',/)'
C
         DO 1 NTEST=1,MAXTEST
C
C***************** GENERATE MATRIX DIMENSIONS NL,NM
C                  AND ALL THE OTHER INPUT PARAMETERS
C
           CALL GENERA(ID,NL,NM,II,JJ,IORD,NSIZE)
C
           PRINT'(6(I7,1X))',ID,NL,NM,II,JJ,IORD
C
C***************** CALL TESTING ROUTINE
C
           CALL TESTFT(ID,NL,NM,NN,II,JJ,IORD,IEXIT,NPRNT)
C
           IF(IEXIT.NE.0)THEN
             IF(IEXIT.EQ.IFFTERR.OR.IEXIT.EQ.DFFTERR.OR.
     $          IEXIT.EQ.TBERR                           )THEN
               PRINT '(//,9X,I7)',IEXIT
C
C***************** IN CASE OF ERROR, PRINT THE CONTENT OF DATA MATRIX
C                  ONLY IF SMALL ENOUGH (I.E. MAX(NL,NM) LESS THAN 32)
C
               IF(MAX(NL,NM).LE.32)THEN
                 CALL TESTFT(ID,NL,NM,NN,II,JJ,IORD,IEXIT,YPRNT)
               ENDIF
               STOP
             ENDIF
           ENDIF
C
1        CONTINUE
C
C
         END
C
C
      SUBROUTINE GENERA(ID,NL,NM,II,JJ,IORD,NSIZE)
C
C***************** THIS ROUTINE PERFORMS A RANDOM GENERATION OF THE
C                  INPUT DATA DIMENSIONS NL,NM, AS WELL AS OF ALL
C                  THE CALLING PARAMETERS OF THE COMPLEX 2-DIM FFT
C                  ROUTINE C2FFT
C
C
      REAL MAXDIM
      PARAMETER (MAXDIM=850)
C
       L2MAX=(LOG(REAL(NSIZE*2)))/LOG(2.0)
       L3MAX=(LOG(REAL(NSIZE*3)))/LOG(3.0)
       L5MAX=(LOG(REAL(NSIZE*5)))/LOG(5.0)
C
10    CONTINUE
      L2=RANF()*L2MAX
      L3=RANF()*L3MAX
      L5=RANF()*L5MAX
C
      M2=RANF()*L2MAX
      M3=RANF()*L3MAX
      M5=RANF()*L5MAX
C
      NL=2**L2*3**L3*5**L5
      IF(NL.EQ.1) NL=2
      NM=2**M2*3**M3*5**M5
      IF(NM.EQ.1) NM=2
C
C***************** GENERATED DIMENSIONS SHOULD NOT EXCEED
C                  MAXIMUM MEMORY STORAGE
C
      IF(NL*NM.GT.NSIZE*NSIZE) GO TO 10
      IF(NL.GT.MAXDIM.OR.NM.GT.MAXDIM) GO TO 10
C
      NPL=RANF()*2
      ID=NL+NPL
C
C***************** GENERATE RANDOM POSITIONS OF INPUT FREQUENCIES
C
      II=RANF()*NL
      JJ=RANF()*NM
C
      IORD=RANF()*2
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT PARAMETERS ARE FORCED TO
C                  ASSUME THE FOLLOWING VALUES..........
C
      NL=10
      NM=4
      ID=11
      II=1
      JJ=1
      IORD=1
C
      RETURN
      END
C
C
      SUBROUTINE TESTFT(ID,NL,NM,NN,II,JJ,IORD,IEXIT,NPRINT)
C
C***************** THIS ROUTINE PERFORMS THE TESTS ON C2FFT
C
      INTEGER SIZERR,DFFTERR,IFFTERR
      PARAMETER (SIZERR=10,DFFTERR=100,IFFTERR=200)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
C
      PARAMETER (NSIZE=513)
      COMPLEX A(NSIZE*NSIZE),B(NSIZE*NSIZE)
      INTEGER WM(4*NSIZE+14),WL(4*NSIZE+14)
      INTEGER IWORK(NSIZE)
C
        DATA EPS/1.E-7/
C
C***************** SET IEXIT TO 1 (MFFTP SHOULD RESET IT TO 0)
C
      IEXIT=1

C
C***************** FILL INPUT DATA MATRICES WITH
C                  MONOCHROMATIC SIGNALS
C
       DO 10 I=1,ID*NM
         A(I)=0
         B(I)=0
10    CONTINUE
      CALL FILL2(A,B,ID,NL,NM,II,JJ)
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT MATRIX IS NOW PRINTED
C
      PRINT '(///,20X,''INPUT DATA MATRIX'',/)'
       CALL SHOWRE2(A,ID,NL,NM)
C
C***************** INITIALIZE FFT PACKAGE
C
       CALL C2FFT(A,ID,NL,NM,WL,WM,0,IORD,IWORK,IEXIT)
       IF(IEXIT.NE.0)RETURN
C
C***************** CALCULATE DIRECT FOURIER TRANSFORM OF A
C
           ISIG=-1
            CALL C2FFT(A,ID,NL,NM,WL,WM,ISIG,IORD,IWORK,IEXIT)
            IF(IEXIT.NE.0)RETURN
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIX IS NOW
C                  PRINTED
C
      PRINT '(/,10X,''DIRECT FOURIER TRANSFORM OF INPUT MATRIX'',/)'
       CALL SHOWRE2(A,ID,NL,NM)
C
C
C***************** COMPARE THE RESULT WITH KNOWN EXACT VALUE
C                  A 2-D DELTA IS SUBTRACTED FROM THE RESULT
C                  THE POSITION OF DELTA IS (II1,JJ1)
C
      IF(IORD.NE.0)THEN
        II1=II
        JJ1=JJ
      ELSE
        II1=IPERM(II,WL,NL)
        JJ1=IPERM(JJ,WM,NM)
      ENDIF
C
      CALL MINDE2(A,ID,NL,NM,II1,JJ1,REAL(NL*NM))
      IF(NPRINT.NE.0) CALL SHOWRE2(A,ID,NL,NM)
C
C***************** THE ERROR IS MAX(ABS(A))
C
          ERROR=XNOR2(A,ID,NL,NM)/REAL(NL*NM)
C
C***************** IF THE ERROR IS LARGE, THEN AN ERROR
C                  MESSAGE IS ISSUED, AND THE TEST IS STOPPED
C
       IF(ERROR.GT.EPS) THEN
          IEXIT=DFFTERR
          RETURN
       ENDIF
C
C***************** ... OTHERWISE MATRIX A IS RESTORED
C                      (I.E., A DELTA IS NOW ADDED)
C
      CALL PLSDE2(A,ID,NL,NM,II1,JJ1,REAL(NL*NM))
C
C***************** INVERSE FFT OF A IS NOW COMPUTED
C
       ISIG=1
       CALL C2FFT(A,ID,NL,NM,WL,WM,ISIG,IORD,IWORK,IEXIT)
       IF(IEXIT.NE.0)RETURN
C
       IF(NPRINT.NE.0) CALL SHOWRE2(A,ID,NL,NM)
C
C***************** NORMALIZE AFTER TWO TRANSFORMS
C
        CALL TMSCO2(A,ID,NL,NM,1.0/REAL(NL*NM))
C
        IF(NPRINT.NE.0) CALL SHOWRE2(A,ID,NL,NM)
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIX IS NOW
C                  PRINTED
C
      PRINT '(/,10X,''DIRECT PLUS INVERSE FOURIER TRANSFORM OF DATA'',
     $        /)'
       CALL SHOWRE2(A,ID,NL,NM)
       PRINT '(/////////)'
C
C****************** CALCULATE THE ERROR ON THE INVERSE TRANSFORM
C
      CALL MINMAT2(A,A,B,ID,NL,NM)
      ERROR=XNOR2(A,ID,NL,NM)
C
C****************** TEST ON THE ERROR AFTER THE INVERSE TRANSFORM
C
       IF(ERROR.GT.EPS) THEN
         IEXIT=IFFTERR
         RETURN
       ENDIF
C
C****************** THE TEST HAS NOW BEEN SATISFACTORILY
C                   COMPLETED.  EVENTUALLY PRINT HERE
C                   HOW HAPPY YOU ARE....
C
        IEXIT=0
C
        END
C
      INTEGER FUNCTION IPERM(I,INDEX,NELEM)
C
C***************** THIS FUNCTION CALCULATES THE POSITION OF
C                  INDEXES WHEN ORDERING IS NOT REQUESTED
C
      INTEGER INDEX(-14:*)
      IPERM=INDEX(3*NELEM+I)
      END
C
      SUBROUTINE FILL2(A,B,ID,NL,NM,II,JJ)
C
C***************** FILL INPUT DATA MATRIX WITH MONOCHROMATIC SIGNAL
C
      COMPLEX A(0:ID-1,0:NM-1),B(0:ID-1,0:NM-1),CI,CJ
C
      PI2=ATAN(1.)*8
      PI2L=PI2/NL
      PI2M=PI2/NM
C
      DO 10 I=0,NL-1
      CI=CMPLX(COS(PI2L*I*II),SIN(PI2L*I*II))
      DO 11 J=0,NM-1
      CJ=CMPLX(COS(PI2M*J*JJ),SIN(PI2M*J*JJ))
C
      A(I,J)=CI*CJ
      B(I,J)=A(I,J)
C
  11  CONTINUE
  10  CONTINUE
      END
C
      SUBROUTINE SHOWRE2(A,ID,NL,NM)
C
C***************** PRINT MATRIX CONTENT
C
      COMPLEX A(ID,NM)
C
      DO 100 JJ=1,NL
 100  WRITE(6,'(1X,16F8.5)') (A(JJ,K),K=1,NM)
      END
C
      SUBROUTINE TMSCO2(A,ID,NL,NM,X)
C
C***************** MULTIPLY MATRIX ELEMENTS BY X
C
      COMPLEX A(ID,NM)
      REAL X
      DO 1 I=1,NM
      DO 1 J=1,ID
      A(J,I)=A(J,I)*X
 1    CONTINUE
      END
C
      FUNCTION XNOR2(A,ID,NL,NM)
C
C***************** CALCULATE MAXIMUM NORM OF MATRIX A
C
      COMPLEX A(ID,NM)
      XNOR2=0.
      DO 1 L=1,NL
      DO 1 M=1,NM
      XNOR2=MAX(XNOR2,ABS(REAL(A(L,M)))+ABS(AIMAG(A(L,M))))
 1    CONTINUE
      END
C
      SUBROUTINE MINDE2(A,ID,NL,NM,II,JJ,X)
C
C***************** NUMBER X IS SUBTRACTED FROM A(II,JJ)
C
      COMPLEX A(0:ID-1,0:NM-1)
      A(II,JJ)=A(II,JJ)-X
      END
      SUBROUTINE MINMAT2(A,B,C,ID,NL,NM)
C
C***************** MATRIX C IS SUBTRACTED FROM MATRIX B
C                  AND STORED IN A
C
      COMPLEX A(ID,NM),B(ID,NM),C(ID,NM)
      DO 1 L=1,NL
      DO 1 M=1,NM
      A(L,M)=B(L,M)-C(L,M)
 1    CONTINUE
      END
      SUBROUTINE PLSDE2(A,ID,NL,NM,II,JJ,X)
C
C***************** NUMBER X IS ADDED TO A(II,JJ)
C
      COMPLEX A(0:ID-1,0:NM-1)
      A(II,JJ)=A(II,JJ)+X
      END
ID       STATION
JOB,JN=TK3FT,MFL,T=10.
ACCOUNT,AC=TK3TSZ13,US=TK3TSZ13,APW=TTT5NHPR,UPW=TTT5NHPR.
ACCESS,DN=NOTE,OWN=TK3TSZ13.
NOTE,TEXT='USR=(vito,cray)'.
CFT,OPT=FASTMD:BTREG,ON=FZ,OFF=TC.
ACCESS,DN=FFTE.
LDR,LIB=FFTE.
EXIT.
DUMPJOB..
DEBUG,BLOCKS.
/EOF
         PROGRAM TEST
C
C***************** TESTING PROGRAM OF THE COMPLEX 3-DIM FFT ROUTINE
C                                       C3FFT
C
      INTEGER NPRNT,YPRNT
      PARAMETER(NPRNT=0,YPRNT=1)
      PARAMETER(MAXTEST=1)
      PARAMETER(NSIZE=65)
C
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACRR=2,TBERR=3)
      INTEGER SIZERR,DFFTERR,IFFTERR
      PARAMETER (SIZERR=10,DFFTERR=100,IFFTERR=200)
C
      PRINT '(/////)'
      PRINT '(6X,''TEST OF COMPLEX THREE DIMENSIONAL MIXED-RADIX FFT ROU
     $TINES'',//)'
      PRINT '(1X,''INPUT PARAMETERS ARE :'',/,
     $T6,''ID'',T14,''NL'',T22,''NM'',T30,''NN'',T38,''II''
     $,T46,''JJ'',T54,''KK'',T62,''IOPT'',T68,''IORD'')'
C
         DO 1 NTEST=1,MAXTEST
C
C***************** GENERATE MATRIX DIMENSIONS NL,NM,NN
C                  AND ALL THE OTHER INPUT PARAMETERS
C
           CALL GENERA(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,NSIZE)
           PRINT'(/9(I7,1X))',ID,NL,NM,NN,II,JJ,KK,IOPT,IORD
C
C
C***************** CALL TESTING ROUTINE
C
           CALL TESTFT(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,IEXIT,NPRNT)
C
           IF(IEXIT.NE.0)THEN
             IF(IEXIT.EQ.IFFTERR.OR.IEXIT.EQ.DFFTERR.OR.
     $          IEXIT.EQ.TBERR                           )THEN
               PRINT '(//,9X,I7)',IEXIT
C
C***************** IN CASE OF ERROR, PRINT THE CONTENT OF DATA MATRIX
C                  ONLY IF SMALL ENOUGH (I.E. MAX(NL,NM) LESS THAN 32)
C
               IF(MAX(NL,NM).LE.32)THEN
                 CALL TESTFT(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,IEXIT,YPRNT)
               ENDIF
               STOP
             ENDIF
           ENDIF
C
1        CONTINUE
C
         END
C
      SUBROUTINE GENERA(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,NSIZE)
C
C***************** THIS ROUTINE PERFORMS A RANDOM GENERATION OF THE
C                  INPUT DATA DIMENSIONS NL,NM,NN AS WELL AS OF ALL
C                  THE CALLING PARAMETERS OF THE COMPLEX 3-DIM FFT
C                  ROUTINE C3FFT
C
      REAL MAXDIM
      PARAMETER (MAXDIM=512)
C
       L2MAX=(LOG(REAL(NSIZE*2)))/LOG(2.0)
       L3MAX=(LOG(REAL(NSIZE*3)))/LOG(3.0)
       L5MAX=(LOG(REAL(NSIZE*5)))/LOG(5.0)
C
10    CONTINUE
C
      L2=RANF()*L2MAX
      L3=RANF()*L3MAX
      L5=RANF()*L5MAX
C
      M2=RANF()*L2MAX
      M3=RANF()*L3MAX
      M5=RANF()*L5MAX
C
      N2=RANF()*L2MAX
      N3=RANF()*L3MAX
      N5=RANF()*L5MAX
C
      NL=2**L2*3**L3*5**L5
      IF(NL.EQ.1) NL=2
      NM=2**M2*3**M3*5**M5
      IF(NM.EQ.1) NM=2
      NN=2**N2*3**N3*5**N5
      IF(NN.EQ.1) NN=2
C
      NPL=RANF()*2
      ID=NL+NPL
C
C***************** GENERATED DIMENSIONS SHOULD NOT EXCEED
C                  MAXIMUM STORAGE CAPACITY
C
      IF(ID*NM*NN.GT.NSIZE*NSIZE*NSIZE) GO TO 10
      IF(NL.GT.MAXDIM.OR.NM.GT.MAXDIM.OR.NN.GT.MAXDIM) GO TO 10
C
      II=RANF()*NL
      JJ=RANF()*NM
      KK=RANF()*NN
C
      IORD=RANF()*2
      IOPT=RANF()*2
C
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT PARAMETERS ARE FORCED TO
C                  ASSUME THE FOLLOWING VALUES..........
C
      NL=9
      NM=5
      NN=2
      ID=10
      II=1
      JJ=1
      KK=1
      IOPT=1
      IORD=1
C
      RETURN
      END
C
      SUBROUTINE TESTFT(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,IEXIT,NPRINT)
C
C***************** THIS ROUTINE PERFORMS THE TESTS ON C3FFT
C
      PARAMETER (NSIZE=65)
      INTEGER SIZERR,DFFTERR,IFFTERR
      PARAMETER (SIZERR=10,DFFTERR=100,IFFTERR=200)
C
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
C
      COMPLEX A(NSIZE*NSIZE*NSIZE),B(NSIZE*NSIZE*NSIZE)
      INTEGER WL(4*NSIZE+14),WM(4*NSIZE+14),WN(4*NSIZE+14)
      INTEGER IWORK(NSIZE)
C
      DATA EPS/1.E-7/
C
C***************** SET IEXIT TO 1 (MFFTP SHOULD RESET IT TO 0)
C
      IEXIT=1
C
C***************** TEST IF INPUT PARAMETERS EXCEED MAXIMUM
C                   WORK AREA SPACE; AN ERROR MESSAGE IS ISSUED
C                   AND EXECUTION IS STOPPED
C
      MXWORK=2*(NM+NM*ID*IOPT)
      IWSIZE=4*NSIZE+14
      IF(MXWORK.GT.IWSIZE) THEN
       IEXIT=SIZERR
       RETURN
      ENDIF
C
C***************** FILL INPUT MATRICES WITH
C                  MONOCHROMATIC SIGNALS
      DO 10 I=1,ID*NM*NN
         A(I)=0
         B(I)=0
10    CONTINUE
      CALL FILL(A,B,ID,NL,NM,NN,II,JJ,KK)
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT MATRIX IS NOW PRINTED
C
      PRINT '(///,26X,''INPUT DATA MATRIX'')'
      CALL SHOWRES(A,ID,NL,NM,NN)
C
C***************** INITIALIZE FFT PACKAGE
C
C
       CALL C3FFT(A,ID,NL,NM,NN,WL,WM,WN,IOPT,0,IORD,IWORK,IEXIT)
       IF(IEXIT.NE.0)RETURN
C
C***************** CALCULATE DIRECT FOURIER TRANSFORM OF A
C
       ISIG=-1
       CALL C3FFT(A,ID,NL,NM,NN,WL,WM,WN,IOPT,ISIG,IORD,IWORK,IEXIT)
       IF(IEXIT.NE.0)RETURN
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIX IS NOW PRINTED
C
      PRINT '(/15X,''DIRECT FOURIER TRANSFORM OF INPUT MATRIX'')'
      CALL SHOWRES(A,ID,NL,NM,NN)
C
C***************** COMPARE THE RESULT WITH KNOWN EXACT VALUE
C                  A 3-D DELTA IS SUBTRACTED FROM THE RESULT
C                  THE POSITION IF DELTA IS (II1,JJ1,KK1)
C
      IF(IORD.NE.0)THEN
        II1=II
        JJ1=JJ
         KK1=KK
      ELSE
        II1=IPERM(II,WL,NL)
        JJ1=IPERM(JJ,WM,NM)
         KK1=IPERM(KK,WN,NN)
      ENDIF
C
      CALL MINDEL(A,ID,NL,NM,NN,II1,JJ1,KK1,REAL(NL*NM*NN))
      IF(NPRINT.NE.0) CALL SHOWRES(A,ID,NL,NM,NN)
C
C***************** ERROR IS MAX(ABS(A))
C
       ERROR=XNORM(A,ID,NL,NM,NN)/REAL(NL*NM*NN)
C
C***************** IF THE ERROR IS LARGE ,THEN AN ERROR
C                  MESSAGE IS ISSUED, AND THE TEST IS STOPPED
C
       IF(ERROR.GT.EPS) THEN
          IEXIT=DFFTERR
          RETURN
       ENDIF
C
C***************** ... OTHERWISE MATRIX A IS RESTORED
C                      (I.E., A DELTA IS NOW ADDED)
C
      CALL PLSDEL(A,ID,NL,NM,NN,II1,JJ1,KK1,REAL(NL*NM*NN))
C
C***************** INVERSE FFT OF A IS NOW COMPUTED
C
       ISIG=1
       CALL C3FFT(A,ID,NL,NM,NN,WL,WM,WN,IOPT,ISIG,IORD,IWORK,IEXIT)
       IF(IEXIT.NE.0)RETURN
C
C***************** NORMALIZE AFTER TWO TRANSFORMS
C
        CALL TMSCON(A,ID,NL,NM,NN,1.0/(REAL(NL*NM*NN)))

        IF(NPRINT.NE.0) CALL SHOWRES(A,ID,NL,NM,NN)
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIX IS NOW PRINTED
C
      PRINT '(/,13X,''DIRECT PLUS INVERSE FOURIER TRANSFORM OF DATA''
     $        )'
      CALL SHOWRES(A,ID,NL,NM,NN)
      PRINT '(/////)'
C
C***************** TEST ON THE ERROR AFTER THE INVERSE TRANSFORM
C
      CALL MINMATR(A,A,B,ID,NL,NM,NN)
      ERROR=XNORM(A,ID,NL,NM,NN)
C
       IF(ERROR.GT.EPS) THEN
         IEXIT=IFFTERR
         RETURN
       ENDIF
        IEXIT=0
C
C***************** THE TEST HAS BEEN NOW SATISFACTORILY
C                  COMPLETED. EVENTUALLY PRINT HERE
C                  HOW HAPPY YOU ARE.........
C
        END
C
      INTEGER FUNCTION IPERM(I,INDEX,NELEM)
C
C***************** THIS FUNCTION CALCULATES THE POSITION OF
C                  INDEXES WHEN ORDERING IS NOT REQUESTED
C
      INTEGER INDEX(-14:*)
      IPERM=INDEX(3*NELEM+I)
      END
C
       SUBROUTINE FILL (A,B,ID,NL,NM,NN,II,JJ,KK)
C
C***************** FILL INPUT MATRIX WITH MONOCHROMATIC
C                  SIGNAL OF FREQUENCY (II,JJ,KK)
C
       COMPLEX A(0:ID-1,0:NM-1,0:NN-1),B(0:ID-1,0:NM-1,0:NN-1),CI,CJ,CK
C
       PI2=ATAN(1.)*8.
       PI2L=PI2/NL
       PI2M=PI2/NM
       PI2N=PI2/NN
C
          DO 10 I=0,NL-1
            CI=CMPLX(COS(PI2L*I*II),SIN(PI2L*I*II))
              DO 11 J=0,NM-1
                CJ=CMPLX(COS(PI2M*J*JJ),SIN(PI2M*J*JJ))
                  DO 12   K=0,NN-1
                   CK= CMPLX(COS(PI2N*K*KK),SIN(PI2N*K*KK))
C
                   A(I,J,K)=CI*CJ*CK
                   B(I,J,K)=A(I,J,K)
C
  12               CONTINUE
  11          CONTINUE
  10       CONTINUE
C
      END
        SUBROUTINE SHOWRES(A,ID,NL,NM,NN)
C
C***************** PRINT THE MARIX ELEMENTS OF A
C
        COMPLEX A(ID,NM,NN)
C
        DO 100 II=1,NN
        WRITE(6,'(/)')
        DO 100 JJ=1,NL
 100    WRITE(6,'(16F7.3)') (A(JJ,K,II),K=1,NM)
        RETURN
        END
        SUBROUTINE TMSCON(A,ID,NL,NM,NN,X)
C
C***************** MULTIPLY MATRIX A BY THE REAL NUMBER X
C
        COMPLEX A(ID,NM,NN)
        REAL X

        DO 1 I=1,NN
        DO 1 J=1,NM
        DO 1 K=1,ID
          A(K,J,I)=A(K,J,I)*X
1       CONTINUE
        END
C
        FUNCTION XNORM(A,ID,NL,NM,NN)
C
C***************** CALCULATE MAXIMUM NORM OF MATRIX A
C
        COMPLEX A(ID,NM,NN)
        XNORM=0
        DO 1 N=1,NN
        DO 1 M=1,NM
        DO 1 L=1,NL
         XNORM=MAX(XNORM,ABS(A(L,M,N)))
1       CONTINUE
        END
C
        SUBROUTINE MINDEL(A,ID,NL,NM,NN,II,JJ,KK,X)
C
C***************** SUBTRACT THE REAL NUMBER X BY MATRIX A
C
        COMPLEX A(0:ID-1,0:NM-1,0:NN-1)
C
        A(II,JJ,KK)=A(II,JJ,KK)-X
C
        END
C
        SUBROUTINE MINMATR(A,B,C,ID,NL,NM,NN)
C
C***************** MATRIX C IS SUBTRACTED BY B TO GIVE A
C
        COMPLEX A(ID,NM,NN),B(ID,NM,NN),C(ID,NM,NN)
        DO 1 N=1,NN
        DO 1 M=1,NM
        DO 1 L=1,NL
           A(L,M,N)=B(L,M,N)-C(L,M,N)
1       CONTINUE
        END
C
        SUBROUTINE PLSDEL(A,ID,NL,NM,NN,II,JJ,KK,X)
C
C***************** THE REAL NUMBER X IS ADDED TO MATRIX A
C
        COMPLEX A(0:ID-1,0:NM-1,0:NN-1)
C
        A(II,JJ,KK)=A(II,JJ,KK)+X
C
        END
ID       STATION
JOB,JN=TK3FT,MFL,T=10.
ACCOUNT,AC=TK3TSZ13,US=TK3TSZ13,APW=TTT5NHPR,UPW=TTT5NHPR.
ACCESS,DN=NOTE,OWN=TK3TSZ13.
NOTE,TEXT='USR=(vito,cray)'.
CFT,OPT=FASTMD:BTREG,ON=FZ,OFF=TC.
ACCESS,DN=FFTE.
LDR,LIB=FFTE.
EXIT.
DUMPJOB..
DEBUG,BLOCKS.
/EOF
         PROGRAM TEST
C
C***************** TESTING PROGRAM OF THE REAL 2-DIM FFT ROUTINE
C                                     R2FFT
      INTEGER NPRNT,YPRNT
      PARAMETER(NPRNT=0,YPRNT=1)
      PARAMETER(MAXTEST=1)
      PARAMETER(NSIZER=520)
C
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      INTEGER SIZERR,DFFTER,IFFTER
      PARAMETER (SIZERR=10,DFFTER=100,IFFTER=200)
C
      PRINT '(/////)'
      PRINT    '(5X,''TEST OF REAL TWO DIMENSIONAL MIXED-RADIX FFT ROUTI
     $NES'',///)'
      PRINT'(T4,''INPUT PARAMETERS ARE :'',//,T6,''ID'',T14,
     $       ''NL'',T23,''NM'',T31,''II'',T39,''JJ'',T46,
     $       ''IORD'',/)'
C
         DO 1 NTEST=1,MAXTEST
C
C***************** GENERATE MATRIX DIMENSIONS NL,NM
C                  AND ALL THE OTHER INPUT PARAMETERS
C
           CALL GENERA(ID,NL,NM,II,JJ,IORD,NSIZER)
C
           PRINT'(6(I7,1X))',ID,NL,NM,II,JJ,IORD
C
C***************** CALL TESTING ROUTINE
C
           CALL TESTFT(ID,NL,NM,II,JJ,IORD,IEXIT,NPRNT)
C
           IF(IEXIT.NE.0)THEN
             IF(IEXIT.EQ.IFFTER.OR.IEXIT.EQ.DFFTER.OR.
     $          IEXIT.EQ.TBERR                           )THEN
               PRINT '(//,9X,I7)',IEXIT
C
C***************** IN CASE OF ERROR, PRINT THE CONTENT OF DATA MATRIX
C                  ONLY IF SMALL ENOUGH (I.E. MAX(NL,NM) LESS THAN 32)
C
               IF(MAX(NL,NM).LE.32)THEN
                 CALL TESTFT(ID,NL,NM,II,JJ,IORD,IEXIT,YPRNT)
               ENDIF
               STOP
             ENDIF
           ENDIF
C
1        CONTINUE
C
C
         END
C
C
      SUBROUTINE GENERA(ID,NL,NM,II,JJ,IORD,NSIZE)
C
C***************** THIS ROUTINE PERFORMS A RANDOM GENERATION OF THE
C                  INPUT DATA DIMENSIONS NL,NM, AS WELL AS OF ALL
C                  THE CALLING PARAMETERS OF THE REAL 2-DIM FFT
C                  ROUTINE R2FFT
C
C
      REAL MAXDIM
      PARAMETER (MAXDIM=850)
C
       L2MAX=(LOG(REAL(NSIZE*2)))/LOG(2.0)
       L3MAX=(LOG(REAL(NSIZE*3)))/LOG(3.0)
       L5MAX=(LOG(REAL(NSIZE*5)))/LOG(5.0)
C
10    CONTINUE
      L2=RANF()*L2MAX
      L3=RANF()*L3MAX
      L5=RANF()*L5MAX
C
      M2=RANF()*L2MAX
      M3=RANF()*L3MAX
      M5=RANF()*L5MAX
C
      NL=2**L2*3**L3*5**L5
      IF(NL.EQ.1) NL=2
      NM=2**M2*3**M3*5**M5
      IF(NM.EQ.1) NM=2
C
C***************** GENERATED DIMENSIONS SHOULD NOT EXCEED
C                  MAXIMUM MEMORY STORAGE
C
      IF(NL*NM.GT.NSIZE*NSIZE) GO TO 10
      IF(NL.GT.MAXDIM.OR.NM.GT.MAXDIM) GO TO 10
C
      ID=NL+2
C
C***************** GENERATE RANDOM POSITIONS OF INPUT FREQUENCIES
C
      II=RANF()*NL
      JJ=RANF()*NM
C
      IORD=RANF()*2
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT PARAMETERS ARE FORCED TO
C                  ASSUME THE FOLLOWING VALUES..........
C
      NL=10
      NM=4
      ID=12
      II=1
      JJ=1
      IORD=1
C
      RETURN
      END
C
C
      SUBROUTINE TESTFT(ID,NL,NM,III,JJJ,IORD,IEXIT,NPRINT)
C
C***************** THIS ROUTINE PERFORMS THE TESTS ON R2FFT
C
      INTEGER SIZERR, DFFTER,IFFTER
      PARAMETER (SIZERR=10,DFFTER=100,IFFTER=200)
C
      PARAMETER (NSIZER=520,NSIZEC=260)
      COMPLEX AC(NSIZEC*NSIZEC),B(NSIZEC*NSIZEC),WEIGHT
      REAL AR(NSIZER*NSIZEC)
C
      INTEGER WL(6*NSIZEC+14),WM(4*NSIZEC+14)
      INTEGER IWORK(NSIZER)
C
      EQUIVALENCE (AR,AC)
C
      DATA EPS/1.E-7/
C
C***************** SET IEXIT TO 1 (MFFTP SHOULD RESET IT TO 0)
C
      IEXIT=1
C
      NL1=NL/2
      NPROD2=NL1*NM
C
C***************** FILL INPUT MATRICES WITH
C                  MONOCHROMATIC SIGNALS
C
      DO 10 I=1,NSIZER*NSIZEC
10    AR(I)=0.
      DO 11 J=1,NSIZEC*NSIZEC
 11   B(J)=0.
      CALL FILL2R(AR,B,ID,NL,NL1,NM,III,JJJ,WEIGHT)
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT MATRIX IS NOW PRINTED
C
      PRINT '(//,18X,''INPUT DATA MATRIX'',/)'
      CALL RSHOW2(AR,ID,NL,NM)
C
C***************** INITIALIZE FFT PACKAGE
C
      CALL R2FFT(AR,ID,NL,NM,WL,WM,0,IORD,IWORK,IEXIT)
      IF(IEXIT.NE.O) RETURN
C
C***************** CALCULATE DIRECT FFT
C
      ISIG=-1
      CALL R2FFT(AR,ID,NL,NM,WL,WM,ISIG,IORD,IWORK,IEXIT)
      IF(IEXIT.NE.0) RETURN
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIIX IS NOW
C                  PRINTED
C
      PRINT '(/,13X,''DIRECT FOURIER TRANSFORM OF INPUT MATRIX'',/)'
      CALL SHOWRE2(AC,ID/2,NL,NM)
C
C***************** COMPARE THE RESULT WITH KNOWN EXACT VALUE
C                  A 2-D DELTA IS SUBTRACTED FROM THE RESULT
C                  THE POSITION OF DELTA IS (IPOSDE,JPOSDE)
C
      II=ABS(III)
      JJ=ABS(JJJ)
      IPOSDE=II
      JPOSDE=JJ
      IF(II.GT.NL1)THEN
        IPOSDE=NL-II
        IF(JJ.NE.0)JPOSDE=NM-JJ
        WEIGHT=CONJG(WEIGHT)
      ENDIF
      IF(IORD.NE.0) THEN
C
      CALL MINDE2(AC,ID/2,NL1,NM,IPOSDE,JPOSDE,NL1*NM*WEIGHT)
      IF(II.EQ.0.OR.II.EQ.NL1) THEN
        IF(JPOSDE.NE.0)THEN
        CALL MINDE2(AC,ID/2,NL1,NM,IPOSDE,NM-JPOSDE,NL1*NM*CONJG(WEIGHT)
     $)
        ELSE
         CALL MINDE2(AC,ID/2,NL1,NM,IPOSDE,JPOSDE,NL1*NM*CONJG(WEIGHT))
        ENDIF
      ENDIF
C
C***************** THE ERROR IS MAX(ABS(AC))
C
      ERROR=XNOR2(AC,ID/2,NL1+1,NM)/REAL(NL1*NM)
C
C***************** IF THE ERROR IS LARGE, THEN AN ERROR
C                  MESSAGE IS ISSUED, AND THE TEST IS STOPPED
C
      IF(ERROR.GT.EPS) THEN
      IEXIT=DFFTER
      RETURN
      ENDIF
C
C***************** ... OTHERWISE MATRIX AC IS RESTORED
C                      (I.E., A DELTA IS NOW ADDED)
C
      CALL PLSDE2(AC,ID/2,NL1,NM,IPOSDE,JPOSDE,NL1*NM*WEIGHT)
      IF(II.EQ.0.OR.II.EQ.NL1) THEN
        IF(JPOSDE.NE.0)THEN
        CALL PLSDE2(AC,ID/2,NL1,NM,IPOSDE,NM-JPOSDE,NL1*NM*CONJG(WEIGHT)
     $)
       ELSE
         CALL PLSDE2(AC,ID/2,NL1,NM,IPOSDE,JPOSDE,NL1*NM*CONJG(WEIGHT))
       ENDIF
      ENDIF
C
      ENDIF
C
C***************** INVERSE FFT IS COMPUTED NOW
C
      ISIG=1
      CALL R2FFT(AC,ID,NL,NM,WL,WM,ISIG,IORD,IWORK,IEXIT)
      IF(IEXIT.NE.0) RETURN
C
C***************** NORMALIZE AFTER TWO TRANSFORMS
C
      CALL TMSCO2(AC,ID/2,NL1,NM,1./(NL*NM))
      IF(NPRINT.NE.0) CALL  RSHOW2(AR,ID,NL,NM)
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIX IS NOW
C                  PRINTED
C
       PRINT '(/,4X,''DIRECT PLUS INVERSE FOURIER TRANSFORM OF DATA''
     $,/)'
      CALL RSHOW2(AR,ID,NL,NM)
      PRINT '(/////)'
C
C***************** CALCULATE THE ERROR ON THE INVERSE TRANSFORM
C
      CALL MINMAT2(AC,AC,B,ID/2,NL1,NM)
      ERROR=XNOR2(AC,ID/2,NL1,NM)/REAL(NL1*NM)
C
C***************** TEST ON THE ERROR AFTER THE INVERSE TRANSFORM
C
      IF(ERROR.GT.EPS) THEN
      IEXIT=IFFTER
      RETURN
C
C
      ENDIF
C
C***************** THE TEST HAS NOW BEEN SATISFACTORILY
C                  COMPLETED. EVENTUALLY PRINT HERE HOW
C                  HAPPY YOU ARE.....
C

      IEXIT=0
      END
C
C
      SUBROUTINE FILL2R(AR,B,ID,NL,NL1,NM,II,JJ,WEIGHT)
C
C***************** FILLS REAL INPUT MATRIX WITH MONOCHROMATIC
C                  SIGNAL, OF FREQUENCY (II,JJ)
C
      REAL AR(0:ID-1,0:*),B(0:ID-1,0:*)
      COMPLEX WEIGHT
C
      PI2=8*ATAN(1.)
C
      IF(II.GE.0.AND.JJ.GE.0)THEN
        WEIGHT=1.
        DO 10 I=0,NL-1
        DO 10 J=0,NM-1
          AR(I,J)=COS(PI2*(I*II/REAL(NL)+J*JJ/REAL(NM)))
          B(I,J)=AR(I,J)
10      CONTINUE
      ELSE
        WEIGHT=(0.,-1.)
        DO 20 I=0,NL-1
        DO 20 J=0,NM-1
          AR(I,J)=SIN(PI2*(I*ABS(II)/REAL(NL)+J*ABS(JJ)/REAL(NM)))
          B(I,J)=AR(I,J)
20      CONTINUE
      ENDIF
C
      END
      SUBROUTINE RSHOW2(A,ID,NL,NM)
C
C***************** PRINTS THE CONTENT OF REAL MATRIX A
C
      REAL A(ID,NM)
C
      DO 100 JJ=1,NL
 100  WRITE(6,'(10X,16F8.5)') (A(JJ,K),K=1,NM)
C
      END
      SUBROUTINE SHOWRE2(A,ID,NL,NM)
C
C***************** PRINT MATRIX CONTENT
C
      COMPLEX A(ID,NM)
C
      DO 100 JJ=1,ID
 100  WRITE(6,'(16F8.5)') (A(JJ,K),K=1,NM)
      END
C
      SUBROUTINE TMSCO2(A,ID,NL,NM,X)
C
C***************** MULTIPLY MATRIX ELEMENTS BY X
C
      COMPLEX A(ID,NM)
      REAL X
      DO 1 I=1,NM
      DO 1 J=1,ID
      A(J,I)=A(J,I)*X
 1    CONTINUE
      END
C
      FUNCTION XNOR2(A,ID,NL,NM)
C
C***************** CALCULATE MAXIMUM NORM OF MATRIX A
C
      COMPLEX A(ID,NM)
      XNOR2=0.
      DO 1 L=1,NL
      DO 1 M=1,NM
      XNOR2=MAX(XNOR2,ABS(REAL(A(L,M)))+ABS(AIMAG(A(L,M))))
 1    CONTINUE
      END
C
      SUBROUTINE MINDE2(A,ID,NL,NM,II,JJ,X)
C
C***************** NUMBER X IS SUBTRACTED FROM A(II,JJ)
C
      COMPLEX A(0:ID-1,0:NM-1)
      A(II,JJ)=A(II,JJ)-X
      END
      SUBROUTINE MINMAT2(A,B,C,ID,NL,NM)
C
C***************** MATRIX C IS SUBTRACTED FROM MATRIX B
C                  AND STORED IN A
C
      COMPLEX A(ID,NM),B(ID,NM),C(ID,NM)
      DO 1 L=1,NL
      DO 1 M=1,NM
      A(L,M)=B(L,M)-C(L,M)
 1    CONTINUE
      END
      SUBROUTINE PLSDE2(A,ID,NL,NM,II,JJ,X)
C
C***************** NUMBER X IS ADDED TO A(II,JJ)
C
      COMPLEX A(0:ID-1,0:NM-1)
      A(II,JJ)=A(II,JJ)+X
      END
ID       STATION
JOB,JN=TK3FT,MFL,T=10.
ACCOUNT,AC=TK3TSZ13,US=TK3TSZ13,APW=TTT5NHPR,UPW=TTT5NHPR.
ACCESS,DN=NOTE,OWN=TK3TSZ13.
NOTE,TEXT='USR=(vito,cray)'.
CFT,OPT=FASTMD:BTREG,ON=FZ,OFF=TC.
ACCESS,DN=FFTE.
LDR,LIB=FFTE.
EXIT.
DUMPJOB..
DEBUG,BLOCKS.
/EOF
         PROGRAM TEST
C
C***************** TESTING PROGRAM OF THE REAL 3-DIM FFT ROUTINE
C                                       R3FFT
C
      INTEGER NPRNT,YPRNT
      PARAMETER(NPRNT=0,YPRNT=1)
      PARAMETER(MAXTEST=1)
      PARAMETER(NSIZE=65)
C
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      INTEGER SIZERR,DFFTER,IFFTER
      PARAMETER (SIZERR=10,DFFTER=100,IFFTER=200)
C
      PRINT '(/////)'
      PRINT '(6X,''TEST OF REAL THREE DIMENSIONAL MIXED-RADIX FFT ROUTIN
     $ES'',//)'
      PRINT '(1X,''INPUT PARAMETERS ARE :'',/,
     $T6,''ID'',T14,''NL'',T22,''NM'',T30,''NN'',T38,''II''
     $,T46,''JJ'',T54,''KK'',T62,''IOPT'',T68,''IORD'')'
C
         DO 1 NTEST=1,MAXTEST
C
C***************** GENERATE MATRIX DIMENSIONS NL,NM,NN
C                  AND ALL THE OTHER INPUT PARAMETERS
C
           CALL GENERA(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,NSIZE)
           PRINT'(/9(I7,1X))',ID,NL,NM,NN,II,JJ,KK,IOPT,IORD
C
C
C***************** CALL TESTING ROUTINE
C
           CALL TESTFT(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,IEXIT,NPRNT)
C
           IF(IEXIT.NE.0)THEN
             IF(IEXIT.EQ.IFFTER.OR.IEXIT.EQ.DFFTER.OR.
     $          IEXIT.EQ.TBERR                           )THEN
               PRINT '(//,9X,I7)',IEXIT
C
C***************** IN CASE OF ERROR, PRINT THE CONTENT OF DATA MATRIX
C                  ONLY IF SMALL ENOUGH (I.E. MAX(NL,NM) LESS THAN 32)
C
               IF(MAX(NL,NM).LE.32)THEN
                 CALL TESTFT(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,IEXIT,YPRNT)
               ENDIF
               STOP
             ENDIF
           ENDIF
C
1        CONTINUE
C
         END
C
      SUBROUTINE GENERA(ID,NL,NM,NN,II,JJ,KK,IOPT,IORD,NSIZE)
C
C***************** THIS ROUTINE PERFORMS A RANDOM GENERATION OF THE
C                  INPUT DATA DIMENSIONS NL,NM,NN AS WELL AS OF ALL
C                  THE CALLING PARAMETERS OF THE REAL 3-DIM FFT
C                  ROUTINE R3FFT
C
      REAL MAXDIM
      PARAMETER (MAXDIM=512)
C
       L2MAX=(LOG(REAL(NSIZE*2)))/LOG(2.0)
       L3MAX=(LOG(REAL(NSIZE*3)))/LOG(3.0)
       L5MAX=(LOG(REAL(NSIZE*5)))/LOG(5.0)
C
10    CONTINUE
C
      L2=RANF()*L2MAX
      L3=RANF()*L3MAX
      L5=RANF()*L5MAX
C
      M2=RANF()*L2MAX
      M3=RANF()*L3MAX
      M5=RANF()*L5MAX
C
      N2=RANF()*L2MAX
      N3=RANF()*L3MAX
      N5=RANF()*L5MAX
C
      NL=2**L2*3**L3*5**L5
      IF(NL.EQ.1) NL=2
      NM=2**M2*3**M3*5**M5
      IF(NM.EQ.1) NM=2
      NN=2**N2*3**N3*5**N5
      IF(NN.EQ.1) NN=2
C
      ID=NL+2
C
C***************** GENERATED DIMENSIONS SHOULD NOT EXCEED
C                  MAXIMUM STORAGE CAPACITY
C
      IF(MOD(NL,2).NE.0) GO TO 10
      IF(ID*NM*NN.GT.NSIZE*NSIZE*NSIZE) GO TO 10
      IF(NL.GT.MAXDIM.OR.NM.GT.MAXDIM.OR.NN.GT.MAXDIM) GO TO 10
C
      II=RANF()*NL
      JJ=RANF()*NM
      KK=RANF()*NN
C
      IORD=RANF()*2
      IOPT=RANF()*2
C
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT PARAMETERS ARE FORCED TO
C                  ASSUME THE FOLLOWING VALUES..........
C
      NL=8
      NM=5
      NN=2
      ID=10
      II=1
      JJ=1
      KK=1
      IOPT=1
      IORD=1
C
      RETURN
      END
C
      SUBROUTINE TESTFT(ID,NL,NM,NN,III,JJJ,KKK,IOPT,IORD,IEXIT,NPRINT)
C
C***************** THIS ROUTINE PERFORMS THE TESTS ON R3FFT
C
      INTEGER SIZERR, DFFTER,IFFTER
      PARAMETER (SIZERR=10,DFFTER=100,IFFTER=200)
C
      PARAMETER (NSIZER=130,NSIZEC=65)
      COMPLEX AC(NSIZEC*NSIZEC*NSIZEC),B(NSIZEC*NSIZEC*NSIZEC),WEIGHT
      INTEGER WL(6*NSIZEC+14),WM(4*NSIZEC+14),WN(4*NSIZEC+14)
      INTEGER IWORK(NSIZEC)
C
      REAL AR(NSIZER*NSIZEC*NSIZEC)
C
      EQUIVALENCE (AR,AC)
C
      DATA EPS/1.E-7/
C
C
C***************** SET IEXIT TO 1 (MFFTP SHOULD RESET IT TO 0)
C
      IEXIT=1
C
      NL1=NL/2
      NPROD3=NL1*NM*NN
C
C***************** FILL INPUT MATRICES WITH
C                  MONOCHROMATIC SIGNALS
C
      DO 10 I=1,NSIZER*NSIZEC*NSIZEC
10    AR(I)=0.
      DO 11 J=1,NSIZEC*NSIZEC*NSIZEC
 11   B(J)=0.
      CALL FILL3R(AR,B,ID,NL,NM,NN,III,JJJ,KKK,WEIGHT)
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE INPUT MATRIX IS NOW PRINTED
C
      PRINT '(///,20X,''INPUT DATA MATRIX'')'
      CALL RSHOW3(AR,ID,NL,NM,NN)
C
C***************** INITIALIZE FFT PACKAGE
C
      CALL R3FFT(AR,ID,NL,NM,NN,WL,WM,WN,IOPT,0,IORD,IWORK,IEXIT)
      IF (IEXIT.NE.0) RETURN
C
C***************** CALCULATE DIRECT FOURIER TRANSFORM OF A
C
      ISIG=-1
      CALL R3FFT(AR,ID,NL,NM,NN,WL,WM,WN,IOPT,ISIG,IORD,IWORK,IEXIT)
      IF(IEXIT.NE.0) RETURN
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIX IS NOW PRINTED
C
      PRINT '(/9X,''DIRECT FOURIER TRANSFORM OF INPUT MATRIX'')'
      CALL SHOWRES(AC,ID/2,NL,NM,NN)
C
C***************** COMPARE THE RESULT WITH KNOWN EXACT VALUE
C                  A 3-D DELTA IS SUBTRACTED FROM THE RESULT
C                  THE POSITION OF DELTA IS (IPOSDE,JPOSDE,KPOSDE)
C
      II=ABS(III)
      JJ=ABS(JJJ)
      KK=ABS(KKK)
C
      IPOSDE=II
      JPOSDE=JJ
      KPOSDE=KK
C
      IF(II.GT.NL1) THEN
      IPOSDE=NL-II
      IF(JJ.NE.0) JPOSDE=NM-JJ
      IF(KK.NE.0) KPOSDE=NN-KK
      WEIGHT=CONJG(WEIGHT)
      ENDIF
C
      IF(IORD.NE.0) THEN
C
C     PRINT*,'-DELTA'
        CALL MINDEL(AC,ID/2,NL1,NM,NN,IPOSDE,JPOSDE,KPOSDE,
     $              NPROD3*WEIGHT)
      IF(II.EQ.0.OR.II.EQ.NL1) THEN
C
        IF(JPOSDE.NE.0.AND.KPOSDE.NE.0) THEN
          CALL MINDEL (AC,ID/2,NL1,NM,NN,IPOSDE,NM-JPOSDE,NN-KPOSDE,
     $                 NPROD3*CONJG(WEIGHT))
        ELSE
          IF(JPOSDE.NE.0) THEN
            CALL MINDEL (AC,ID/2,NL1,NM,NN,IPOSDE,NM-JPOSDE,KPOSDE,
     $                   NPROD3*CONJG(WEIGHT))
          ENDIF
C
          IF(KPOSDE.NE.0) THEN
            CALL MINDEL (AC,ID/2,NL1,NM,NN,IPOSDE,JPOSDE,NN-KPOSDE,
     $                   NPROD3*CONJG(WEIGHT))
          ENDIF
C
          IF(JPOSDE.EQ.0.AND.KPOSDE.EQ.0) THEN
            CALL MINDEL(AC,ID/2,NL1,NM,NN,IPOSDE,JPOSDE,KPOSDE,
     $                  NPROD3*CONJG(WEIGHT))
          ENDIF
C
        ENDIF
C
      ENDIF
C
C***************** ERROR IS MAX(ABS(AC))
C
      ERROR=XNORM(AC,ID/2,NL1+1,NM,NN)/REAL(NPROD3)
C
C***************** IF THE ERROR IS LARGE, THEN AN ERROR
C                  MESSAGE IS ISSUED, AND THE TEST IS STOPPED
C
      IF(ERROR.GT.EPS) THEN
        IEXIT=DFFTER
        RETURN
      ENDIF
C
C***************** ... OTHERWISE MATRIX AC IS RESTORED
C                      (I.E., A DELTA IS NOW ADDED)
C
      CALL PLSDEL(AC,ID/2,NL1,NM,NN,IPOSDE,JPOSDE,KPOSDE,
     $              NPROD3*WEIGHT)
      IF(II.EQ.0.OR.II.EQ.NL1) THEN
        IF(JPOSDE.NE.0.AND.KPOSDE.NE.0) THEN
          CALL PLSDEL(AC,ID/2,NL1,NM,NN,IPOSDE,NM-JPOSDE,NN-KPOSDE,
     $                NPROD3*CONJG(WEIGHT))
C
        ELSE
C
          IF(JPOSDE.NE.0) THEN
            CALL PLSDEL(AC,ID/2,NL1,NM,NN,IPOSDE,NM-JPOSDE,KPOSDE,
     $                  NPROD3*CONJG(WEIGHT))
          ENDIF
C
          IF(KPOSDE.NE.0) THEN
            CALL PLSDEL(AC,ID/2,NL1,NM,NN,IPOSDE,JPOSDE,NN-KPOSDE,
     $                  NPROD3*CONJG(WEIGHT))
          ENDIF
C
          IF(JPOSDE.EQ.0.AND.KPOSDE.EQ.0) THEN
            CALL PLSDEL(AC,ID/2,NL1,NM,NN,IPOSDE,JPOSDE,KPOSDE,
     $                  NPROD3*CONJG(WEIGHT))
          ENDIF
C
      ENDIF
C
      ENDIF
C
      ENDIF
C
C***************** INVERSE FFT IS NOW COMPUTED
C
      ISIG=1
      CALL R3FFT(AC,ID,NL,NM,NN,WL,WM,WN,IOPT,ISIG,IORD,IWORK,IEXIT)
      IF(IEXIT.NE.0) RETURN
C
C***************** NORMALIZE AFTER TWO TRANSFORMS
C
      CALL TMSCON(AC,ID/2,NL1,NM,NN,1./(NL*NM*NN))
C
C***************** THE PRESENT VERSION IS SUITED TO PUBLICATION
C                  THEREFORE THE TRANSFORMED MATRIX IS NOW PRINTED
C
      PRINT '(/,6X,''DIRECT PLUS INVERSE FOURIER TRANSFORM OF DATA''
     $       )'
      CALL RSHOW3(AR,ID,NL,NM,NN)
      PRINT '(/////)'
C
C***************** TEST ON THE ERROR AFTER THE INVERSE TRANSFORM
C
      CALL MINMATR(AC,AC,B,ID/2,NL1,NM,NN)
      ERROR=XNORM(AC,ID/2,NL1,NM,NN)/REAL(NL*NM*NN)
C
      IF(ERROR.GT.EPS) THEN
        IEXIT=IFFTER
        RETURN
      ENDIF
C
C***************** THE TEST HAS NOW BEEN SATISFACTORILY
C                  COMPLETED. EVENTUALLY PRINT HERE
C                  HOW HAPPY YOU ARE...................
C
C
      IEXIT=0
      END
C
C
      SUBROUTINE FILL3R(AR,B,ID,NL,NM,NN,II,JJ,KK,WEIGHT)
C
C****************** FILL INPUT MATRIX WITH MONOCHROMATIC
C                   SIGNAL, OF FREQUENCY (II,JJ,KK)
C
      REAL AR(0:ID-1,0:NM-1,0:NN-1),B(0:ID-1,0:NM-1,0:NN-1)
      COMPLEX WEIGHT
C
      PI2=8*ATAN(1.)
C
      IF(II.GE.0.AND.JJ.GE.0.AND.KK.GE.0) THEN
      WEIGHT=1.
      DO 10 I=0,NL-1
      DO 10 J=0,NM-1
      DO 10 K=0,NN-1
C
      AR(I,J,K)=COS(PI2*(I*II/REAL(NL)+J*JJ/REAL(NM)+K*KK/REAL(NN)))
      B(I,J,K)=AR(I,J,K)
C
 10   CONTINUE
C
        ELSE
C
      WEIGHT=(0.,-1.)
      DO 20 I=0,NL-1
      DO 20 J=0,NM-1
      DO 20 K=0,NN-1
C
      AR(I,J,K)=SIN(PI2*(I*ABS(II)/REAL(NL)+J*ABS(JJ)/REAL(NM)+
     $                   K*ABS(KK)/REAL(NN)))
      B(I,J,K)=AR(I,J,K)
 20   CONTINUE
C
      ENDIF
C
      END
C
      SUBROUTINE RSHOW3(A,ID,NL,NM,NN)
C
C***************** PRINT THE CONTENT OF REAL MATRIX A
C
      REAL A(ID,NM,NN)
C
      DO 100 KK=1,NN
C
      WRITE(6,'(/)')
C
      DO 100 II=1,NL
 100  WRITE (6,'(8X,16F8.5)') (A(II,M,KK),M=1,NM)
C
C
      END
        SUBROUTINE SHOWRES(A,ID,NL,NM,NN)
C
C***************** PRINT THE MARIX ELEMENTS OF A
C
        COMPLEX A(ID,NM,NN)
C
        DO 100 II=1,NN
        WRITE(6,'(/)')
        DO 100 JJ=1,ID
 100    WRITE(6,'(16F7.3)') (A(JJ,K,II),K=1,NM)
        RETURN
        END
        SUBROUTINE TMSCON(A,ID,NL,NM,NN,X)
C
C***************** MULTIPLY MATRIX A BY THE REAL NUMBER X
C
        COMPLEX A(ID,NM,NN)
        REAL X
        DO 1 I=1,NN
        DO 1 J=1,NM
        DO 1 K=1,ID
          A(K,J,I)=A(K,J,I)*X
1       CONTINUE
        END
C
        FUNCTION XNORM(A,ID,NL,NM,NN)
C
C***************** CALCULATE MAXIMUM NORM OF MATRIX A
C
        COMPLEX A(ID,NM,NN)
        XNORM=0
        DO 1 N=1,NN
        DO 1 M=1,NM
        DO 1 L=1,NL
         XNORM=MAX(XNORM,ABS(A(L,M,N)))
1       CONTINUE
        END
C
        SUBROUTINE MINDEL(A,ID,NL,NM,NN,II,JJ,KK,X)
C
C***************** SUBTRACT THE REAL NUMBER X BY MATRIX A
C
        COMPLEX A(0:ID-1,0:NM-1,0:NN-1)
C
        A(II,JJ,KK)=A(II,JJ,KK)-X
C
        END
C
        SUBROUTINE MINMATR(A,B,C,ID,NL,NM,NN)
C
C***************** MATRIX C IS SUBTRACTED BY B TO GIVE A
C
        COMPLEX A(ID,NM,NN),B(ID,NM,NN),C(ID,NM,NN)
        DO 1 N=1,NN
        DO 1 M=1,NM
        DO 1 L=1,NL
           A(L,M,N)=B(L,M,N)-C(L,M,N)
1       CONTINUE
        END
C
        SUBROUTINE PLSDEL(A,ID,NL,NM,NN,II,JJ,KK,X)
C
C***************** THE REAL NUMBER X IS ADDED TO MATRIX A
C
        COMPLEX A(0:ID-1,0:NM-1,0:NN-1)
C
        A(II,JJ,KK)=A(II,JJ,KK)+X
C
        END

.