C
C     THIS SET OF PROGRAM AND SUBROUTINE UNITS IS TO SUPPORT
C
C       'AN OPTIMIZED MASS STORAGE FFT',  BY DONALD FRASER
C
C           REVISION DATE: JULY 1978 (MINOR REV JUNE 79).
C
C     THE SET INCLUDES A UNIVERSAL TEST DRIVER PROGRAM, WHICH SIMULATES
C  MASS STORE THROUGH A FORTRAN ARRAY, SAMPLE PROGRAMS AND I/O
C  SUBROUTINES FOR CONTROL DATA 6000 AND CYBER COMPUTERS AND FOR
C  DEC PDP 11 MINICOMPUTERS.
C
C     I/O SUBROUTINES FOR OTHER SYSTEMS SHOULD BE EASILY CONSTRUCTED
C  FROM THE EXAMPLES AND WITH REFERENCE TO THE FORMAL PAPER.
C  BUT CARE SHOULD BE TAKEN WITH SOME FUSSY COMPILERS SINCE FFT
C  SUBROUTINES ASSIGN EITHER TYPE REAL OR TYPE COMPLEX TO THE SAME
C  ARRAY AS IT IS PASSED AS A FORMAL PARAMETER.  NOTE THAT COMPLEX
C  DATA IS ASSUMED TO BE STORED  REAL/IMAG/REAL/IMAG... IN BOTH
C  MASS STORE AND CORE STORE.
C
C     THE PROGRAM UNITS APPEAR IN THE FOLLOWING ORDER:
C
C         FIRST, THE FFT SUBROUTINE SET:
C  1  RMFFT        OPTIMIZED MASS STORAGE FFT (REAL DATA OR REAL RESULT)
C  2  CMFFT        CALLED BY 1, OR MASS STORAGE FFT (ALL COMPLEX DATA)
C  3  MFCOMP       IN-CORE FFT
C  4  MFSORT       MASS STORE SORTING
C  5  MFREV        IN-CORE SORTING OR WHOLE BLOCK SORTING
C  6  MFLOAD       LOADING/UNLOADING CORE STORE
C  7  MFINDX       BLOCK INDEXING ALGORITHM (VIRTUAL PERMUTATION)
C  8  MFSUM        MEXA() EXPONENT SUMMATIONS
C  9  MFRCMP       REAL-COMPLEX UNSCRAMBLING/SCRAMBLING
C 10  MFRLOD       LOADING/UNLOADING CORE STORE FOR MFRCMP
C 11  MFPAR        HELPER ROUTINE (NOT ESSENTIAL, BUT RECOMMENDED)
C 12  DMPERM       MASS STORE DIMENSION SHIFTING (BONUS SUBROUTINE)
C
C           THEN, TEST PROGRAMS AND SAMPLE I/O SUBROUTINES
C
C 13  UNIVERSAL TEST PROGRAM (SIMULATED MASS STORE, NEEDS 14 TO 17 ALSO)
C 14  RANMF        RANDOM NUMBER GENERATOR
C 15  NAIVE        DISCRETE FOURIER TRANSFORM COMPUTED NAIVELY
C 16  MFREAD       DUMMY I/O ROUTINE USING SIMULATED MASS STORE
C 17  MFWRIT          AS ABOVE
C
C 18  CYBER MASS STORE SAMPLE PROGRAM
C 19  MFREAD/MFWRIT  CYBER MASS STORE I/O ROUTINE
C
C 20  CYBER EXTENDED CORE/LCM SAMPLE PROGRAM
C 21  MFREAD/MFWRIT  CYBER EXTENDED CORE/LCM I/O ROUTINES
C
C 22  PDP 11 MASS STORE SAMPLE PROGRAM
C 23  MFREAD       PDP 11 STANDARD FORTRAN MASS STORE I/O ROUTINES
C 24  MFWRIT          AS ABOVE
C
C 25  PDP 11 FAST MACRO I/O SAMPLE PROGRAM (NEEDS MACRO OPEN SUBRTN)
C 26  MFREAD/MFWRIT  PDP 11 FAST MACRO I/O ROUTINE FOR RSX11M/RT11.
C
      END
      SUBROUTINE RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C
C REAL-TO-COMPLEX FFT (OR VICE-VERSA) OF MULTI-DIMENSD MASS STORE ARRAY
C  (FRASER, ACM TOMS - 1978/79, AND J.ACM, V.23,N.2, APRIL 76, PP. 298-3
C  MASS STORE ARRAY IS EITHER REAL OR COMPLEX DATA (SEE NOTE BELOW)
C
C  NOTE WELL THAT TYPE COMPLEX DATA MUST EXIST AS ALTERNATING
C  REAL/IMAG/REAL/IMAG... ELEMENTS, BOTH IN MASS STORE AND IN FORTRAN
C  WORKING ARRAY BUFA;  IN THIS FFT, DIFFERENT SUBROUTINES WILL SET
C  DIFFERENT TYPE (REAL OR COMPLEX) FOR ARRAY BUFA.
C
C MEXA(J) LIST OF DIMENSION SIZE EXPONS (BASE 2), ADJACENT VARIABLES FIR
C NDIM IS NUMBER OF EXPONENTS IN LIST AND THUS THE NUMBER OF DIMENSIONS
C  SUM TO NDIM OF MEXA(J) = M, WHERE 2**M IS EFFECT SIZE OF MASS STORE A
C     THUS, 2**M PACKED REAL VALUES,
C     OR,   2**M COMPLEX VALUES, IF COMPLEX RESULT WITH IPAK=1 (SEE BELO
C RMFFT ALWAYS REVERSES DIMENSION ORDER AND MEXA LIST (TRANSPOSED, IF 2
C
C ISGN GIVES SIGN OF COMPLEX EXPONENT OF TRANSFORM (+ OR -), AND
C IDIR DETERMINES DIRECTION OF TRANSFORM, THUS:
C     IDIR=-1,  REAL-TO-COMPLEX
C     IDIR=+1,  COMPLEX-TO-REAL
C SCAL IS REAL MULTIPLIER OF RESULT (EG. SET SCAL=1. FWD, 1./2**M INV)
C
C BUFA IS CORE STORE WORKING ARRAY BASE ADDRESS (SEE NOTE ABOVE)
C IBEX, ICEX ARE BLOCK AND CORE SIZE EXPONENTS, THUS
C  2**IBEX IS NUMBER OF REAL ELEMENTS IN MASS STORE BLOCK
C  2**ICEX IS NUMBER OF REAL ELEMENTS IN CORE STORE BUFA
C
C IPAK IS ARRAY PACKING DETERMINATOR, THUS:
C  IPAK=+1 EXPANDS COMPLEX ARRAY TO FULL REDUNDANCY (SAME AS SUBRTN CMFF
C  IPAK=0  COMPUTES COMPLEX ARRAY OF JUST OVER HALF SIZE (COMMON METHOD)
C  IPAK=-1 HOLDS COMPLEX ARRAY AT EXACTLY HALF SIZE (2**(M-1) CMPLX ELMT
C
C MASS STORE ARRAY MUST BE OPEN FOR ACCESS BY USER SUBRTNS MFREAD/MFWRIT
C  EG. SUBRTN MFREAD(BUFA,NB,JB) AND MFWRIT(BUFA,NB,JB) TRANSFER ONE
C  BLOCK, INDEX JB, BETWEEN MASS STORE AND CORE STORE BUFA, NB REALS
C  (1.LE.JB.LE.2**(M-IBEX) IF REAL, OR 2**(M-IBEX+1) IF CMPLX AND IPAK=1
C
      COMPLEX BUFA(1)
      INTEGER B,MEXA(1)
C
      MH=MFSUM(MEXA,NDIM,99)-1
      B=IBEX-1
      IF(IDIR.GT.0)GO TO 10
C
C BELOW, REAL-TO-COMPLEX TRANSFORM
      MEXA(1)=MEXA(1)-1
      CALL CMFFT  (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
      CALL MFRCMP (MEXA,NDIM,ISGN,IDIR,IPAK, BUFA,B,MH)
      MEXA(NDIM)=MEXA(NDIM)+1
      RETURN
C
C BELOW, COMPLEX-TO-REAL TRANSFORM
10    MEXA(NDIM)=MEXA(NDIM)-1
      CALL MFRCMP (MEXA,NDIM,ISGN,IDIR,IPAK, BUFA,B,MH)
      CALL CMFFT  (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
      MEXA(1)=MEXA(1)+1
      RETURN
C
      END
      SUBROUTINE CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C
C COMPLEX FFT OF MULTI-DIMENSD MASS STORE ARRAY, CALLED BY USER OR RMFFT
C  FOR COMMENTS, SEE SUBRTN RMFFT; ARGUMNTS HAVE SAME MEANING EXCEPT FOR
C  IDIR=+1 OR -1, ARRAY ALWAYS COMPLEX, DIMENSION ORDER REVERSED
C  WHILE IDIR=0, DIMENSION ORDER (AND MEXA LIST) ARE NOT REVERSED
C
      COMPLEX BUFA(1)
      INTEGER B,C,MEXA(1)
      DATA LSET/1/,LREAD/1/,LWRIT/2/
C
      M=MFSUM(MEXA,NDIM,99)
      MREAL=M+1
      B=IBEX-1
      C=ICEX-1
      NC=2**C
      IPAS=0
      MPAS=M-C
      NPAS=C-B
C MOST EFFICIENT USE OF CORE STORE - TRIES TO DO C-B PASSES PER LOAD
      IDUM=MFINDX(LSET,B,M,M,NPAS)
C DUMMY CALL TO MFINDX TO SPECIFY VIRTUAL (B.'S'.M)**(C-B) PERMUTATION
C
C FIRST, PIECE-MEAL ATTACK ON FFT COMPUTATION FOLLOWS
10    IF(MPAS.LE.0) GO TO 40
      IF(MPAS.LT.NPAS)NPAS=MPAS
20    CALL MFLOAD(LREAD,BUFA,IBEX,ICEX,IFLG)
C LOAD CORE WORKING SPACE WITH 2**(C-B) BLOCKS ACCORDING TO MFINDX
      IF(IFLG.LT.0) GO TO 30
      CALL MFCOMP (MEXA,NDIM,ISGN, BUFA,B,C,M, NPAS,IPAS)
C DO MODIFIED IN-CORE FFT, REQUIRING NPAS PASSES STARTING WITH IPAS
      CALL MFLOAD(LWRIT,BUFA,IBEX,ICEX,IFLG)
C UNLOAD CORE AREA, WRITING BLOCKS BACK IN-PLACE TO MASS STORE
      GO TO 20
30    IPAS=IPAS+NPAS
      MPAS=MPAS-NPAS
      GO TO 10
C END OF FIRST PART
C
C SPECIFY BLOCKS TO BE READ IN NEXT PART IN TRUE ORDER (NO PERM)
40    IDUM=MFINDX(LSET,B,M,M,0)
C FINAL, CONCLUDING ATTACK ON FFT COMPUTATION FOLLOWS
50    CALL MFLOAD(LREAD,BUFA,IBEX,ICEX,IFLG)
      IF(IFLG.LT.0)GO TO 80
      CALL MFCOMP (MEXA,NDIM,ISGN, BUFA,C,C,M, C,IPAS)
C DO FINAL, C-PASS IN-CORE FFT OF EACH CORE-LOAD
      IF(SCAL.EQ.1.)GO TO 70
      DO 60 J=1,NC
60    BUFA(J)=BUFA(J)*SCAL
70    CALL MFLOAD(LWRIT,BUFA,IBEX,ICEX,IFLG)
      GO TO 50
C
C BELOW, SORT ARRAY (FULL BIT-REVERSAL AND DIMEN REVERSAL IF IDIR.NE.0)
80    IF(IDIR.EQ.0)GO TO 90
      CALL MFSORT(BUFA,IBEX,ICEX,1,MREAL,MREAL)
      M=MFSUM(MEXA,NDIM,-1)
C DO FULL BIT-REVERSAL OF M BITS (AND REVERSE MEXA LIST)
      RETURN
C
C BELOW, REVERSE BITS OF EACH DIMEN SEPARATELY (NO DIMEN REVERSAL)
90    IH=1
      DO 100 J=1,NDIM
      IG=IH
      IH=IH+MEXA(J)
100   CALL MFSORT(BUFA,IBEX,ICEX,IG,IH,MREAL)
      RETURN
C
      END
      SUBROUTINE MFCOMP (MEXA,NDIM,ISGN, BUFA,B,C,M, NPAS,IPAS)
C
C MODIFIED, IN-CORE FFT OF 2**C ELEMENTS, NPAS PASSES STARTING WITH IPAS
C  MEXA, NDIM, ISGN ,BUFA AND M HAVE SAME MEANING AS IN RMFFT COMMENTS
C  B,C EQUIVALENT TO IBEX,ICEX EXCEPT HERE REFER TO NUM COMPLEX ELMTS
C  (2**C CMPLX ELMTS IN CORE STORE BUFA, IN BLOCKS OF 2**B CMPLX ELMTS)
C
C FFT W PHASE FACTOR COMPUTED RECURSIVELY EXCEPT ON BLOCK BOUNDS
C  MULTIDIMEN FFT ACHIEVED BY REPEATING W SEQUENCES
C
      INTEGER B,C,SPAN,STEP,MEXA(1)
      COMPLEX TEMP,W,D,BUFA(1)
      DATA LINDX/0/,LREST/4/
C
      DATA PI/3.141592653589793/
      PIMOD=PI*2.0**(1-M)
      IF(ISGN.LT.0)PIMOD=-PIMOD
      NC=2**C
C
C BELOW, NPAS COMPUTATION PASSES WHILE DATA IN-CORE
      DO 50 JPAS=1,NPAS
      KPAS=IPAS+JPAS-1
C KPAS IS GLOBAL COMPUTING PASS NUMBER (JPAS-1 IS LOCAL)
      KDIFF=M-MFSUM(MEXA,NDIM,KPAS)
      KPEFF=KPAS+KDIFF
C KPEFF IS EFFECTIVE GLOBAL PASS FOR MULTIDIMEN. FFT W GENERATION
      D=CEXP(CMPLX(0.,PIMOD*2.0**KPEFF))
C D IS USED FOR RECURSIVE MODIFICATION OF W PHASE FACTOR
C
      ITEM=C-JPAS
      SPAN=2**ITEM
      STEP=2*SPAN
C SPAN SEPARATES VALUES IN FFT KERNEL, STEP TO NEXT PAIR, SAME W
C
      IF(B.LT.ITEM)ITEM=B
      NB=2**ITEM
      IF(ITEM.GT.KDIFF)ITEM=KDIFF
      NRPT=2**ITEM
C NRPT COUNTS REPETITION OF W FOR MULTIDIMEN. FFT
C
      IMOD=2**(M-B-KPAS-1)
      IF(IMOD.LE.1)GO TO 20
      IDUM=MFINDX(LREST,0,0,0,0)
      MEXP=KPEFF
      ITEM=KDIFF-B
      IF(ITEM.GT.0)GO TO 10
      MEXP=B+KPAS
      ITEM=0
10    NCLR=2**ITEM
      PIMOD2=PIMOD*2.0**MEXP
C IMOD AND NCLR ARE USED TO COMPUTE W WITHOUT EXCEEDING SMALL INTEGER
C
C BELOW, START OF ONE PASS THROUGH CORE, NOTING BLOCK BOUNDARIES
20    DO 50 I1=1,SPAN,NB
      W=(1.,0.)
      IF(IMOD.LE.1)GO TO 30
      INDWM=MOD(MFINDX(LINDX,0,0,0,0)-1,IMOD)/NCLR
      ANDWM=INDWM
      W=CEXP(CMPLX(0.,PIMOD2*ANDWM))
C NEW W COMPUTED DIRECTLY AT BEGINNING OF NEW BLOCK AREA
C
C BELOW, COMPUTATIONS WITHIN EACH BLOCK OF 2**B CMPLX ELMTS
30    DO 50 I2=1,NB,NRPT
C
C BELOW, REPETITION OF SAME W DUE TO MULTIDIMEN FFT
      DO 40 I3=1,NRPT
      I4=I1+I2+I3-2
C
C BELOW, STEPPING THOUGH INDICES HAVING SAME W IN ONE DIMEN FFT
      DO 40 J=I4,NC,STEP
      K=J+SPAN
      TEMP=(BUFA(J)-BUFA(K))*W
      BUFA(J)=BUFA(J)+BUFA(K)
40    BUFA(K)=TEMP
C FFT 2-POINT KERNEL ARITHMETIC (ALGORITHM BIT-REVERSAL FOLLOWS COMPUT)
C
      W=W*D
C RECURSIVE MODIFICATION OF W WITHIN BLOCK BOUNDARIES
50    CONTINUE
      RETURN
C
      END
      SUBROUTINE MFSORT(BUFA,IBEX,ICEX,IG,IH,M)
C
C BIT-REVERSED PERMUTATION (IG.'R'.IH) OF MASS STORE REAL ARRAY
C  REVERSES IH-IG BITS IN INDEX (M-1,...,IH-1,...,IG,...,0)
C  NOTE THAT THIS IS MORE GENERAL THAN THE FULL M-BIT REVERSAL
C  OF REFERENCE (FRASER, J.ACM, V.23, N.2, APR. 76, P. 306),
C  BUT ALGORITHM IS LOGICALLY THE SAME, WITH ALTERED BIT LIMITS.
C BUFA,IBEX,ICEX AND M HAVE SAME MEANING AS IN COMMENTS IN RMFFT
C  (BLOCKS 2**IBEX, CORE BUFA 2**ICEX, TOTAL 2**M, ALL REAL)
C
      REAL BUFA(1)
      DATA LSET/1/,LPERM/2/,LREAD/1/,LWRIT/2/
C
      IDUM=MFINDX(LSET,IBEX,M,M,0)
C DUMMY CALL TO INITIALISE MFINDX (INITIALLY UNPERMUTED ARRAY)
      IF(IH-IG.LE.1)RETURN
      IF(IG.GE.IBEX)GO TO 50
      IF(IH.LE.ICEX)GO TO 60
C CHECK FOR SPECIAL CASES, REQUIRING SIMPLER TREATMENT
C
C BELOW, MIXED PERMUTATION OF BOTH ELEMENTS AND BLOCKS
      IPAS=0
      NPAS=ICEX-IBEX
C MOST EFFIENT USE OF CORE STORE - TRIES TO DO ICEX-IBEX PASSES PER LOAD
      MPAS=IH-IBEX
      IF(IBEX-IG.LT.MPAS)MPAS=IBEX-IG
C
C BELOW, FIRST VIRTUAL 'S' PERMUTATIONS
10    IF(MPAS.LE.0)GO TO 40
      IF(MPAS.LT.NPAS) NPAS=MPAS
      IGCOR=IG+IPAS
      IF((IGCOR.GT.IBEX-1).AND.(IGCOR.GT.IBEX+NPAS-1))GO TO 30
      IF((IPAS.EQ.0).AND.(IGCOR.GT.IBEX+NPAS-1))GO TO 30
C BYPASS UNNECESSARY CORE LOAD IF TRIVIAL CASES
      IDUM=MFINDX(LPERM,IBEX,IH,M,NPAS)
C DUMMY CALL TO MFINDX TO SPECIFY VIRTUAL (IBEX.S.IH)**NPAS PERM
C
C BELOW, LOAD CORE ACCORDING TO VIRTUAL PERMUTATION AND PERM ELMTS
20    CALL MFLOAD(LREAD,BUFA,IBEX,ICEX,IFLG)
      IF(IFLG.LT.0)GO TO 30
      IF(IPAS.NE.0) CALL MFREV(BUFA,IGCOR,IBEX,ICEX,-1)
      CALL MFREV(BUFA,IGCOR,IBEX+NPAS,ICEX,-1)
C CARRY OUT IN-CORE, SYMMETRIC R PERMS. (ONE ONLY ON FIRST PASS)
      CALL MFLOAD(LWRIT,BUFA,IBEX,ICEX,IFLG)
C UNLOAD CORE AREA, WRITING BLOCKS BACK IN-PLACE TO MASS STORE
      GO TO 20
C
30    IPAS=IPAS+NPAS
      MPAS=MPAS-NPAS
      GO TO 10
C END OF FIRST PART
C
C BELOW, FINAL 'R' PERM. OF BLOCKS IN MASS STORE, IF (IH-2*IBEX+IG).GT.1
40    CALL MFREV(BUFA,0,IH-2*IBEX+IG,M-IBEX,IBEX)
      RETURN
C
C BELOW, PERMUTATION OF BLOCKS ONLY REQUIRED
50    CALL MFREV(BUFA,IG-IBEX,IH-IBEX,M-IBEX,IBEX)
      RETURN
C
C BELOW, PERMUTATION OF ELEMENTS IN CORE ONLY REQUIRED
60    CALL MFLOAD(LREAD,BUFA,IBEX,ICEX,IFLG)
      IF(IFLG.LT.0)RETURN
      CALL MFREV(BUFA,IG,IH,ICEX,-1)
      CALL MFLOAD(LWRIT,BUFA,IBEX,ICEX,IFLG)
      GO TO 60
C
      END
      SUBROUTINE MFREV(BUFA,IG,IH,M,IBEX)
C
C BIT-REVERSED PERMUTATION OF RANDOMLY ADDRESSABLE ELEMENTS
C  REVERSES IH-IG BITS IN INDEX (M-1,...,IH-1,...,IG,...,0)
C  (GENERAL PERM IG.'R'.IH, FRASER, J.ACM, V.23, N.2, APR 1976, P. 300)
C IF IBEX.LT.0, SORTS 2**M REAL ELMTS IN CORE BUFA,
C IF IBEX.GE.0, SORTS BLOCKS IN MASS STORE
C (2**IBEX REAL ELMTS PER BLOCK AND 2**M BLOCKS IN SECOND CASE)
C
C  THE ALGORITHM MAINTAINS A SET OF 'REVERSED' INTEGERS IN ARRAY IRA()
C  OF INCREASING NUMBER OF BITS, UP TO 2 LESS THAN (IH-IG) BITS.
C  INCREMENTING A REVERSED INTEGER THEN REQUIRES THE ALTERNATE
C  ADDITION OF NRA() TO IRA(), OR REPLACEMENT BY THE NEXT LOWER
C  INCREMENTED REVERSED INTEGER IN THE HIERARCHY, RECURSIVELY.
C  THIS IN ITSELF IS FAST, AS RECURSION DEPTHS ARE ON AVERAGE SMALL.
C  BUT, IN ADDITION, ONLY QUARTER LENGTH SERIES ARE GENERATED (-2 BITS)
C  AND THE FULL LENGTH DERIVED BY SCALING BY 2 AND ADDING OFFSETS.
C  IN THIS FINAL STAGE, ONLY VALID SWAP PAIRS ARE GENERATED (1 OR 3 EACH
C
C  WITHIN THE INNER LOOPS, GROUPS OF 2**IG ELMTS ARE MOVED TOGETHER
C  WHILE THIS IS REPEATED OVER 2**(M-IH) PARTS OF THE ARRAY,
C  CORRESPONDING TO THE UNPERMUTED BITS M TO IH AND IG-1 TO 0.
C
      REAL BUFA(1),TEMP
      INTEGER IRA(16),NRA(16),IFOFA(3),IROFA(3)
C
      IHG=IH-IG-3
      IF(IHG.LE.-2)RETURN
C NO PERMUTATION REQUIRED
C
      NB=2**IBEX
      NB1=NB+1
      NG=2**IG
      NGDB=NG*2
      NH=2**IH
      NHHF=NH/2
C NG IS MOVEMENT GROUP SIZE, NH IS PERMUTATION REPLICATION SIZE
      NM=2**M
      NPARS=NM-NH+1
      NREV=NH/4
C
      DO 10 J=1,IHG
      IRA(J)=0
      NREV=NREV/2
10    NRA(J)=NREV
C REVERSED INTEGER RECURSION SETS INITIALISED
C
      NREV=NH/4
      IFOFA(1)=NG-1
      IROFA(1)=NH/2-1
      IFOFA(2)=-1
      IROFA(2)=-1
      IFOFA(3)=NH/2+NG-1
      IROFA(3)=NH/2+NG-1
C THREE PAIRS OF OFFSETS TO CONVERT QUARTER TO FULL LENGTH SERIES
C
      IFOR=0
      IREV=0
C
C BELOW, GENERATE INDEX PAIRS AND SWAP (IREV IS 'TOP' OF IRA() SET)
20    NOF=3
      IF(IFOR.GE.IREV)NOF=1
C SELECTS ONCE-ONLY SWAP PAIRS (EITHER 1 OR 3 PAIRS)
      DO 40 JOF=1,NOF
      IFOF=IFOFA(JOF)
      IROF=IROFA(JOF)
      DO 40 I1=1,NG
C REPETITION OVER GROUP OR SUPER ELEMENT OF NG ACTUAL ELEMENTS
      IN2F=IFOR+IFOF+I1
      IN2R=IREV+IROF+I1
      DO 40 I2=1,NPARS,NH
C REPETITION OF SAME PERMUTATION OVER ARRAY PARTS
      IN3F=IN2F+I2
      IN3R=IN2R+I2
      IF(IBEX.GE.0)GO TO 30
C
C BELOW, IN-CORE ELEMENT SORTING
      TEMP=BUFA(IN3R)
      BUFA(IN3R)=BUFA(IN3F)
      BUFA(IN3F)=TEMP
      GO TO 40
C
C BELOW, SORTING WHOLE BLOCKS IN MASS STORE
30    CALL MFREAD(BUFA,NB,IN3F)
      CALL MFREAD(BUFA(NB1),NB,IN3R)
      CALL MFWRIT(BUFA,NB,IN3R)
      CALL MFWRIT(BUFA(NB1),NB,IN3F)
C
40    CONTINUE
C END OF INNER, REPETITION LOOPS
C
      IFOR=IFOR+NGDB
C INCREMENT FORWARD QUARTER-LENGTH INTEGER (ALREADY SCALED BY NG*2)
      IF(IFOR.GE.NHHF)RETURN
C RETURN FORM SUBROUTINE
C
      IF(IREV.GE.NREV)GO TO 50
C TEST FOR ALTERNATE METHODS OF REVERSE-INCREMENTING (SIMPLE BELOW)
C NOTE THAT REVERSE QUARTER-LENGTH INTEGER IS ALREADY SCALED BY NG*2
C
      IREV=IREV+NREV
      GO TO 20
C
C ALTERNATE RECURSIVE ALTERATION TO QUARTER-LENGTH REVERSED SERIES
50    DO 60 J=1,IHG
      IF(IRA(J).LT.NRA(J))GO TO 70
60    CONTINUE
C
C BELOW, SIMPLE INCREMENT OF REVERSE INTEGER, LOWER IN HIERARCHY
70    IRA(J)=IRA(J)+NRA(J)
      IREV=IRA(J)
80    IF(J.EQ.1)GO TO 20
      J=J-1
      IRA(J)=IREV
      GO TO 80
C
      END
      SUBROUTINE MFLOAD(LOAD,BUFA,IBEX,ICEX,IFLG)
C
C LOADS, UNLOADS CORE STORE ARRAY BUFA, 2**ICEX REALS, 2**IBEX PER BLOCK
C  RETURNS IFLG=+1 NORMALLY, IFLG=-1 WHEN FINISHED ONE PASS OF MASS STOR
C BLOCKS INDEXED ACCORDING TO VIRTUAL PERMUTATION FUNCTION MFINDX
C LOAD=1 (LREAD) READS BLOCKS FROM MASS STORE INTO CORE STORE BUFA
C LOAD=2 (LWRIT) WRITES BLOCKS BACK IN-PLACE TO MASS STORE
C
      REAL BUFA(1)
      DATA LINDX/0/,LHOLD/3/,LREST/4/
C
      NB=2**IBEX
      NCB=2**(ICEX-IBEX)
      IF(LOAD.EQ.2)GO TO 30
      IFLG=+1
      IDUM=MFINDX(LHOLD,0,0,0,0)
C HOLDS CURRENT MFINDX VALUE FOR ENTRY 2 AND SUBRTN MFCOMP
      DO 10 J=1,NCB
      K=(J-1)*NB
      JB=MFINDX(LINDX,0,0,0,0)
      IF(JB.LT.0)GO TO 20
      CALL MFREAD(BUFA(K+1),NB,JB)
C READS BLOCK WITH NEXT VIRTUAL MFINDX INDEX
10    CONTINUE
      RETURN
20    IFLG=-1
      RETURN
C
30    IDUM=MFINDX(LREST,0,0,0,0)
C RESETS MFINDX TO START OF IN-PLACE BLOCK
      DO 40 J=1,NCB
      K=(J-1)*NB
      JB=MFINDX(LINDX,0,0,0,0)
      CALL MFWRIT(BUFA(K+1),NB,JB)
C WRITES BLOCK WITH NEXT VIRTUAL MFINDX INDEX (REPEAT MFREAD SEQUENCE)
40    CONTINUE
      RETURN
C
      END
      FUNCTION MFINDX(LSPEC,B,H,M,N)
C
C VIRTUAL 'S' PERMUTATION (FRASER, J.ACM, V.23, N.2, APR. 76, P.303)
C  CYCLIC SHIFTS H-B BITS IN INDEX (M-1,...,H-1,...,B,...,0)
C COMPUTES NEXT INDEX FOR SEQUENTIAL CORE LOAD, PERM (B.'S'.H)**N
C BLOCK SIZE EXPON B, MASS STORE EXPON M (0.LE.B.LE.H.LE.M)
C N IS EFFECTIVE NUMBER OF LEFT SHIFTS PER I/O PASS (-N RIGHT SHIFTS)
C
C NOTE VARIABLE NAMES AS FOLLOWS:
C     IPERM IS 'P' OF ALGORITHM
C     NPERM=N (ARGUMENT) IS 'N' OF ALGORITHM
C     ISTEP IS 'Q**P' OF ALGORITHM
C     JAY AND KAY ARE 'J' AND 'K' OF ALGORITHM
C NOTE UPPER BOUND H INSTEAD OF M, REQUIRING 2**(M-H) REPEATS
C
C LSPEC=0 (LINDX) RETURNS MFINDX FOR INDEX (B,H,M,N DUMMIES HERE)
C LSPEC=1 (LSET) SETS IPERM=0 (UNPERMED), ENTERS B,H,M,N PARAMS
C LSPEC=2 (LPERM) CHANGES THE B,H,M,N PARAMETERS
C LSPEC=3 (LHOLD) HOLDS CURRENT INDEXING STATE (B,H,M,N DUMMIES HERE)
C LSPEC=4 (LREST) RESTORES STATE TO LAST LHOLD (B,H,M,N DUMMIES HERE)
C
      INTEGER B,H
C
      IF(LSPEC.EQ.1)GO TO 100
      IF(LSPEC.EQ.2)GO TO 200
      IF(LSPEC.EQ.3)GO TO 300
      IF(LSPEC.EQ.4)GO TO 400
C
10    IF(ISTEP.NE.0)GO TO 20
C BELOW, PRECEDES FIRST MFINDX OF A PASS
      IF(NPERM.GT.0)IPERM=MOD(IPERM-NPERM,IHB)
      IF(IPERM.LT.0)IPERM=IHB+IPERM
      ISTEP=2**IPERM
C
C 20 BELOW, NORMAL GENERATION OF NEXT MFINDX
20    MFINDX=JAY+JOFF
      IF(MFINDX.GT.NMB)GO TO 40
      KAY=JAY+ISTEP
      JAY=MOD(KAY,NHB)
      IF(KAY.GE.NHB)JAY=JAY+1
      NRPT=NRPT-1
      IF(NRPT.GT.0)RETURN
C
C NRPT,JOFF REQUIRED TO REPEAT SEQUENCE ON 2**(M-H) PARTS OF ARRAY
      JOFF=JOFF+NHB
30    NRPT=NHB
      JAY=0
      RETURN
C
C 40 BELOW, END OF ONE PASS, PARS RESET, IPERM ALTERED IF INVERSE
40    IF(NPERM.LT.0)IPERM=MOD(IPERM-NPERM,IHB)
      MFINDX=-1
50    JOFF=1
      ISTEP=0
      GO TO 30
C
C LSPEC=1 (LSET) SETS IPERM=0 (UNPERMED), ENTERS B,H,M,N PARAMS
100   IPERM=0
C
C LSPEC=2 (LPERM) CHANGES THE B,H,M,N PARAMETERS (DUMMIES ELSEWHERE)
200   IHB=H-B
      NPERM=N
      NHB=2**IHB
      NMB=2**(M-B)
      MFINDX=IPERM
      GO TO 50
C
C LSPEC=3 (LHOLD) HOLDS CURRENT MFINDX INDEXING PARAMETERS
300   JAYH=JAY
      JOFH=JOFF
      NRPTH=NRPT
310   MFINDX=IPERM
      RETURN
C
C LSPEC=4 (LREST) RESTORES PARAMETERS TO INDEX MFINDX AT LAST LHOLD
400   JAY=JAYH
      JOFF=JOFH
      NRPT=NRPTH
      GO TO 310
C
      END
      FUNCTION MFSUM(MEXA,NDIM,MLIM)
C
C SCANS MEXA LIST IN REVERSE ORDER, RETURNING (MFSUM.JUST GT.MLIM)
C  (IF MLIM LARGE ENOUGH, RETURNS M TOTAL FOR NDIM VALUES)
C  (IF MLIM NEGATIVE, RETURNS M TOTAL, REVERSES ORDER OF MEXA LIST)
C
      INTEGER MEXA(1)
C
      MFSUM=0
      IF(NDIM.LE.0)RETURN
C
      DO 10 J=1,NDIM
      I=NDIM+1-J
      MFSUM=MFSUM+MEXA(I)
      IF((MLIM.GE.0).AND.(MLIM.LT.MFSUM))RETURN
10    CONTINUE
      IF(MLIM.GE.0)RETURN
C
C BELOW, REVERSE ORDER OF MEXA LIST
      NDIMH=NDIM/2
      DO 20 J=1,NDIMH
      K=NDIM+1-J
      MTEM=MEXA(J)
      MEXA(J)=MEXA(K)
20    MEXA(K)=MTEM
      RETURN
C
      END
      SUBROUTINE MFRCMP (MEXA,NDIM,ISGN,IDIR,IPAK, BUFA,B,M)
C
C UNSCRAMBLES REAL-TO-COMPLEX FFT OR VICE-VERSA, CALLED BY SUBRTN RMFFT
C  MOST ARGUMENTS HAVE SAME MEANING AS IN RMFFT COMMENTS
C  BUT 2**B COMPLEX ELMTS IN MASS STORE BLOCK,
C  USES (2**B)*4 CMPLX IN BUFA, 'LOWER', 'UPPER' PLUS EXPANSION AREAS
C  TOTAL MASS STORE ARRAY SIZE OF 2**M COMPLEX ELMTS.
C
      COMPLEX ATEM,BTEM,TEMP,W,D,BUFA(1)
      INTEGER B,JAYA(4),KAYA(4),JWKA(4),KWKA(4),MEXA(1)
C JAYA,KAYA,JWKA,KWKA ALLOW UP TO 4 DIMENSIONS - INCREASE IF REQUIRED
      DATA PI/3.141592653589793/
      DATA LOWER/1/,LUPPR/2/,LCLR/-1/
C
      DO 10 IDIM=1,NDIM
      JAYA(IDIM)=0
      KAYA(IDIM)=0
10    CONTINUE
C MULTIDIMEN. CONJUGATE-SYMMETRIC INDICES ZEROED
C
      IEXPND=1
      NB=2**B
      NBDB=NB*2
      JBOF=2**(M-B)
      MAX=M-B-MEXA(NDIM)
      IF(IDIR*IPAK.LT.0)MAX=M-B
      JBMAX=2**MAX
      IF(MAX.LT.0)JBMAX=1
      IF(IPAK.LT.0)JBMAX=0
C JBMAX IS MAXIMUM BLOCK INDEX REQUIRED (DEPENDS ON IPAK)
C
      IWFG=0
      W=(1.,0.)
      D=CEXP(CMPLX(0.,PI*2.0**(-MEXA(NDIM))))
      IF(ISGN.LT.0)D=CONJG(D)
C W IS COMPLEX PHASE FACTOR, D IS RECURSIVE MODIFIER OF W
C
20    JAY=0
      KAY=0
      IDIM=NDIM
30    IF(IDIM.EQ.1)GO TO 40
      NUMD=2**MEXA(IDIM)
      JWKA(IDIM)=JAY
      KWKA(IDIM)=KAY
      JAY=JAY*NUMD+JAYA(IDIM)
      KAY=KAY*NUMD+KAYA(IDIM)
      IDIM=IDIM-1
      GO TO 30
C CONJUGATE-SYMMETRIC BASE INDICES COMPUTED FROM MULTIDIMEN. SET
C
40    MEX1=MEXA(1)
C 2**MEX1 IS NUMBER OF VALUES ADJACENT IN FIRST DIMENSION
      IFLG=-1
      IF(MEX1.LE.B)GO TO 140
C
C BELOW, FIRST DIMEN. GREATER THAN BLOCK SIZE, MULTIPLE BLOCKS
      NBPD1=2**(MEX1-B)
      NBLCNT=NBPD1
      KINC=NB
      KBINC=NBLCNT-1
      IF(JAY.EQ.KAY)NBLCNT=NBLCNT/2
      JB=JAY*NBPD1+1
      KB=KAY*NBPD1+1
C JB AND KB ARE BLOCK INDEX PAIRS CONTAINING CONJUGATE ELEMENTS
      J1=0
      K1=0
C
50    NCNT=NB+1
60    CALL MFRLOD(LOWER,IOF,BUFA,NB,JB,JBMAX,JBOF,IDIR,NCNT)
C LOWER BLOCK LOADED
      IF(JB.GT.JBMAX)IEXPND=-1
      J2=J1+IOF
      JB=JB+1
      J3=0
      K3=0
      NEWBLK=NCNT-NB
      IF(IFLG.GE.0)GO TO 80
C
70    CALL MFRLOD(LUPPR,IOF,BUFA,NB,KB,JBMAX,JBOF,IDIR,IFLG)
C UPPER BLOCK LOADED
      IFLG=IFLG+1
      K2=K1+IOF
      KB=KB+KBINC
C FIRST TIME, UPPER BLOCK STEPS HIGH, FOLLOWING STEPS SMALL NEGATIVE
      KBINC=-1
C
80    J=J2+J3
      K=K2+K3
C J AND K INDEX CONJUGATE-SYMMETRIC PAIRS IN CORE
      JJ=J+NBDB
      KK=K+NBDB
      IF(IDIR.GT.0)GO TO 180
C
C BELOW, UNSCRAMBLING FOR REAL-TO-COMPLEX FFT
      TEMP=(BUFA(J)+CONJG(BUFA(K)))*0.5
      BTEM=BUFA(K)-CONJG(BUFA(J))
      BTEM=(CMPLX(AIMAG(BTEM),REAL(BTEM)))*0.5*W
      ATEM=TEMP+BTEM
      BTEM=TEMP-BTEM
      BUFA(J)=ATEM
      IF(IEXPND.GT.0)BUFA(JJ)=BTEM
      IF(IWFG.EQ.0)GO TO 150
      BUFA(K)=CONJG(BTEM)
      IF(IEXPND.GT.0)BUFA(KK)=CONJG(ATEM)
C
90    J3=J3+1
      K3=KINC-J3
C IN-CORE INDEX PAIRS STEPPED IN OPPOSING DIRECTIONS
      IF(IDIM.NE.NDIM)GO TO 95
      IWFG=1
      W=W*D
C RECURSIVE MODIFICATION OF W IF UNIDIMEN. TRANSFORM
95    NCNT=NCNT-1
      IF(NCNT.LE.0)GO TO 100
C ENTER RECURSION ROUTINE IF OPERATION COMPLETE IN CURRENT DIMEN
      IF(J3.EQ.1)GO TO 70
      IF(NCNT.GT.NEWBLK)GO TO 80
C END OF INNER LOOP (NOTE SPECIAL CASE WHEN J3=1 ABOVE)
C
C BELOW, MAY REQUIRE TO READ NEW BLOCKS
      NBLCNT=NBLCNT-1
      IF(NBLCNT.GT.0)GO TO 50
      IF(JAY.EQ.KAY)GO TO 60
C JAY.EQ.KAY NEEDS SYMMETRICAL MIDDLE, OTHERWISE CURRENT DIMEN COMPLT
C
C 100 BELOW, RECURSION TO COMPUTE MULTIDIMEN. CONJUGATE-SYMMETRY
100   JAYA(IDIM)=0
      KAYA(IDIM)=0
      IDIM=IDIM+1
      IF(IDIM.GT.NDIM)GO TO 120
      NUMD=2**MEXA(IDIM)
      IF(NUMD.LE.1)GO TO 120
      IF(IDIM.NE.NDIM)GO TO 105
      IWFG=1
      W=W*D
C RECURSIVE MODIFICATION OF W IF MULTIDIMEN. FFT
105   IF(JAYA(IDIM).EQ.0)GO TO 110
      IF((JWKA(IDIM)*NUMD+JAYA(IDIM)).EQ.
     X  (KWKA(IDIM)*NUMD+KAYA(IDIM)))GO TO 100
      IF(KAYA(IDIM).EQ.1)GO TO 100
C
110   JAYA(IDIM)=JAYA(IDIM)+1
      KAYA(IDIM)=NUMD-JAYA(IDIM)
C RECURSIVE STEPPING OF MULTIDIMEN. CONJUGATE-SYMMETRIC INDEX PAIRS
      GO TO 20
C
C BELOW, OPERATION COMPLETE, TIDY UP AND RETURN FROM SUBROUTINE
120   DO 130 IAREA=1,2
      CALL MFRLOD(IAREA,IOF,BUFA,NB,LCLR,JBMAX,JBOF,IDIR,IFLG)
C DUMMY CALL TO MFRLOD TO WRITE ANY UNWRITTEN BLOCKS TO MASS STORE
130   CONTINUE
      RETURN
C RETURN FROM SUBROUTINE
C
C
C 140 BELOW, FIRST DIMEN. LESS THAN BLOCK SIZE, INDEX PAIRS ALL IN-CORE
140   NUMD1=2**MEX1
      NBPD1=NB/NUMD1
      NCNT=NUMD1
      KINC=NCNT
      KBINC=0
      IF(JAY.EQ.KAY)NCNT=NCNT/2+1
      JB=JAY/NBPD1+1
      KB=KAY/NBPD1+1
C JB AND KB ARE BLOCK INDEX PAIRS CONTAINING CONJUGATE ELEMENTS
      J1=(JAY-(JB-1)*NBPD1)*NUMD1
      K1=(KAY-(KB-1)*NBPD1)*NUMD1
      GO TO 60
C
C 150 BELOW, UNSCRAMBLING WITH W0 (IWFG=0) MUST BE TREATED DIFFERENTLY
150   IF(IEXPND.LT.0)GO TO 160
C BELOW, ARRAY EXPANSION (EITHER IPAK=+1 OR IPAK=0 AND STILL REDUNDANT)
      BUFA(K)=CONJG(ATEM)
      BUFA(KK)=CONJG(BTEM)
      GO TO 90
C 160 BELOW, NO ARRAY EXPANSION (EITHER IPAK=-1 OR IPAK=0 NOT REDUNDANT)
160   IF(J.EQ.K)GO TO 170
      BUFA(K)=CONJG(BTEM)
      GO TO 90
C 170 BELOW, SPECIAL CASE IF IPAK=-1 AND ELEMENTS ARE SAME
170   BUFA(J)=CMPLX(REAL(ATEM),REAL(BTEM))
      GO TO 90
C
C 180 BELOW, SCRAMBLING FOR COMPLEX-TO-REAL FFT
180   IF(IWFG.EQ.0)GO TO 200
      BTEM=CONJG(BUFA(K))
190   ATEM=(BUFA(J)+BTEM)
      BTEM=(BUFA(J)-BTEM)*W
      BTEM=CMPLX(AIMAG(BTEM),REAL(BTEM))
      BUFA(J)=ATEM-CONJG(BTEM)
      BUFA(K)=CONJG(ATEM)+BTEM
      GO TO 90
C
C 200 BELOW, SCRAMBLING WITH W0 (IWFG=0) MUST BE TREATED DIFFERENTLY
200   IF(IEXPND.LT.0)GO TO 210
      BTEM=BUFA(JJ)
      GO TO 190
C
C 210 BELOW, NO REDUNDANCY (EITHER IPAK=-1 OR IPAK=0 OR 1 NOT REDUND)
210   IF(J.EQ.K)GO TO 220
      BTEM=CONJG(BUFA(K))
      GO TO 190
C
C 220 BELOW, SPECIAL CASE IF IPAK=-1 AND ELEMENTS ARE SAME
220   BTEM=CMPLX(AIMAG(BUFA(J)),0.)
      BUFA(J)=CMPLX(REAL(BUFA(J)),0.)
      GO TO 190
C
      END
      SUBROUTINE MFRLOD(IAREA,IOF,BUFA,NB,JB,JBMAX,JBOF,IDIR,NCNT)
C
C LOADS, UNLOADS CORE STORE ARRAY BUFA, FOR REAL FFT UNSCRAMBLING ROUTIN
C  BLOCK SIZE NB CMPLX, BLOCK NUMBER JB (JB=-1 DOES FINAL TIDY O/P)
C  JBMAX IS MAX BLOCK INDX FOR EXPANSN, JBOF OFFSET TO EXPANDING BLOCKS
C  IDIR=-1 DIRECTION REAL/CMPLX, +1 CMPLX/REAL
C IAREA=1 (LOWER) OR 2 (UPPER) OF TWO AREAS IN LOGICAL UNSCRAMBLING
C  NOTE THAT BLOCK NORMALLY PHYSICALLY LOADED IN THESE AREAS,
C  BUT, IF BLOCK ALREADY RESIDENT, MAY BE IN DIFFERENT AREA, SO
C IOF RETURNED AS ACTUAL OFFSET IN BUFA TO LOADED BLOCK.
C  USES (2**B)*4 CMPLX IN BUFA, 'LOWER', 'UPPER' PLUS EXPANSION AREAS
C  RETURNS IOF AS BUFFER OFFSET TO AREA (MAY NOT BE SAME, IF BLOCK RESID
C
C NCNT IS COUNT OF ELEMENTS TO BE ACCESSED IN THIS LOAD, TO ALLOW
C  NOTE TO BE TAKEN OF ANY PARTLY FILLED BLOCKS DURING EXPANSION,
C  PREVENTING THE READING OF 'NON-EXISTENT' BLOCKS.
C LISTS JBPFA(), NCPFA() OF SIZE NPARF HOLD THIS INFORMATION, DEFAULTS
C  TO ALL-READ IF EXCEEDED, BUT INCREASE NPARF ETC, IF PROBLEM.
C
      COMPLEX BUFA(1)
      INTEGER JBAREA(2),NCAREA(2),JBPFA(5),NCPFA(5)
      DATA JBAREA(1)/-1/,JBAREA(2)/-1/,NPARF/5/
C
      IF(JBAREA(1).GE.0)GO TO 20
C JBAREA() HOLDS INDEX OF BLOCK LOADED IN AREA 1 OR 2 (FIRST TIME -1 BEL
C
      IEXIST=-1
      DO 10 I=1,NPARF
10    JBPFA(I)=-1
C PRE-CLEARS PARTLY FILLED BLOCK LIST (ONCE BLOCK FILLED, ALSO CLEARED)
C
20    NBDB=NB*2
      IF(JB.LT.0)GO TO 50
      IF(IAREA.EQ.2)GO TO 30
      NCLOW=NCNT
      IF(MOD(NCNT,2).NE.0)NCLOW=(NCLOW-1)*2
30    NCHLD=NCLOW
      IF(IAREA.EQ.2)NCHLD=NCLOW-1
      IF(NCNT.LT.0)NCHLD=1
C NCHLD IS THE NUMBER OF ELEMENTS TO BE ACCESSED IN CURRENT READ
C
      DO 40 I=1,2
      IF(JB.EQ.JBAREA(I))GO TO 140
40    CONTINUE
C TEST DONE TO SEE IF REQUIRED BLOCK ALREADY IN CORE (TRIVIAL IF SO)
C
C OTHERWISE BELOW, FIRST WRITE OUT RESIDENT BLOCK, THEN READ IN NEW
50    IOF=(IAREA-1)*NB+1
C IOF IS BASE OFFSET OF CORE AREA WHERE BLOCK IS TO BE FOUND
      IOFDB=IOF+NBDB
      IF(JBAREA(IAREA).LT.0)GO TO 90
      CALL MFWRIT(BUFA(IOF),NBDB,JBAREA(IAREA))
C WRITE OUT BLOCK BEFORE READING NEW BLOCK
      IF(IDIR.GT.0)GO TO 90
      IF(JBAREA(IAREA).GT.JBMAX)GO TO 90
      IF(NCAREA(IAREA).GE.NB)GO TO 80
C
C BELOW, IF LAST BLOCK ONLY PART-FILLED, INDEX, ELMTS ACCESSED NOTED
      DO 60 I=1,NPARF
      IF(JBPFA(I).LT.0)GO TO 70
60    CONTINUE
C NO ROOM IN LISTS, DEFAULTS TO ALL READ
      I=NPARF
      IEXIST=I
70    JBPFA(I)=JBAREA(IAREA)
      NCPFA(I)=NCAREA(IAREA)
C
80    CALL MFWRIT(BUFA(IOFDB),NBDB,JBOF+JBAREA(IAREA))
C SIMILARLY, WRITE OUT BLOCK PAIR IF EXPANDING
C
C BELOW, READ BLOCK NOTING BLOCK INDX (READ EXPANDED, IF PART FILLED)
90    JBAREA(IAREA)=JB
      IF(JB.LT.0)GO TO 130
      CALL MFREAD(BUFA(IOF),NBDB,JB)
C READ REQUIRED BLOCK AND NOTE ACCESS COUNT
      NCAREA(IAREA)=NCHLD
      IF(JB.GT.JBMAX)GO TO 130
      IF(IDIR.GT.0)GO TO 120
C
C BELOW, EXPANSION - DOES BLOCK EXIST TO READ
      DO 100 I=1,NPARF
      IF(JB.EQ.JBPFA(I))GO TO 110
100   CONTINUE
      IF(IEXIST.LT.0)GO TO 130
110   JBPFA(I)=-1
      NCAREA(IAREA)=NCPFA(I)+NCHLD
C IF BLOCK TO BE READ WAS ONLY PART FILLED, THEN IT EXISTS TO READ
120   CALL MFREAD(BUFA(IOFDB),NBDB,JBOF+JB)
C READ EXPANDED BLOCK IF REQUIRED
C
130   RETURN
C RETURN FROM SUBROUTINE
C
C BELOW, TRIVIAL CASE - BLOCK ALREADY LOADED
140   IOF=(I-1)*NB+1
C IOF IS BASE OFFSET OF CORE AREA WHERE BLOCK IS TO BE FOUND
      IF(I.NE.IAREA)GO TO 150
      NCAREA(I)=NCAREA(I)+NCHLD
C INCREASE ACCESS COUNT IF CURRENT IAREA MATCHES ORIGINAL IAREA
150   RETURN
C
      END
      FUNCTION MFPAR(IRMF,ICOMP)
C
C  HELPER ROUTINE TO CROSS-COMPUTE MASS STORE FFT FILE PARAMETERS
C  PARAMETERS ARE HELD AND COMPUTED IN 3 COMMON AREAS (SEE BELOW)
C  MFPAR RETURNS 0 NORMALLY, -1 IF NOT ALL MFINT CORRECT, +1 IBEX ERROR
C
C  COMMON/MFARG/ HOLDS ARGUMENTS AS USED IN FFT CALLS, AS FOLLOWS:
C  VARIABLE NAMES HAVE SAME MEANING AS COMMENTS, SUBROUTINE RMFFT
C     MEXA() HOLDS EXPONENTS FOR UP TO 4 DIMENSIONS (R/T ZEROS EXCESS)
C     NDIM NUM DIMENS, IBEX,ICEX BLOCK AND CORE EXPONS, IPAK RMFFT PACKI
C     ISGN,IDIR,SCAL ARE IGNORED HERE, BUT INCLUDED FOR COMPLETENESS
C
C  COMMON/MFVAL/,/MFINT/ RETURN COMPUTED VALUES, USEFUL FOR FILE ACCESS
C     NDMA(4),DIMA(4) HOLD DIMENSION SIZES CORRESPONDING TO MEXA()
C     (EG. NDMA(1)=DIMA(1)=2.**MEXA(1), ETC. AND =1. BEYOND NDIM)
C
C     NTD1,TDM1  IS CURRENT TOTAL NUM OF 'RECDS' OF SIZE NRD1,RDM1
C     NRD1,RDM1  IS NUM OF REALS IN CURRENT FIRST DIMENSION
C     (USEFUL FOR ACCESSING DATA BY MFREAD/MFWRIT, NRD1,RDM1 REALS,
C      ASSUMING THAT MFREAD/MFWRIT CAN HANDLE 'RECDS' OF DIFFERENT SIZES
C
C     NFBK,FBLK  IS MAXIMUM FILE SIZE OF 'RECDS' OF SIZE NRBK,RBLK
C     NTBK,TBLK  IS CURRENT TOTAL NUM OF 'RECDS' OF SIZE NRBK,RBLK
C     NRBK,RBLK  IS NUM OF REALS IN FFT WORKING BLOCK (2**IBEX REALS)
C     (GIVES MAX AND CURRENT FILE SIZE AND ACCESS BY FFT ROUTINES,
C   NFBK.GT.NTBK ONLY WITH PACKED REAL DATA WHEN EXPANDING, IPAK=0 OR 1)
C
C     NRCR,RCOR  IS NUM OF REALS IN FFT WORKING CORE  (2**ICEX REALS)
C     NSZE=SIZE=2.**M, WHICH IS THE EFFECTIVE TOTAL SIZE OF TRANSFORM,
C        WHERE M IS SUM TO NDIM OF MEXA() (SEE RMFFT COMMENTS)
C
C  NOTE THAT ALL /MFARG/ ARE INTGS (EXCEPT SCAL), ALL /MFVAL/ REALS
C  (/MFINT/ IS INTEGER CONVERSION OF /MFVAL/, ANY VALUE OF MFINT
C   IS SET -1 IF TOO LARGE, BY FIXMAX, AND MFPAR RETURNED -1 AS FLAG,
C   FIXMAX SET BY DATA STATEMENT TO 32767. HERE, BUT ALTER TO SUIT)
C
C ROUTINE ARGUMENTS HAVE THE FOLLOWING EFFECT:
C  IRMF=-1, DATA IS PACKED REAL, +1 DATA IS COMPLEX, ROUTINE RMFFT,
C  IRMF=0,  DATA IS COMPLEX, ROUTINE CMFFT
C
C  ICOMP=0, COMPUTES VALUES IN /MFVAL/ FROM VALUES GIVEN IN /MFARG/
C  ICOMP=1 , REVERSE COMPUTES EXPONENTS IN /MFARG/ FROM /MFVAL/
C     (DIMA(),RBLK,RCOR GIVEN INSTEAD OF MEXA(),IBEX,ICEX)
C
C  NOTE, ROUTINE FORCES ICEX, IBEX TO CORRECT RANGE, MFPAR=+1 IF CANNOT
C
      COMMON/MFARG/MEXA(4),NDIM,ISGN,IDIR,SCAL,IBEX,ICEX,IPAK
      COMMON/MFVAL/DIMA(4),TDM1,RDM1,FBLK,TBLK,RBLK,RCOR,SIZE
      COMMON/MFINT/NDMA(4),NTD1,NRD1,NFBK,NTBK,NRBK,NRCR,NSZE
      REAL VAL(11)
      INTEGER INT(11)
      EQUIVALENCE (VAL(1),DIMA(1)),(INT(1),NDMA(1))
C
      DATA FIXMAX/32767./,NMAX/4/
C
      MFPAR=0
      IF(ICOMP.EQ.0)GO TO 20
      ALG2=ALOG(2.)
      IBEX=IFIX(ALOG(RBLK)/ALG2+0.5)
      ICEX=IFIX(ALOG(RCOR)/ALG2+0.5)
C
      DO 10 I=1,NDIM
10    MEXA(I)=IFIX(ALOG(DIMA(I))/ALG2+0.5)
C
20    M=0
      DO 30 I=1,NMAX
      IF(I.GT.NDIM)MEXA(I)=0
      M=M+MEXA(I)
30    DIMA(I)=2.**MEXA(I)
      SIZE=2.**M
C
      IF(IRMF.EQ.0)GO TO 90
      IF(ICEX.GT.M)ICEX=M
      IF(IBEX.GT.ICEX-2)IBEX=ICEX-2
      IF(IBEX.LT.2)MFPAR=1
C FORCES ICEX.NGT.M AND IBEX.NGT.ICEX-2, OR MFPAR=1 (IRMF=+/- 1)
C
40    RBLK=2.**IBEX
      RCOR=2.**ICEX
C
      FADD=SIZE
      IF(IRMF.EQ.0.OR.IPAK.GT.0)GO TO 50
C FADD IS ADDITIONAL FILE SIZE IN REALS, 'SIZE' IF CMFFT OR IPAK=1
C
      FADD=0.
      IF(IPAK.LT.0)GO TO 50
C IPAK=-1 REQUIRES NO FILE EXPANSION
C
      IDIM=NDIM
      IF(IRMF.LT.0)IDIM=1
      FADD=SIZE*2./DIMA(IDIM)
C FADD COMPUTED FOR PARTICULAR CASE OF IPAK=0, WHEN COMPLX
C
50    FSIZ=SIZE+FADD
      ITEM=IFIX(FSIZ/RBLK+0.5)
      IF(FLOAT(ITEM)*RBLK+0.5 .LT.FSIZ)ITEM=ITEM+1
      FBLK=FLOAT(ITEM)
C FBLK IS MAXIMUM NUMBER OF 'RECDS', SIZE RBLK, POSSIBLE
      TBLK=FBLK
      IF(IRMF.GE.0)GO TO 60
C GENERALLY TBLK=FBLK, BUT FOR PACKED REAL NOT SO, BELOW
      FSIZ=SIZE
      TBLK=FSIZ/RBLK
C
60    TDM1=1.
      RDM1=FSIZ
      IF(NDIM.EQ.1)GO TO 70
C JOB COMPLETED IF NDIM=1
C
      RDM1=DIMA(1)
      IF(IRMF.GE.0)RDM1=RDM1*2.
      TDM1=FSIZ/RDM1
C OTHERWISE COMPUTE TDM1 AS NUMBER OF 'RECDS', SIZE RDM1 REALS
C
70    DO 80 I=1,11
      INT(I)=-1
      IF(VAL(I).LE.FIXMAX)INT(I)=IFIX(VAL(I)+0.5)
      IF(INT(I).LT.0.AND.MFPAR.EQ.0)MFPAR=-1
80    CONTINUE
C CONVERT VALUES IN /MFVAL/ TO INTEGERS IN /MFINT/ (-1 IF TOO LARGE)
      RETURN
C
90    IF(ICEX.GT.M+1)ICEX=M+1
      IF(IBEX.GT.ICEX-1)IBEX=ICEX-1
      IF(IBEX.LT.1)MFPAR=1
      GO TO 40
C  FORCES ICEX.NGT.M+1 AND IBEX.NGT.ICEX-1, OR MFPAR=1 (IRMF=0)
C
      END
      SUBROUTINE DMPERM (MEXA,NDIM,NSHFT,IREX, BUFA,IBEX,ICEX)
C
C SHIFTS ORDER OF DIMENSIONS OF REAL OR COMPLEX MASS STORE ARRAY
C NOTE, THIS IS NOT USED BY FFT SUBRTNS BUT IS INCLUDED FOR COMPLETENESS
C  (FRASER, ACM TOMS - 1978/79, AND J.ACM, V.23,N.2, APRIL 76, PP. 298-3
C
C MEXA(J) LIST OF DIMENSION SIZE EXPONS (BASE 2), ADJACENT VARIABLES FIR
C NDIM IS NUMBER OF EXPONENTS IN LIST AND THUS THE NUMBER OF DIMENSIONS
C  SUM TO NDIM: MEXA(J)=M, WHERE 2**M IS SIZE OF MASS STORE ARRAY (SEE B
C NSHFT IS DIMENSION SHIFT COUNT, THUS:
C  NSHFT=0,  NO SHIFT OR CHANGE OCCURS
C  NSHFT=1,2 ETC., FIRST TO NEXT DIMENSION, CIRC NSHFT PLACE SHIFT (MOD
C  NSHFT=-1, REVERSES THE ORDER OF DIMENSIONS
C IREX=0 REAL, 1 COMPLEX (THAT IS, MOVEMENT GROUP IS 2**IREX REALS,
C    AND TOTAL MASS STORE SIZE IS 2**(M+IREX) REAL ELEMENTS)
C
      REAL BUFA(1)
      INTEGER MEXA(1)
C
      NS=MOD(NSHFT,NDIM)
      IF(NS.EQ.0)RETURN
      M=MFSUM(MEXA,NDIM,-1)+IREX
C FINDS M TOTAL AND REVERSES MEXA LIST
      CALL MFSORT(BUFA,IBEX,ICEX,IREX,M,M)
C INITIAL OVERALL BIT-REVERSAL M BITS ABOVE IREX BITS
      IF(NSHFT.LT.0)GO TO 10
C
C BELOW, REVERSAL OF TWO PARTS, TO FORM REQUIRED SHIFT
      IH=MFSUM(MEXA,NS,-1)+IREX
      CALL MFSORT(BUFA,IBEX,ICEX,IREX,IH,M)
C REVERSE LOWER PART OF MEXA LIST AND LOWER PART OF ARRAY BITS
      CALL MFSORT(BUFA,IBEX,ICEX,IH,M,M)
C SEPARATELY REVERSE UPPER PART OF ARRAY BITS
      IH=MFSUM(MEXA(NSHFT+1),NDIM-NS,-1)
C REVERSE UPPER PART OF MEXA LIST
      RETURN
C RETURN FROM SUBROUTINE AFTER CYCLIC SHIFTS
C
C BELOW, SEPARATELY REVERSE OVER EACH DIMENSION (DIMEN REVERSAL)
10    IH=IREX
      DO 20 J=1,NDIM
      IG=IH
      IH=IH+MEXA(J)
20    CALL MFSORT(BUFA,IBEX,ICEX,IG,IH,M)
      RETURN
C
      END
C
C THIS PROGRAM TESTS THE MASS STORE FFT BY COMPARISON WITH NAIVE DFT
C     FFT PARAMETERS MAY BE ALTERED AT WILL (SEE COMMENTS)
C MASS STORE IS SIMULATED BY FORTRAN ARRAYS (SEE DUMMY I/O SUBROUTINES)
C PRINTING MAY BE COPIOUS, OR ONLY MAX DIFFERENCES (SEE COMMENTS)
C TEST OK IF MAX DIFFERENCES ARE NEAR ORDER OF MACHINE ROUND-OFF
C
C  NOTE WELL THAT TYPE COMPLEX DATA MUST EXIST AS ALTERNATING
C  REAL/IMAG/REAL/IMAG... ELEMENTS, BOTH IN MASS STORE AND IN FORTRAN
C  WORKING ARRAY BUFA;  IN THIS FFT, DIFFERENT SUBROUTINES WILL SET
C  DIFFERENT TYPE (REAL OR COMPLEX) FOR ARRAY BUFA.
C
      COMMON/MFARG/MEXA(4),NDIM,ISGN,IDIR,SCAL,IBEX,ICEX,IPAK
      COMMON/MFVAL/DIMA(4),TDM1,RDM1,FBLK,TBLK,RBLK,RCOR,SIZE
      COMMON/MFINT/NDMA(4),NTD1,NRD1,NFBK,NTBK,NRBK,NRCR,NSZE
C
C COMMON AREAS /MFARG/,/MFVAL/,/MFINT/ USEFUL FOR RUNNING MASS STORE FFT
C  /MFARG/ HOLDS ARGUMENTS USED IN FFT CALLES, MOSTLY EXPONENTS
C    HELPER ROUTINE MFPAR COMPUTES VALUES FROM /MFARG/ INTO
C  /MFVAL/ (REALS AS SOME LARGE), /MFINT/ (INTEGER EQUIVALENTS IF POSS)
C    OR CAN REVERSE-COMPUTE SOME /MFARG/ FROM /MFVAL/
C       SEE COMMENTS, ROUTINE MFPAR.
C
      COMMON/MASS/RMAS(1024)
      COMPLEX ANAIV(512),BNAIV(512),CBUFA(512),CDIF
      REAL RANDA(1024),BUFA(1024)
      EQUIVALENCE (CBUFA(1),BUFA(1))
C
C COMMON/MASS/RMAS() USED BY ROUTINES MFREAD/MFWRIT TO SIMULATE MASS STORE
C  ANAIV(), BNAIV()  HOLD RESULT OF NAIVE DFT FOR COMPARISON
C  BUFA()=CBUFA()    IS WORKING AREA IN CORE STORE FOR FFT AND PROGRAM
C  RANDA  HOLDS PSEUDO RANDOM DATA USED IN TEST
C
      LP=5
      IPRINT=0
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 FOR COPIOUS PRINT,
C  IPRINT=0 FOR MAX DIFFERENCES ONLY, -1 OVERALL MAX DIFFERENCE ONLY
C
      M=5
C M SETS THE OVERALL ARRAY SIZE FOR AUTO IBEX,ICEX,MEXA,NDIM STEPPING
C  FIXED VALUES CAN BE USED (SEE COMMENTS BELOW DO 100 STATEMENTS)
C
      IRMF=-1
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN=-1
      IDIR=-1
      IPAK=1
C FFT ARGS, ISGN=+/- 1, IDIR=-1 (RMFFT), IDIR=1 OR 0 (CMFFT)
C IPAK=1 OR 0 (RMFFT), -1 GIVES APPARENT FAILURES DUE TO SQUEESED RESULT
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      IF(IPRINT.LE.0)WRITE(LP,910)IPRINT,IRMF
910   FORMAT(31H1MASS STORE FFT TEST - IPRINT =,I3,
     X  36H (1=COPIOUS,0=MAX DIFFS,-1=OVERALL),,
     X  9H   IRMF =,I3,28H (0=CMFFT TEST,1=RMFFT TEST) /)
      DIFMG=0.
C
C BELOW, DO 100 COMPUTES ALL POSSIBLE MEXA FOR 1,2 AND 3 DIMENSIONS
      DO 100 NDIM=1,3
      M2M=M-NDIM+1
      IF(NDIM.EQ.1)M2M=1
      DO 100 M2=1,M2M
      M3M=M2M-M2+1
      IF(NDIM.NE.3)M3M=1
      DO 100 M3=1,M3M
      IF(NDIM.EQ.1)MEXA(1)=M
      IF(NDIM.EQ.2)MEXA(1)=M-M2
      IF(NDIM.EQ.3)MEXA(1)=M-M2-M3
      MEXA(2)=M2
      MEXA(3)=M3
C REPLACE DO 100 SET BY FIXED MEXA LIST, IF DESIRED
C
      IBEX=2
      ICEX=4
C DUMMY IBEX, ICEX SO NO ERROR IN MFPAR BELOW
      IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
C CALL HELPER ROUTINE TO COMPUTE SIZES USED BY NAIVE DFT SUBRTN
C
      M=MFSUM(MEXA,NDIM,99)
      AR=RANMF(1)
C RESET RANDOM NUMBER GENERATOR AND COMPUTE M IN CASE NOT GIVEN
      DO 10 J=1,NSZE
      RANDA(J)=RANMF(-1)
10    ANAIV(J)=CMPLX(RANDA(J),0.)
C LOAD RANDOM NUMBERS FOR FFT AND FOR NAIVE SUBROUTINE
C
      CALL NAIVE(ANAIV,BNAIV,NDMA(1),NDMA(2),NDMA(3),ISGN)
C NAIVE DFT SUBROUTINE CALLED TO COMPUTE 'SLOW' FOURIER TRANSFORM
C
C
      IF(IRMF.EQ.0)GO TO 20
C BELOW, SET LIMTS FOR STEPPING IBEX, ICEX, WHEN CALLING RMFFT
      IBEXL=2
      ICEXL=4
      ICEXM=M
      GO TO 30
C
C SETS DIFFERENT LIMITS FOR IBEX, ICEX FOR CMFFT BELOW (IRMF.EQ.0)
20    IBEXL=1
      ICEXL=2
      ICEXM=M+1
C
30    IF(ICEXL.GT.ICEXM)GO TO 1000
      DO 100 ICEX=ICEXL,ICEXM
      IBEXM=ICEX-2
      IF(IRMF.EQ.0)IBEXM=ICEX-1
      DO 100 IBEX=IBEXL,IBEXM
C IBEX AND ICEX COMPUTED; REPLACE DO 100 ABOVE BY FIXED IBEX,ICEX IF REQUD.
C
      IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1
C HELPER ROUTN, NTD1 TOTAL NUMB OF NRD1 REALS IN FIRST DIMENSION
C   (NCNT IS NUMBER OF REAL ELMTS IN FIRST DIMENSION)
      SCAL=1.
      IF(IRMF.EQ.0)GO TO 50
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 40 JB=1,NTD1
      K=(JB-1)*NCNT
      DO 35 I=1,NCNT
      J=K+I
35    BUFA(I)=RANDA(J)
      CALL MFWRIT(BUFA,NRD1,JB)
40    CONTINUE
C LOAD RANDOM NUMBERS IN REAL ARRAY IN 'RECDS' OF FIRST DIMEN LENGTH
C  (NOTE, THIS REQUIRES MFREAD/MFWRIT TO BE ABLE TO ACCEPT 'RECORDS'
C   OF DIFFERENT LENGTH; OTHERWISE MUST USE NTBK AND NRBK HERE)
C
      CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
      GO TO 70
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FIRST DIMEN
50    NCNT=NCNT/2
      DO 60 JB=1,NTD1
      K=(JB-1)*NCNT
      DO 55 I=1,NCNT
      J=K+I
55    CBUFA(I)=CMPLX(RANDA(J),0.)
      CALL MFWRIT(BUFA,NRD1,JB)
60    CONTINUE
C LOAD REAL VALUES IN CMPLX ARRAY IN 'RECDS' OF FIRST DIMEN LENGTH
C  (NOTE, THIS REQUIRES MFREAD/MFWRIT TO BE ABLE TO ACCEPT 'RECORDS'
C   OF DIFFERENT LENGTH; OTHERWISE MUST USE NTBK AND NRBK HERE)
C
      CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
70    IF(IPRINT.LE.0)GO TO 75
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
920   FORMAT(8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      WRITE(LP,930)
930   FORMAT(8H   INDEX,12X,3HFFT,24X,5HNAIVE,24X,4HDIFF)
C
75    IERR=MFPAR(-IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1/2
C HELPER ROUTN, NTD1 TOTAL NUMB OF NRD1 REALS IN FIRST DIMENSION
C  (NOTE IRMF=+1 FOR DATA IN COMPLEX STATE, NCNT ELMTS)
C
C BELOW, COMPARES FFT COMPLEX RESULT WITH NAIVE RESULT
      DIFM=0.
      DO 80 JB=1,NTD1
      K=(JB-1)*NCNT
      CALL MFREAD(BUFA,NRD1,JB)
C READ 'RECDS' OF FIRST DIMEN LENGTH (RECOMPUTED ABOVE BY MFPAR)
C  (NOTE, THIS REQUIRES MFREAD/MFWRIT TO BE ABLE TO ACCEPT 'RECORDS'
C   OF DIFFERENT LENGTH; OTHERWISE MUST USE NTBK AND NRBK HERE)
C
      DO 80 I=1,NCNT
      J=K+I
      INDEX=J-1
      IF(IDIR.NE.0)CDIF=CBUFA(I)-BNAIV(J)
      IF(IDIR.EQ.0)CDIF=CBUFA(I)-ANAIV(J)
      DIF=ABS(REAL(CDIF))
      IF(DIF.GT.DIFM)DIFM=DIF
      DIF=ABS(AIMAG(CDIF))
      IF(DIF.GT.DIFM)DIFM=DIF
      IF(IPRINT.LE.0)GO TO 80
      IF(IDIR.NE.0)WRITE(LP,940)INDEX,CBUFA(I),BNAIV(J),CDIF
      IF(IDIR.EQ.0)WRITE(LP,940)INDEX,CBUFA(I),ANAIV(J),CDIF
940   FORMAT(1X,I5,3(2X,2E13.4))
80    CONTINUE
C BELOW, PRINT INTERMEDIATE DIFFERENCES
C
      IF(IPRINT.GE.0)WRITE(LP,950)DIFM,NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
950   FORMAT(10H MAX DIFF ,E11.4,1X,
     X  8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      IF(DIFM.GT.DIFMG)DIFMG=DIFM
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
      ISGN=-ISGN
      IDIR=-IDIR
      SCAL=1./SIZE
      IF(IRMF.NE.0)
     X  CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF(IRMF.EQ.0)
     X  CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF(IPRINT.LE.0)GO TO 85
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
      WRITE(LP,960)
960   FORMAT(8H   INDEX,4X,8HFFT/2**M,6X,5HINPUT,10X,4HDIFF)
C
85    IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1
      IF(IRMF.EQ.0)NCNT=NCNT/2
C HELPER ROUTN, NTD1 TOTAL NUMB OF NRD1 REALS IN FIRST DIMENSION
C   (NCNT IS NUMBER OF ELMTS IN FIRST DIMENSION, REAL OR CMPLX)
C
C BELOW, COMPARES FFT INVERSE RESULT WITH INITIAL RANDOM INPUT
      DIFM=0.
      DO 90 JB=1,NTD1
      K=(JB-1)*NCNT
      CALL MFREAD(BUFA,NRD1,JB)
C READ 'RECDS' OF FIRST DIMEN LENGTH (RECOMPUTED ABOVE BY MFPAR)
C  (NOTE, THIS REQUIRES MFREAD/MFWRIT TO BE ABLE TO ACCEPT 'RECORDS'
C   OF DIFFERENT LENGTH; OTHERWISE MUST USE NTBK AND NRBK HERE)
C
      DO 90 I=1,NCNT
      J=K+I
      INDEX=J-1
      IF(IRMF.NE.0)RM=BUFA(I)
      IF(IRMF.EQ.0)RM=REAL(CBUFA(I))
      DIF=ABS(RM-RANDA(J))
      IF(DIF.GT.DIFM)DIFM=DIF
      IF(IPRINT.GT.0)WRITE(LP,940)INDEX,RM,RANDA(J),DIF
90    CONTINUE
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RANDOM SET)
C
      IF(IPRINT.GE.0)WRITE(LP,950)DIFM,NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
      IF(DIFM.GT.DIFMG)DIFMG=DIFM
      ISGN=-ISGN
      IDIR=-IDIR
C RESTORE ISGN AND IDIR FOR FORWARD TRANSFORM
100   CONTINUE
C
C BELOW, PRINT OVERALL MAXIMUM DIFFERENCE
      WRITE(LP,970)DIFMG,M,ISGN,IDIR,IPAK
970   FORMAT(18H OVERALL MAX DIFF ,E11.4,3X,
     X 48H FOR ALL MEXA(1 TO 3 DIM),IBEX,ICEX,MEXA FOR M =,I3,
     X 19H,  ISGN,IDIR(+/-) =,2I4,8H  IPAK =,I4)
      STOP
C
1000  WRITE(LP,980)IERR
980   FORMAT(25H IBEX FORCED TOO SMALL OR,
     X   30H NOT ALL /MFINT/ CORRECT, IERR,I4)
      STOP
      END
      FUNCTION RANMF(J)
C
C RANDOM NUMBER GENERATOR FOR MASS STORE FFT TEST
C
      IF(J.GE.0)GO TO 20
C POSITIVE J CAUSES RESET OF INITIAL K
C NEGATIVE J MUST BE USED NORMALLY
C
      MODULO=2048
      FLMOD=2048.0
      DO 10 I=1,15
10    K=MOD(5*K,MODULO)
      Z=FLOAT(K)/FLMOD
      RANMF=Z
      RETURN
C
20    K=J
      RANMF=J
      RETURN
C
      END
      SUBROUTINE NAIVE(ANAIV,BNAIV,NJ,NK,NL,ISGN)
C
C NAIVE DISCRETE FOURIER TRANSFORM - 1 TO 3 DIMENSIONS
C  USED TO TEST MASS STORE FFT, INPUT ARRAY  ANAIV(NJ,NK,NL)
C RESULT RETURNED IN BOTH ARRAYS ANAIV AND BNAIV
C  ANAIV DIMENSIONS IN INITIAL ORDER, BNAIV REVERSED ORDER
C  NJ,NK,NL DIMENSIONING, ISGN SIGN OF COMPLEX EXPONENT OF FOURIER
C
      COMPLEX TEMP,ANAIV(NJ,NK,NL),BNAIV(NL,NK,NJ)
      DATA PI/3.141592653589793/
C
      PI2=PI*2.0
      IF(ISGN.LT.0)PI2=-PI2
C
      DO 20 JB=1,NJ
      AJB=FLOAT(JB-1)/FLOAT(NJ)
      DO 20 KB=1,NK
      AKB=FLOAT(KB-1)/FLOAT(NK)
      DO 20 LB=1,NL
      ALB=FLOAT(LB-1)/FLOAT(NL)
C
      TEMP=(0.,0.)
      DO 10 JA=1,NJ
      AJA=FLOAT(JA-1)*AJB
      DO 10 KA=1,NK
      AKA=FLOAT(KA-1)*AKB
      DO 10 LA=1,NL
      ALA=FLOAT(LA-1)*ALB
10    TEMP=TEMP+ANAIV(JA,KA,LA)*CEXP(CMPLX(0.,PI2*(AJA+AKA+ALA)))
C
20    BNAIV(LB,KB,JB)=TEMP
C
      DO 30 JA=1,NJ
      DO 30 KA=1,NK
      DO 30 LA=1,NL
30    ANAIV(JA,KA,LA)=BNAIV(LA,KA,JA)
      RETURN
C
      END
      SUBROUTINE MFREAD(BUFA,NB,JB)
C
C DUMMY SUBROUTINE TO SIMULATE RANDOM ACCESS MASS STORE READ
C READ BLOCK, INDEX JB, FROM MASS STORE TO BUFA, NB REAL VALUES
C COMMON ARRAY RMAS SIMULATES MASS STORE ARRAY
C
      COMMON/MASS/RMAS(1024)
      REAL BUFA(NB)
C
      IOF=(JB-1)*NB
      DO 10 I=1,NB
      K=IOF+I
10    BUFA(I)=RMAS(K)
      RETURN
C
      END
      SUBROUTINE MFWRIT(BUFA,NB,JB)
C
C DUMMY SUBROUTINE TO SIMULATE RANDOM ACCESS MASS STORE WRITE
C WRITE BLOCK, INDEX JB, FORM BUFA TO MASS STORE, NB REAL VALUES
C COMMON ARRAY RMAS SIMULATES MASS STORE ARRAY
C
      COMMON/MASS/RMAS(1024)
      REAL BUFA(NB)
C
      IOF=(JB-1)*NB
      DO 10 I=1,NB
      K=IOF+I
10    RMAS(K)=BUFA(I)
      RETURN
C
      END
C     PROGRAM MASTOM(TAPE1,INPUT,OUTPUT,TAPE60=INPUT,TAPE5=OUTPUT)
C
C  CONTROL DATA 6000 AND CYBER MASS STORE I/O FFT TEST PROGRAM
C
C  NOTE WELL THAT TYPE COMPLEX DATA MUST EXIST AS ALTERNATING
C  REAL/IMAG/REAL/IMAG... ELEMENTS, BOTH IN MASS STORE AND IN FORTRAN
C  WORKING ARRAY BUFA;  IN THIS FFT, DIFFERENT SUBROUTINES WILL SET
C  DIFFERENT TYPE (REAL OR COMPLEX) FOR ARRAY BUFA.
C
C
      COMMON/FFTCOM/LUN,MINDX(512)
C
C COMMON /FFTCOM/ HOLDS LOGICAL UNIT NUMBER FOR MASS STORE I/O
C ARRAY MINDX HOLDS RECORD INDICES FOR CYBER MASS STORE I/O
C
      COMMON/MFARG/MEXA(4),NDIM,ISGN,IDIR,SCAL,IBEX,ICEX,IPAK
      COMMON/MFVAL/DIMA(4),TDM1,RDM1,FBLK,TBLK,RBLK,RCOR,SIZE
      COMMON/MFINT/NDMA(4),NTD1,NRD1,NFBK,NTBK,NRBK,NRCR,NSZE
C
C COMMON AREAS /MFARG/,/MFVAL/,/MFINT/ USEFUL FOR RUNNING MASS STORE FFT
C  /MFARG/ HOLDS ARGUMENTS USED IN FFT CALLES, MOSTLY EXPONENTS
C    HELPER ROUTINE MFPAR COMPUTES VALUES FROM /MFARG/ INTO
C  /MFVAL/ (REALS AS SOME LARGE), /MFINT/ (INTEGER EQUIVALENTS IF POSS)
C    OR CAN REVERSE-COMPUTE SOME /MFARG/ FROM /MFVAL/
C       SEE COMMENTS, ROUTINE MFPAR.
C
      COMPLEX CBUFA(4096)
      REAL BUFA(8192)
      EQUIVALENCE (CBUFA(1),BUFA(1))
C
C  BUFA()=CBUFA()    IS WORKING AREA IN CORE STORE FOR FFT AND PROGRAM
C
      LUN=1
C LUN IS LOGICAL UNIT FOR CYBER MASS STORE I/O
C
      LP=5
      IPRINT=0
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF=-1
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN=-1
      IDIR=-1
      IPAK=1
C FFT ARGS, ISGN=+/- 1, IDIR=-1 (RMFFT), IDIR=1 OR 0 (CMFFT)
C IPAK=1, 0 OR -1 (RMFFT), NOT USED (CMFFT)
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      IF(IPRINT.LE.0)WRITE(LP,910)IPRINT,IRMF
910   FORMAT(31H1MASS STORE FFT TEST - IPRINT =,I3,
     X  25H (1=COPIOUS,0=MAX DIFFS),,
     X  9H   IRMF =,I3,28H (0=CMFFT TEST,1=RMFFT TEST) /)
      DIFMG=0.
C
      MEXA(1)=8
      MEXA(2)=6
      NDIM=2
      IBEX=9
      ICEX=13
C MORE FFT ARGS, 2**8 ROWS OF 2**6 ELMTS, 2 DIMEN,
C  FFT MASS STORE BLOCKS 2**IBEX REALS, BUFA 2**ICEX REALS
C
      IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      CALL OPENMS(LUN,MINDX,NFBK+1,0)
C CYBER 'READMS/WRITMS' MASS STORE OPENED WITH 'NUM REC'+1=NFBK+1
C  (NOTE HELPER ROUTN MFPAR RETURNS NFBK AS MAXIMUM FILE SIZE)
C
      NCNT=NRBK
C HELPER ROUTINE, NCNT NUM OF ELMTS IN FFT BLOCK, NTBK TOTAL BLOCKS
C
      M=MFSUM(MEXA,NDIM,99)
      SCAL=1.
      VALU=0.
      VINC=1./SIZE
      IF(IRMF.EQ.0)GO TO 50
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 40 JB=1,NTBK
      DO 35 I=1,NCNT
      BUFA(I)=VALU
35    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRBK,JB)
40    CONTINUE
C LOAD RAMP FUNCTN IN REAL ARRAY IN 'RECDS' OF NRBK LENGTH
C  (NOTE, CYBER MFREAD/MFWRIT (USING READMS/WRITMS) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
      GO TO 70
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FFT BLOCK
50    NCNT=NCNT/2
      DO 60 JB=1,NTBK
      DO 55 I=1,NCNT
      CBUFA(I)=CMPLX(VALU,0.)
55    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRBK,JB)
60    CONTINUE
C LOAD RAMP FUNCTN IN COMPLEX ARRAY IN 'RECDS' OF NRBK LENGTH
C  (NOTE, CYBER MFREAD/MFWRIT (USING READMS/WRITMS) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
70    IF(IPRINT.LE.0)GO TO 75
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
920   FORMAT(8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      WRITE(LP,930)
930   FORMAT(8H   INDEX,12X,3HFFT)
C
      IERR=MFPAR(-IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRBK/2
C HELPER ROUTN, NTBK TOTAL NUMB OF NRBK REALS IN FFT BLOCK
C  (NOTE IRMF=+1 FOR DATA IN COMPLEX STATE, NCNT ELMTS)
C
C BELOW, PROGRAM CAN ACCESS FFT COMPLEX RESULT (AND FORM COMPARISON, IF KNOWN)
      DO 80 JB=1,NTBK
      K=(JB-1)*NCNT
      CALL MFREAD(BUFA,NRBK,JB)
C READ 'RECDS' OF NRBK LENGTH, TOTALLING NTBK (RECOMPUTED BY MFPAR)
C  (NOTE, CYBER MFREAD/MFWRIT (USING READMS/WRITMS) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      DO 80 I=1,NCNT
      INDEX=J-1
      IF(IPRINT.LE.0)GO TO 80
      WRITE(LP,940)INDEX,CBUFA(I)
940   FORMAT(1X,I5,2X,2E13.4)
80    CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
75    ISGN=-ISGN
      IDIR=-IDIR
      SCAL=1./SIZE
      IF(IRMF.NE.0)
     X  CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF(IRMF.EQ.0)
     X  CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF(IPRINT.LE.0)GO TO 85
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
      WRITE(LP,960)
960   FORMAT(8H   INDEX,4X,8HFFT/2**M,6X,5HINPUT,10X,4HDIFF)
C
85    IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRBK
      IF(IRMF.EQ.0)NCNT=NCNT/2
C HELPER ROUTN, NTBK TOTAL NUMB OF NRBK REALS IN FFT BLOCK
C   (NCNT IS NUMBER OF ELMTS IN FFT BLOCK, REAL OR CMPLX)
C
C BELOW, COMPARES FFT INVERSE RESULT WITH INITIAL RANDOM INPUT
      DIFM=0.
      VALU=0.
      DO 90 JB=1,NTBK
      CALL MFREAD(BUFA,NRBK,JB)
C READ 'RECDS' OF NRBK LENGTH, TOTALLING NTBK (RECOMPUTED BY MFPAR)
C  (NOTE, CYBER MFREAD/MFWRIT (USING READMS/WRITMS) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      DO 90 I=1,NCNT
      IF(IRMF.NE.0)RM=BUFA(I)
      IF(IRMF.EQ.0)RM=REAL(CBUFA(I))
      DIF=ABS(RM-VALU)
      IF(DIF.GT.DIFM)DIFM=DIF
      IF(IPRINT.GT.0)WRITE(LP,990)I,RM,VALU,DIF
990   FORMAT(1X,I5,3(2X,E13.4))
      VALU=VALU+VINC
90    CONTINUE
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF(IPRINT.GE.0)WRITE(LP,950)DIFM,NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
950   FORMAT(10H MAX DIFF ,E11.4,1X,
     X  8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      IF(DIFM.GT.DIFMG)DIFMG=DIFM
C
      STOP
C
1000  WRITE(LP,980)IERR
980   FORMAT(25H IBEX FORCED TOO SMALL OR,
     X   30H NOT ALL /MFINT/ CORRECT, IERR,I4)
      STOP
      END
      SUBROUTINE MFREAD(BUFA,NB,JB)
C
C  CONTROL DATA 6000 AND CYBER MASS STORE I/O ROUTINES FOR FFT
C
C  (LOGICAL UNIT LUN IN COMMON/FFTCOM/ MUST HAVE BEEN OPENED
C  PREVIOUSLY; EG.      CALL OPENMS(LUN,INDXARRAY,NREC+1,0)
C  WHERE NREC=2**(M-IBEX+1) IF IPAK=1 OR LESS IF IPAK=0 OR -1)
C
C  SEE ALSO ALTERNATIVE SUBROUTINES USING EXTENDED CORE OR LCM
C  (IN ADDITION, GETW, PUTW OR READM, WRITEM MACROS CAN BE USED)
C
C READ BLOCK, INDEX JB, FROM MASS STORE TO BUFA, NB REAL VALUES
C
      COMMON/FFTCOM/LUN,MINDX(512)
      REAL BUFA(NB)
C
      CALL READMS(LUN,BUFA,NB,JB)
      RETURN
C
      ENTRY MFWRIT
C WRITE BLOCK, INDEX JB, FROM BUFA TO MASS STORE, NB REAL VALUES
C
      CALL WRITMS(LUN,BUFA,NB,JB,-1,0)
      RETURN
C
      END
C     PROGRAM LCMTOM(INPUT,OUTPUT,TAPE60=INPUT,TAPE5=OUTPUT)
C
C  CONTROL DATA 6000 AND CYBER EXTENDED CORE/LCM FFT TEST PROGRAM
C
C  NOTE WELL THAT TYPE COMPLEX DATA MUST EXIST AS ALTERNATING
C  REAL/IMAG/REAL/IMAG... ELEMENTS, BOTH IN MASS STORE AND IN FORTRAN
C  WORKING ARRAY BUFA;  IN THIS FFT, DIFFERENT SUBROUTINES WILL SET
C  DIFFERENT TYPE (REAL OR COMPLEX) FOR ARRAY BUFA.
C
C
      LEVEL 3, LBUFA
      COMMON/FFTCOM/LBUFA(32768)
C
C  LBUFA IN EXTENDED CORE/LCM SIMULATES FAST MASS STORE
C  DIMENSION LBUFA AS LARGE AS NECESSARY FOR RESULTANT ARRAY
C
      COMMON/MFARG/MEXA(4),NDIM,ISGN,IDIR,SCAL,IBEX,ICEX,IPAK
      COMMON/MFVAL/DIMA(4),TDM1,RDM1,FBLK,TBLK,RBLK,RCOR,SIZE
      COMMON/MFINT/NDMA(4),NTD1,NRD1,NFBK,NTBK,NRBK,NRCR,NSZE
C
C COMMON AREAS /MFARG/,/MFVAL/,/MFINT/ USEFUL FOR RUNNING MASS STORE FFT
C  /MFARG/ HOLDS ARGUMENTS USED IN FFT CALLES, MOSTLY EXPONENTS
C    HELPER ROUTINE MFPAR COMPUTES VALUES FROM /MFARG/ INTO
C  /MFVAL/ (REALS AS SOME LARGE), /MFINT/ (INTEGER EQUIVALENTS IF POSS)
C    OR CAN REVERSE-COMPUTE SOME /MFARG/ FROM /MFVAL/
C       SEE COMMENTS, ROUTINE MFPAR.
C
      COMPLEX CBUFA(4096)
      REAL BUFA(8192)
      EQUIVALENCE (CBUFA(1),BUFA(1))
C
C  BUFA()=CBUFA()    IS WORKING AREA IN CORE STORE FOR FFT AND PROGRAM
C
      LP=5
      IPRINT=0
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF=-1
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN=-1
      IDIR=-1
      IPAK=1
C FFT ARGS, ISGN=+/- 1, IDIR=-1 (RMFFT), IDIR=1 OR 0 (CMFFT)
C IPAK=1, 0 OR -1 (RMFFT), NOT USED (CMFFT)
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      IF(IPRINT.LE.0)WRITE(LP,910)IPRINT,IRMF
910   FORMAT(31H1MASS STORE FFT TEST - IPRINT =,I3,
     X  25H (1=COPIOUS,0=MAX DIFFS),,
     X  9H   IRMF =,I3,28H (0=CMFFT TEST,1=RMFFT TEST) /)
      DIFMG=0.
C
      MEXA(1)=8
      MEXA(2)=6
      NDIM=2
      IBEX=7
      ICEX=13
C MORE FFT ARGS, 2**8 ROWS OF 2**6 ELMTS, 2 DIMEN,
C  FFT MASS STORE BLOCKS 2**IBEX REALS, BUFA 2**ICEX REALS
C
      IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1
C HELPER ROUTINE, NCNT NUM OF ELMTS IN FIRST DIMEN, NTD1 TOTAL
C
      M=MFSUM(MEXA,NDIM,99)
      SCAL=1.
      VALU=0.
      VINC=1./SIZE
      IF(IRMF.EQ.0)GO TO 50
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 40 JB=1,NTD1
      DO 35 I=1,NCNT
      BUFA(I)=VALU
35    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRD1,JB)
40    CONTINUE
C LOAD RAMP FUNCTN IN REAL ARRAY IN 'RECDS' OF NRD1 LENGTH
C  (NOTE, 'LCM' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN CYBER MASS STORE)
C
      CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
      GO TO 70
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FFT BLOCK
50    NCNT=NCNT/2
      DO 60 JB=1,NTD1
      DO 55 I=1,NCNT
      CBUFA(I)=CMPLX(VALU,0.)
55    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRD1,JB)
60    CONTINUE
C LOAD RAMP FUNCTN IN COMPLEX ARRAY IN 'RECDS' OF NRD1 LENGTH
C  (NOTE, 'LCM' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN CYBER MASS STORE)
C
      CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
70    IF(IPRINT.LE.0)GO TO 75
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
920   FORMAT(8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      WRITE(LP,930)
930   FORMAT(8H   INDEX,12X,3HFFT)
C
      IERR=MFPAR(-IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1/2
C HELPER ROUTINE, NCNT NUM OF ELMTS IN FIRST DIMEN, NTD1 TOTAL
C  (NOTE IRMF=+1 FOR DATA IN COMPLEX STATE, NCNT ELMTS)
C
C BELOW, PROGRAM CAN ACCESS FFT COMPLEX RESULT (AND FORM COMPARISON, IF KNOWN)
      DO 80 JB=1,NTD1
      K=(JB-1)*NCNT
      CALL MFREAD(BUFA,NRD1,JB)
C READ 'RECDS' OF NRD1 LENGTH, TOTALLING NTD1 (RECOMPUTED BY MFPAR)
C  (NOTE, 'LCM' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN CYBER MASS STORE)
C
      DO 80 I=1,NCNT
      INDEX=J-1
      IF(IPRINT.LE.0)GO TO 80
      WRITE(LP,940)INDEX,CBUFA(I)
940   FORMAT(1X,I5,2X,2E13.4)
80    CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
75    ISGN=-ISGN
      IDIR=-IDIR
      SCAL=1./SIZE
      IF(IRMF.NE.0)
     X  CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF(IRMF.EQ.0)
     X  CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF(IPRINT.LE.0)GO TO 85
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
      WRITE(LP,960)
960   FORMAT(8H   INDEX,4X,8HFFT/2**M,6X,5HINPUT,10X,4HDIFF)
C
85    IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1
      IF(IRMF.EQ.0)NCNT=NCNT/2
C HELPER ROUTINE, NCNT NUM OF ELMTS IN FIRST DIMEN, NTD1 TOTAL
C   (NCNT IS NUMBER OF ELMTS IN FFT BLOCK, REAL OR CMPLX)
C
C BELOW, COMPARES FFT INVERSE RESULT WITH INITIAL RANDOM INPUT
      DIFM=0.
      VALU=0.
      DO 90 JB=1,NTD1
      CALL MFREAD(BUFA,NRD1,JB)
C READ 'RECDS' OF NRD1 LENGTH, TOTALLING NTD1 (RECOMPUTED BY MFPAR)
C  (NOTE, 'LCM' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN CYBER MASS STORE)
C
      DO 90 I=1,NCNT
      IF(IRMF.NE.0)RM=BUFA(I)
      IF(IRMF.EQ.0)RM=REAL(CBUFA(I))
      DIF=ABS(RM-VALU)
      IF(DIF.GT.DIFM)DIFM=DIF
      IF(IPRINT.GT.0)WRITE(LP,990)I,RM,VALU,DIF
990   FORMAT(1X,I5,3(2X,E13.4))
      VALU=VALU+VINC
90    CONTINUE
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF(IPRINT.GE.0)WRITE(LP,950)DIFM,NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
950   FORMAT(10H MAX DIFF ,E11.4,1X,
     X  8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      IF(DIFM.GT.DIFMG)DIFMG=DIFM
C
      STOP
C
1000  WRITE(LP,980)IERR
980   FORMAT(25H IBEX FORCED TOO SMALL OR,
     X   30H NOT ALL /MFINT/ CORRECT, IERR,I4)
      STOP
      END
      SUBROUTINE MFREAD(BUFA,NB,JB)
C
C  CONTROL DATA 6000 AND CYBER EXTENDED CORE STORE I/O ROUTINES FOR FFT
C
C  (EXTENDED CORE OR LCM TAKES PLACE OF MASS STORE -
C   DIMENSION LBUFA AS LARGE AS NECESSARY FOR LARGE ARRAYS)
C
C  ALTERNATIVELY, EXTENDED CORE USE COULD BE COMBINED WITH
C  READM/WRITEM MACROS TO ACCESS LARGER FILE WHEN NECESSARY
C  (USE VIRTUAL MEMORY ALGORITHM - MANY SECTORS HELD IN
C    EXTENDED CORE AND ONLY ACCESSED FROM DISC IF NOT PRESENT).
C
C READ BLOCK, INDEX JB, FROM EXTENDED CORE TO BUFA, NB REAL VALUES
C
      LEVEL 3, LBUFA
      COMMON/FFTCOM/LBUFA(32768)
      REAL BUFA(NB)
C
      CALL MOVLEV(LBUFA((JB-1)*NB+1),BUFA,NB)
      RETURN
C
      ENTRY MFWRIT
C WRITE BLOCK, INDEX JB, FROM BUFA TO EXTENDED CORE, NB REAL VALUES
C
      CALL MOVLEV(BUFA,LBUFA((JB-1)*NB+1),NB)
      RETURN
C
      END
C
C PDP 11 MASS STORE I/O FFT SAMPLE TEST PROGRAM
C
C  NOTE WELL THAT TYPE COMPLEX DATA MUST EXIST AS ALTERNATING
C  REAL/IMAG/REAL/IMAG... ELEMENTS, BOTH IN MASS STORE AND IN FORTRAN
C  WORKING ARRAY BUFA;  IN THIS FFT, DIFFERENT SUBROUTINES WILL SET
C  DIFFERENT TYPE (REAL OR COMPLEX) FOR ARRAY BUFA.
C
C
      COMMON/FFTCOM/LUN
C
C COMMON /FFTCOM/ HOLDS LOGICAL UNIT NUMBER FOR MASS STORE I/O
C
      COMMON/MFARG/MEXA(4),NDIM,ISGN,IDIR,SCAL,IBEX,ICEX,IPAK
      COMMON/MFVAL/DIMA(4),TDM1,RDM1,FBLK,TBLK,RBLK,RCOR,SIZE
      COMMON/MFINT/NDMA(4),NTD1,NRD1,NFBK,NTBK,NRBK,NRCR,NSZE
C
C COMMON AREAS /MFARG/,/MFVAL/,/MFINT/ USEFUL FOR RUNNING MASS STORE FFT
C  /MFARG/ HOLDS ARGUMENTS USED IN FFT CALLES, MOSTLY EXPONENTS
C    HELPER ROUTINE MFPAR COMPUTES VALUES FROM /MFARG/ INTO
C  /MFVAL/ (REALS AS SOME LARGE), /MFINT/ (INTEGER EQUIVALENTS IF POSS)
C    OR CAN REVERSE-COMPUTE SOME /MFARG/ FROM /MFVAL/
C       SEE COMMENTS, ROUTINE MFPAR.
C
      COMPLEX CBUFA(1024)
      REAL BUFA(2048)
      EQUIVALENCE (CBUFA(1),BUFA(1))
C
C  BUFA()=CBUFA()    IS WORKING AREA IN CORE STORE FOR FFT AND PROGRAM
C
      LUN=1
C LUN IS LOGICAL UNIT FOR PDP 11 MASS STORE I/O
C
      LP=5
      IPRINT=0
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF=-1
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN=-1
      IDIR=-1
      IPAK=1
C FFT ARGS, ISGN=+/- 1, IDIR=-1 (RMFFT), IDIR=1 OR 0 (CMFFT)
C IPAK=1, 0 OR -1 (RMFFT), NOT USED (CMFFT)
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      IF(IPRINT.LE.0)WRITE(LP,910)IPRINT,IRMF
910   FORMAT(31H1MASS STORE FFT TEST - IPRINT =,I3,
     X  25H (1=COPIOUS,0=MAX DIFFS),,
     X  9H   IRMF =,I3,28H (0=CMFFT TEST,1=RMFFT TEST) /)
      DIFMG=0.
C
      MEXA(1)=8
      MEXA(2)=6
      NDIM=2
      IBEX=7
      ICEX=11
C MORE FFT ARGS, 2**8 ROWS OF 2**6 ELMTS, 2 DIMEN,
C  FFT MASS STORE BLOCKS 2**IBEX REALS, BUFA 2**ICEX REALS
C
      IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      CALL ASSIGN(LUN,'FFTEST.DAT',0)
      NWD=NRBK*2
      DEFINE FILE LUN(NFBK,NWD,U,INDX)
C PDP 11 MASS STORE OPENED WITH NUM REC=NFBK, NRBK*2 WORDS PER REC
C  (NOTE HELPER ROUTN MFPAR RETURNS NFBK AS MAXIMUM FILE SIZE)
C
      NCNT=NRBK
C HELPER ROUTINE, NCNT NUM OF ELMTS IN FFT BLOCK, NTBK TOTAL BLOCKS
C
      M=MFSUM(MEXA,NDIM,99)
      SCAL=1.
      VALU=0.
      VINC=1./SIZE
      IF(IRMF.EQ.0)GO TO 50
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 40 JB=1,NTBK
      DO 35 I=1,NCNT
      BUFA(I)=VALU
35    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRBK,JB)
40    CONTINUE
C LOAD RAMP FUNCTN IN REAL ARRAY IN 'RECDS' OF NRBK LENGTH
C  (NOTE, PDP 11 MFREAD/MFWRIT (USING FORTRAN I/O) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
      GO TO 70
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FFT BLOCK
50    NCNT=NCNT/2
      DO 60 JB=1,NTBK
      DO 55 I=1,NCNT
      CBUFA(I)=CMPLX(VALU,0.)
55    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRBK,JB)
60    CONTINUE
C LOAD RAMP FUNCTN IN COMPLEX ARRAY IN 'RECDS' OF NRBK LENGTH
C  (NOTE, PDP 11 MFREAD/MFWRIT (USING FORTRAN I/O) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
70    IF(IPRINT.LE.0)GO TO 75
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
920   FORMAT(8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      WRITE(LP,930)
930   FORMAT(8H   INDEX,12X,3HFFT)
C
      IERR=MFPAR(-IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRBK/2
C HELPER ROUTN, NTBK TOTAL NUMB OF NRBK REALS IN FFT BLOCK
C  (NOTE IRMF=+1 FOR DATA IN COMPLEX STATE, NCNT ELMTS)
C
C BELOW, PROGRAM CAN ACCESS FFT COMPLEX RESULT (AND FORM COMPARISON, IF KNOWN)
      DO 80 JB=1,NTBK
      K=(JB-1)*NCNT
      CALL MFREAD(BUFA,NRBK,JB)
C READ 'RECDS' OF NRBK LENGTH, TOTALLING NTBK (RECOMPUTED BY MFPAR)
C  (NOTE, PDP 11 MFREAD/MFWRIT (USING FORTRAN I/O) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      DO 80 I=1,NCNT
      INDEX=J-1
      IF(IPRINT.LE.0)GO TO 80
      WRITE(LP,940)INDEX,CBUFA(I)
940   FORMAT(1X,I5,2X,2E13.4)
80    CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
75    ISGN=-ISGN
      IDIR=-IDIR
      SCAL=1./SIZE
      IF(IRMF.NE.0)
     X  CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF(IRMF.EQ.0)
     X  CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF(IPRINT.LE.0)GO TO 85
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
      WRITE(LP,960)
960   FORMAT(8H   INDEX,4X,8HFFT/2**M,6X,5HINPUT,10X,4HDIFF)
C
85    IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRBK
      IF(IRMF.EQ.0)NCNT=NCNT/2
C HELPER ROUTN, NTBK TOTAL NUMB OF NRBK REALS IN FFT BLOCK
C   (NCNT IS NUMBER OF ELMTS IN FFT BLOCK, REAL OR CMPLX)
C
C BELOW, COMPARES FFT INVERSE RESULT WITH INITIAL RANDOM INPUT
      DIFM=0.
      VALU=0.
      DO 90 JB=1,NTBK
      CALL MFREAD(BUFA,NRBK,JB)
C READ 'RECDS' OF NRBK LENGTH, TOTALLING NTBK (RECOMPUTED BY MFPAR)
C  (NOTE, PDP 11 MFREAD/MFWRIT (USING FORTRAN I/O) CANNOT ACCEPT 'RECDS'
C   OF DIFFERENT LENGTH; OTHERWISE COULD USE NTD1 'RECDS' OF NRD1 HERE)
C
      DO 90 I=1,NCNT
      IF(IRMF.NE.0)RM=BUFA(I)
      IF(IRMF.EQ.0)RM=REAL(CBUFA(I))
      DIF=ABS(RM-VALU)
      IF(DIF.GT.DIFM)DIFM=DIF
      IF(IPRINT.GT.0)WRITE(LP,990)I,RM,VALU,DIF
990   FORMAT(1X,I5,3(2X,E13.4))
      VALU=VALU+VINC
90    CONTINUE
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF(IPRINT.GE.0)WRITE(LP,950)DIFM,NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
950   FORMAT(10H MAX DIFF ,E11.4,1X,
     X  8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      IF(DIFM.GT.DIFMG)DIFMG=DIFM
C
      STOP
C
1000  WRITE(LP,980)IERR
980   FORMAT(25H IBEX FORCED TOO SMALL OR,
     X   30H NOT ALL /MFINT/ CORRECT, IERR,I4)
      STOP
      END
      SUBROUTINE MFREAD(BUFA,NB,JB)
C
C  PDP 11 FORTRAN DIRECT ACCESS READ ROUTINE FOR FFT
C
C  (LOGICAL UNIT LUN IN COMMON/FFTCOM/ MUST HAVE BEEN OPENED
C  PREVIOUSLY; EG.      CALL ASSIGN(LUN,'FILENAME',0)
C  AND                  DEFINE FILE LUN(NREC,NWD,U,INDX)
C  WHERE NWD=2**(IBEX+1) AND NREC=2**(M-IBEX+1) OR LESS IF IPAK=0 OR -1)
C
C  SEE ALSO ALTERNATIVE, FAST MACRO I/O SUBROUTINES
C  (THESE ARE TO BE PREFERRED, SINCE SPEED 2 TO 10 TIMES BETTER)
C
C READ BLOCK, INDEX JB, FROM MASS STORE TO BUFA, NB REAL VALUES
C
      COMMON/FFTCOM/LUN
      REAL BUFA(NB)
C
      READ(LUN'JB)BUFA
      RETURN
C
      END
      SUBROUTINE MFWRIT(BUFA,NB,JB)
C
C  PDP 11 FORTRAN DIRECT ACCESS WRITE ROUTINE FOR FFT
C
C  (LOGICAL UNIT LUN IN COMMON/FFTCOM/ MUST HAVE BEEN OPENED
C  PREVIOUSLY; EG.      CALL ASSIGN(LUN,'FILENAME',0)
C  AND                  DEFINE FILE LUN(NREC,NWD,U,INDX)
C  WHERE NWD=2**(IBEX+1) AND NREC=2**(M-IBEX+1) OR LESS IF IPAK=0 OR -1)
C
C  SEE ALSO ALTERNATIVE, FAST MACRO I/O SUBROUTINES
C  (THESE ARE TO BE PREFERRED, SINCE SPEED 2 TO 10 TIMES BETTER)
C
C WRITE BLOCK, INDEX JB, FROM BUFA TO MASS STORE, NB REAL VALUES
C
      COMMON/FFTCOM/LUN
      REAL BUFA(NB)
C
      WRITE(LUN'JB)BUFA
      RETURN
C
      END
C
C PDP 11 FAST MACRO I/O FFT SAMPLE TEST PROGRAM
C (REQUIRES CALL TO SYSTEM MACRO TO OPEN FILE FOR I/O - SEE LINE 70)
C
C  NOTE THAT FAST MACRO I/O IS MORE EFFEICIENT THAN STANDARD FORTRAN I/O
C
C  NOTE WELL THAT TYPE COMPLEX DATA MUST EXIST AS ALTERNATING
C  REAL/IMAG/REAL/IMAG... ELEMENTS, BOTH IN MASS STORE AND IN FORTRAN
C  WORKING ARRAY BUFA;  IN THIS FFT, DIFFERENT SUBROUTINES WILL SET
C  DIFFERENT TYPE (REAL OR COMPLEX) FOR ARRAY BUFA.
C
C
      COMMON/FFTCOM/IFDB,IOERR
C
C COMMON /FFTCOM/ HOLDS FDB ADDRESS FOR MACRO I/O (RSX11M), CHANNEL (RT11)
C     (SEE COMMENTS IN FAST MACRO ROUTINE MFREAD/MFWRIT)
C
      COMMON/MFARG/MEXA(4),NDIM,ISGN,IDIR,SCAL,IBEX,ICEX,IPAK
      COMMON/MFVAL/DIMA(4),TDM1,RDM1,FBLK,TBLK,RBLK,RCOR,SIZE
      COMMON/MFINT/NDMA(4),NTD1,NRD1,NFBK,NTBK,NRBK,NRCR,NSZE
C
C COMMON AREAS /MFARG/,/MFVAL/,/MFINT/ USEFUL FOR RUNNING MASS STORE FFT
C  /MFARG/ HOLDS ARGUMENTS USED IN FFT CALLES, MOSTLY EXPONENTS
C    HELPER ROUTINE MFPAR COMPUTES VALUES FROM /MFARG/ INTO
C  /MFVAL/ (REALS AS SOME LARGE), /MFINT/ (INTEGER EQUIVALENTS IF POSS)
C    OR CAN REVERSE-COMPUTE SOME /MFARG/ FROM /MFVAL/
C       SEE COMMENTS, ROUTINE MFPAR.
C
      COMPLEX CBUFA(1024)
      REAL BUFA(2048)
      EQUIVALENCE (CBUFA(1),BUFA(1))
C
C  BUFA()=CBUFA()    IS WORKING AREA IN CORE STORE FOR FFT AND PROGRAM
C
      LUN=1
C LUN IS LOGICAL UNIT FOR PDP 11 MASS STORE I/O
C
      LP=5
      IPRINT=0
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF=-1
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN=-1
      IDIR=-1
      IPAK=1
C FFT ARGS, ISGN=+/- 1, IDIR=-1 (RMFFT), IDIR=1 OR 0 (CMFFT)
C IPAK=1, 0 OR -1 (RMFFT), NOT USED (CMFFT)
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      IF(IPRINT.LE.0)WRITE(LP,910)IPRINT,IRMF
910   FORMAT(31H1MASS STORE FFT TEST - IPRINT =,I3,
     X  25H (1=COPIOUS,0=MAX DIFFS),,
     X  9H   IRMF =,I3,28H (0=CMFFT TEST,1=RMFFT TEST) /)
      DIFMG=0.
C
      MEXA(1)=8
      MEXA(2)=6
      NDIM=2
      IBEX=7
      ICEX=11
C MORE FFT ARGS, 2**8 ROWS OF 2**6 ELMTS, 2 DIMEN,
C  FFT MASS STORE BLOCKS 2**IBEX REALS, BUFA 2**ICEX REALS
C
      IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
C  (NOTE HELPER ROUTN MFPAR RETURNS NFBK AS MAXIMUM FILE SIZE)
C
      STOP 'ASSIGN AND OPEN FILE HERE'
C DELETE ABOVE LINE AND INVOKE SYSTEM MACRO 'OPEN$' (RSX) OR '.OPEN (RT11)
C
C FOR EXAMPLE, A FORTRAN CALLABLE SUBROUTINE 'MFOPEN' COULD OPEN FILE
C  'FFTTEST.DAT', RETURNING IFDB=FDB ADDRESS (RSX11M) IN COMMON/FFTCOM/
C  OF CHANNEL NUMBER (RT11), FOR USE BY MACRO MFREAD/MFWRIT, THUS:
C
C     CALL MFOPEN(LUN,'FFTTEST.DAT',IFDB)
C
      IOERR=0
C PRESET IOERR IN /FFTCOM/; ZERO IF NO ERRORS IN I/O
C
      NCNT=NRD1
C HELPER ROUTINE, NCNT NUM OF ELMTS IN 'FIRST' DIMEN, NTD1 TOTAL BLOCKS
C
      M=MFSUM(MEXA,NDIM,99)
      SCAL=1.
      VALU=0.
      VINC=1./SIZE
      IF(IRMF.EQ.0)GO TO 50
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 40 JB=1,NTD1
      DO 35 I=1,NCNT
      BUFA(I)=VALU
35    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRD1,JB)
40    CONTINUE
C LOAD RAMP FUNCTN IN REAL ARRAY IN 'RECDS' OF NRD1 LENGTH
C  (NOTE, 'MACRO' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN FORTRAN MASS STORE)
C
      CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
      GO TO 70
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN 'FIRST' DIMEN
50    NCNT=NCNT/2
      DO 60 JB=1,NTD1
      DO 55 I=1,NCNT
      CBUFA(I)=CMPLX(VALU,0.)
55    VALU=VALU+VINC
      CALL MFWRIT(BUFA,NRD1,JB)
60    CONTINUE
C LOAD RAMP FUNCTN IN COMPLEX ARRAY IN 'RECDS' OF NRD1 LENGTH
C  (NOTE, 'MACRO' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN FORTRAN MASS STORE)
C
      CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
70    IF(IPRINT.LE.0)GO TO 75
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
920   FORMAT(8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      WRITE(LP,930)
930   FORMAT(8H   INDEX,12X,3HFFT)
C
      IERR=MFPAR(-IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1/2
C HELPER ROUTN, NTD1 TOTAL NUMB OF NRD1 REALS IN 'FIRST' DIMEN
C  (NOTE IRMF=+1 FOR DATA IN COMPLEX STATE, NCNT ELMTS)
C
C BELOW, PROGRAM CAN ACCESS FFT COMPLEX RESULT (AND FORM COMPARISON, IF KNOWN)
      DO 80 JB=1,NTD1
      K=(JB-1)*NCNT
      CALL MFREAD(BUFA,NRD1,JB)
C READ 'RECDS' OF NRD1 LENGTH, TOTALLING NTD1 (RECOMPUTED BY MFPAR)
C  (NOTE, 'MACRO' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN FORTRAN MASS STORE)
C
      DO 80 I=1,NCNT
      INDEX=J-1
      IF(IPRINT.LE.0)GO TO 80
      WRITE(LP,940)INDEX,CBUFA(I)
940   FORMAT(1X,I5,2X,2E13.4)
80    CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
75    ISGN=-ISGN
      IDIR=-IDIR
      SCAL=1./SIZE
      IF(IRMF.NE.0)
     X  CALL RMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX, IPAK)
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF(IRMF.EQ.0)
     X  CALL CMFFT (MEXA,NDIM,ISGN,IDIR,SCAL, BUFA,IBEX,ICEX)
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF(IPRINT.LE.0)GO TO 85
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
      WRITE(LP,910)IPRINT,IRMF
      WRITE(LP,920)NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
      WRITE(LP,960)
960   FORMAT(8H   INDEX,4X,8HFFT/2**M,6X,5HINPUT,10X,4HDIFF)
C
85    IERR=MFPAR(IRMF,0)
      IF(IERR.NE.0)GO TO 1000
      NCNT=NRD1
      IF(IRMF.EQ.0)NCNT=NCNT/2
C HELPER ROUTN, NTD1 TOTAL NUMB OF NRD1 REALS IN 'FIRST' DIMEN
C   (NCNT IS NUMBER OF ELMTS IN 'FIRST' DIMEN, REAL OR CMPLX)
C
C BELOW, COMPARES FFT INVERSE RESULT WITH INITIAL RANDOM INPUT
      DIFM=0.
      VALU=0.
      DO 90 JB=1,NTD1
      CALL MFREAD(BUFA,NRD1,JB)
C READ 'RECDS' OF NRD1 LENGTH, TOTALLING NTD1 (RECOMPUTED BY MFPAR)
C  (NOTE, 'MACRO' MFREAD/MFWRIT CAN ACCEPT 'RECORDS' OF DIFFERENT LENGTH;
C   SO CAN USE NTD1,NRD1 HERE INSTEAD OF NTBK,NRBK AS IN FORTRAN MASS STORE)
C
      DO 90 I=1,NCNT
      IF(IRMF.NE.0)RM=BUFA(I)
      IF(IRMF.EQ.0)RM=REAL(CBUFA(I))
      DIF=ABS(RM-VALU)
      IF(DIF.GT.DIFM)DIFM=DIF
      IF(IPRINT.GT.0)WRITE(LP,990)I,RM,VALU,DIF
990   FORMAT(1X,I5,3(2X,E13.4))
      VALU=VALU+VINC
90    CONTINUE
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF(IPRINT.GE.0)WRITE(LP,950)DIFM,NDIM,ISGN,IDIR,IBEX,ICEX,IPAK,
     X  (MEXA(J),J=1,NDIM)
950   FORMAT(10H MAX DIFF ,E11.4,1X,
     X  8H  NDIM =,I3,8H  ISGN =,I3,8H  IDIR =,I3,
     X  8H  IBEX =,I3,8H  ICEX =,I3,8H  IPAK =,I3,
     X  10H  MEXA() =,3I5)
      IF(DIFM.GT.DIFMG)DIFMG=DIFM
C
      STOP
C
1000  WRITE(LP,980)IERR
980   FORMAT(25H IBEX FORCED TOO SMALL OR,
     X   30H NOT ALL /MFINT/ CORRECT, IERR,I4)
      STOP
      END
;
;  PDP11 RSX11M OR RT11 FAST MACRO I/O ROUTINES FOR FFT
;
;     SUBROUTINE MFREAD(BUFA,NB,JB)
;     SUBROUTINE MFWRIT(BUFA,NB,JB)
;
;C READ  BLOCK, INDEX JB, FROM MASS STORE TO BUFA, NB REAL VALUES
;C WRITE BLOCK, INDEX JB, FROM BUFA TO MASS STORE, NB REAL VALUES
;
;     COMMON/FFTCOM/IFDB,IOERR  (RSX11 FDB ADDRESS, ERROR)
;OR   COMMON/FFTCOM/ICHAN,IOERR (RT11  CHANNEL NUM, ERROR)
;
;NOTE1: FILE MUST BE PROPERLY OPEN FOR READ/WRITE ACCESS AND WITH
;       RSX11M FDB ADDRESS OR RT11 CHANNEL NUMBER IN FIRST WORD
;       OF COMMON/FFTCOM/ (IE. INVOKE RSX11M OR RT11 OPEN MACROS).
;
;NOTE2: I/O ERRORS, IF THEY OCCUR, ARE RETURNED IN SECOND WORD
;       OF COMMON/FFTCOM/.  THIS WORD (IOERR) REMAINS UNCHANGED
;       IF NO ERROR OCCURS - THUS, IF PRESET TO 0 BEFORE CALLING
;       MFREAD/MFWRIT (OR MASS STORE FFT), WILL FINALLY BE
;       0 IF NO ERRORS, OR NON ZERO=LAST ERROR (SEE FOLLOWING LABEL L4:).
;
;NOTE3: NB MUST BE AN INTEGRAL POWER OF 2; 128 OR MORE MOST EFFICIENT
;       (NB.LE.64 USES LOCAL 256 WORD JBUF TO HOLD SECTOR;
;        WRITE ALWAYS WRITES THIS SECTOR TO FILE; SLOW BUT SAFE).
;
;NOTE4: READY TO ASSEMBLE FOR RSX11M. FOR RT11, ERASE FOLLOWING LINE:
        RSX=0
;
        .IF     DF,RSX
; CONDITIONAL SECTION FOR RSX11M (V3)
        .TITLE  MFREAD(RSX11M)
        .MCALL  READ$,WRITE$,WAIT$
        .PSECT  FFTCOM,RW,D,GBL,REL,OVR
IFDB:   .WORD   0       ;FDB ADDRESS
IOERR:  .WORD   0       ;ERROR FLAG (UNCHANGED IF NO ERROR)
        .PSECT
        .ENDC
;
        .IF     NDF,RSX
; CONDITIONAL SECTION FOR RT11 (V3)
        .TITLE  MFREAD(RT11)
        .MCALL  .READW,.WRITW
        ERRBYT=52
        .PSECT  FFTCOM,RW,D,GBL,REL,OVR
ICHAN:  .WORD   0       ;CHANNEL NUMBER
IOERR:  .WORD   0       ;ERROR FLAG (UNCHANGED IF NO ERROR)
        .PSECT  MFREAD,RW,I,LCL,REL,OVR
        .ENDC
;
        .GLOBL  MFREAD,MFWRIT
;
MFREAD:  MOV     #1,R0           ;MFREAD SETS JSW=1
        BR      L1
;
MFWRIT:  MOV     #-1,R0          ;MFWRIT SETS JSW=-1
;
L1:     MOV     R0,JSW          ;HOLD JSW
        TST     (R5)+           ;IGNORE ARGUMENT COUNT
        MOV     (R5)+,R3        ;R3=ADDRESS OF 'BUFA'
        MOV     @(R5)+,R2       ;R2='NB'
        MOV     @(R5)+,R1       ;R1='JB'
; FINISH ACCESSING SUBROUTINE ARGUMENTS
;
        DEC     R1              ;R1=JB-1
        ASL     R2              ;R2=WORD COUNT, NB*2
        MOV     R2,R5           ;COPY TO R5
        BEQ     L6              ;FINISH IF COUNT ZERO
        CLR     R0              ;CLEAR R0 FOR WORD OFFSET
        CMP     R5,#256.        ;CHECK WHETHER WHOLE SECTORS
        BLT     L12             ;GO TO L12 IF NOT (MORE COMPLEX)
;
; BELOW, SCALE R1 TO SECTOR OFFSET (ASSUMES NWD POWER OF 2)
L10:    BIT     #256.,R5        ;LOOK FOR BIT AT 256.
        BNE     L11             ;R1 COMPLETE WHEN FOUND
        ASL     R1              ;SCALE R1 ACCORDING TO R5
        ASR     R5
        BR      L10
;
; NOW HAVE      R1=SECTOR START (FIRST=0)
;               R2=WORD COUNT
;               R3=BUFA ADDRESS
;
L11:    MOV     JSW,JSWIO       ;SET UP FOR READ/WRITE
        BGT     L13
        MOV     #-1,SECHLD      ;DIRECT WRITE, ENSURE JBUF LOOKS EMPTY
L13:    JSR     PC,INOUT        ;CALL DIRECT INPUT/OUTPUT
L6:     RTS     PC              ;RETURN FROM SUBROUTINE
; END OF SIMPLE DIRECT INPUT/OUTPUT
;
        .IF     DF,RSX
; CONDITIONAL SECTION FOR RSX11M
INOUT:  ASL     R2              ;RSX11M NEEDS BYTE COUNT=NB*4
        INC     R1              ;RSX11M SECTOR STARTS WITH 1
        MOV     R1,SECTL        ;HOLD AS LOWER SECTOR VALUE
        MOV     IFDB,R4         ;R4=FDB ADDRESS
        MOVB    F.LUN(R4),R1    ;USE LUN FOR EVENT FLAG
;
        TST     JSWIO
        BLT     L5              ;BRANCH IF WRITE
;
        MOV     F.EFBK+2(R4),R0 ;GET LOW SECTOR OF EOF
        DEC     R0              ;ALLOW FOR HEADER SECTOR
        CMP     R0,SECTL        ;TEST IF SECTOR EXISTS YET
        BLT     L7              ;NO READ IF DOES NOT EXIST
;
        READ$   R4,R3,R2,#SECT,R1 ;FDB,BUF,BYTCNT,SECT,EVFLG
;
L4:     MOV     IFDB,R4         ;R4=FDB ADDRESS
        MOVB    F.LUN(R4),R1    ;USE LUN FOR EVENT FLAG
        WAIT$   R4,R1           ;WAIT FOR I/O FINISH
        BCC     L7              ;NO ERROR IF CARRY CLEAR
;
        MOV     IFDB,R4         ;R4=FDB ADDRESS
        MOVB    F.ERR(R4),R4    ;HOLD NEGATIVE ERROR CODE
        MOV     R4,IOERR        ;IN IOERR IN COMMON/FFTCOM/
L7:     RTS     PC
;
L5:     WRITE$  R4,R3,R2,#SECT,R1 ;FDB,BUF,BYTCNT,SECT,EVFLG
        BR      L4
;
SECT:   .WORD   0       ;HOLDS SECTOR VALUE (HIGHER, ALWAYS 0)
SECTL:  .WORD   0       ;(LOWER, COMPUTED)
        .ENDC
;
        .IF     NDF,RSX
;CONDITIONAL SECTION FOR RT11
INOUT:  MOV     ICHAN,R4        ;R4=CHANNEL NUMBER
;
        TST     JSWIO
        BLT     L5              ;BRANCH IF WRITE
;
        .READW  #AREA,R4,R3,R2,R1       ;AREA,CHAN,BUF,WDCNT,SECT
;
L4:     BCC     L7              ;NO ERROR IF CARRY CLEAR
;
        MOVB    ERRBYT,R4       ;HOLD ERROR CODE + 1 (TO MAKE NON ZERO)
        INC     R4               ;ERRORS 0,1,2 GO TO 1,2,3
        MOV     R4,IOERR        ;IN IOERR IN COMMON/FFTCOM/
L7:     RTS     PC
;
L5:     .WRITW  #AREA,R4,R3,R2,R1       ;AREA,CHAN,BUF,WDCNT,SECT
        BR      L4
AREA:   .BLKW   10      ;AREA FOR RT11 I/O MACROS USE
        .ENDC
;
; BELOW, MORE COMPLEX HANDLING OF NWD.LT.256 (PART SECTORS)
L12:    CLC
        ROR     R1              ;SCALE R1/R0 RIGHT ACCORDING TO R5
        ROR     R0              ;R0 HOLDS FLOW OUT OF R1
        ASL     R5
        BIT     #256.,R5        ;LOOK FOR BIT 256
        BNE     L20             ;R1/R0 COMPLETE IF FOUND
        BR      L12
;
L20:    SWAB    R0
        ASL     R0
;
; NOW HAVE      R0=BYTE OFFSET IN LOCAL BUFFER
;               R1=SECTOR REQUIRED
;               R2=WORD COUNT
;               R3=USER BUFA ADDRESS
;
; BELOW, TRANSFER BETWEEN LOCAL AND USER BUFFER
        ADD     #JBUF,R0        ;R0=JBUF ADDRESS + OFFSET
        CMP     SECHLD,R1               ;CHECK CURRENT SECTOR LOADED
        BNE     L24             ;MUST FIRST LOAD SECTO AT L24 IF NEC
L21:    TST     JSW             ;TEST WHETHER READ/WRITE
        BLT     L23
;
; BELOW, 'READ' PART SECTOR
L22:    MOV     (R0)+,(R3)+     ;COPY LOCAL TO USER BUFA
        DEC     R2
        BGT     L22
        JMP     L6              ;FINISH
;
; BELOW, 'WRITE' PART SECTOR
L23:    MOV     (R3)+,(R0)+     ;COPY USER TO LOCAL BUF
        DEC     R2
        BGT     L23
        MOV     #-1,JSWIO               ;SET UP I/O FOR WRITE
        JSR     PC,LINOUT       ;WRITE FROM LOCAL JBUF (SLOW BUT SAFE)
        JMP     L6              ;FINISH
;
; BELOW, LOAD LOCAL JBUF, IF POSSIBLE
L24:    MOV     R1,SECHLD       ;HOLD CURRENT SECTOR NUM
        MOV     #1,JSWIO        ;SET UP FOR READ
        JSR     PC,LINOUT       ;DO LOCAL READ TO JBUF
        BR      L21
;
; BELOW, LOCAL INPUT/OUTPUT ROUTINE
LINOUT: MOV     R0,-(SP)        ;SAVE REGISTERS R0 TO R3
        MOV     R1,-(SP)
        MOV     R2,-(SP)
        MOV     R3,-(SP)
        MOV     #256.,R2        ;R2=256. FOR WHOLE SECTOR
        MOV     #JBUF,R3        ;R3 IS LOCAL JBUF ADDRESS
        JSR     PC,INOUT        ;CALL INOUT ROUTINE
        MOV     (SP)+,R3        ;RESTORE REGISTERS R0 TO R3
        MOV     (SP)+,R2
        MOV     (SP)+,R1
        MOV     (SP)+,R0
        RTS     PC
;
; BELOW, LOCAL VARIABLES
JSW:    .WORD   0       ;JSW=+1 MFREAD, -1 MFWRIT
JSWIO:  .WORD   0       ;ACTUAL JSW USED FOR INOUT ROUTINE
SECHLD: .WORD   -1      ;CURRENT SECTOR LOADED IN JBUF
JBUF:   .BLKW   256.    ;JBUF 256 WORD LOCAL BUFFER FOR NWD.LT.256
;
        .END