Make your own free website on Tripod.com
Optimized Mass Storage FFT Program

Return to Main Software Page

C
C-----------------------------------------------------------------------
C  PROGRAM:  AN OPTIMIZED MASS STORAGE FFT
C  AUTHOR:   DONALD FRASER
C            CSIRO, DIVISION OF COMPUTING RESEARCH
C            PO BOX 1800
C            CANBERRA CITY, ACT 2601, AUSTRALIA
C
C           REVISION DATE: JULY 1978.
C
C     THE SET INCLUDES A UNIVERSAL TEST 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 MAY 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 RESULT)
C  2  CMFFT        CALLED BY 1, OR MASS STORAGE FFT (ALL COMPLEX)
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-----------------------------------------------------------------------
C
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  RMFFT
C REAL-TO-COMPLEX FFT (OR VICE-VERSA) OF MULTI-DIMENSD MASS STORE ARRAY
C (FRASER, ACM TOMS - 1979, AND J.ACM, V.23,N.2, APRIL 76, PP. 298-309)
C MASS STORE ARRAY IS EITHER REAL OR COMPLEX DATA (SEE NOTE BELOW)
C-----------------------------------------------------------------------
C
      SUBROUTINE RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA, IBEX, ICEX,
     *    IPAK)
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 FIRST
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 TRANS.
C     THUS, 2**M PACKED REAL VALUES,
C     OR,   2**M COMPLEX VALUES, IF COMPLEX RESULT WITH IPAK=1
C RMFFT ALWAYS REVERSES DIMENSION ORDER AND MEXA LIST
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 CMFFT)
C  IPAK=0  COMPUTES COMPLEX ARRAY OF JUST OVER HALF SIZE
C  IPAK=-1 HOLDS COMPLEX ARRAY AT EXACTLY HALF SIZE (2**(M-1) CMPLX
C
C MASS STORE ARRAY MUST BE OPEN FOR ACCESS BY 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 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
C
      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
C
  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
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CMFFT
C COMPLEX FFT OF MULTI-DIMENSD MASS STORE ARRAY, CALLED BY USER OR RMFFT
C FOR COMMENTS, SEE SUBRTN RMFFT; ARGUMENTS HAVE SAME MEANING EXCEPT
C FOR IDIR=+1 OR -1, ARRAY ALWAYS COMPLEX, DIMENSION ORDER REVERSED
C WHILE IDIR=0, DIMENSION ORDER (AND MEXA LIST) ARE NOT REVERSED
C-----------------------------------------------------------------------
C
      SUBROUTINE CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA, IBEX, ICEX)
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
C MOST EFFICIENT USE OF CORE STORE - TRIES TO DO C-B PASSES PER LOAD
C
      IDUM = MFINDX(LSET,B,M,M,NPAS)
C
C DUMMY CALL TO MFINDX TO SPECIFY VIRTUAL (B.'S'.M)**(C-B) PERMUTATION
C
C FIRST, PIECE-MEAL ATTACK ON FFT COMPUTATION FOLLOWS
C
  10  IF (MPAS.LE.0) GO TO 40
      IF (MPAS.LT.NPAS) NPAS = MPAS
  20  CALL MFLOAD(LREAD, BUFA, IBEX, ICEX, IFLG)
C
C LOAD CORE WORKING SPACE WITH 2**(C-B) BLOCKS ACCORDING TO MFINDX
C
      IF (IFLG.LT.0) GO TO 30
      CALL MFCOMP(MEXA, NDIM, ISGN, BUFA, B, C, M, NPAS, IPAS)
C
C DO MODIFIED IN-CORE FFT, REQUIRING NPAS PASSES STARTING WITH IPAS
C
      CALL MFLOAD(LWRIT, BUFA, IBEX, ICEX, IFLG)
C
C UNLOAD CORE AREA, WRITING BLOCKS BACK IN-PLACE TO MASS STORE
C
      GO TO 20
  30  IPAS = IPAS + NPAS
      MPAS = MPAS - NPAS
      GO TO 10
C
C END OF FIRST PART
C
C SPECIFY BLOCKS TO BE READ IN NEXT PART IN TRUE ORDER (NO PERM)
C
  40  IDUM = MFINDX(LSET,B,M,M,0)
C
C FINAL, CONCLUDING ATTACK ON FFT COMPUTATION FOLLOWS
C
  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
C DO FINAL, C-PASS IN-CORE FFT OF EACH CORE-LOAD
C
      IF (SCAL.EQ.1.) GO TO 70
      DO 60 J=1,NC
        BUFA(J) = BUFA(J)*SCAL
  60  CONTINUE
  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)
C
  80  IF (IDIR.EQ.0) GO TO 90
      CALL MFSORT(BUFA, IBEX, ICEX, 1, MREAL, MREAL)
      M = MFSUM(MEXA,NDIM,-1)
C
C DO FULL BIT-REVERSAL OF M BITS (AND REVERSE MEXA LIST)
C
      RETURN
C
C BELOW, REVERSE BITS OF EACH DIMEN SEPARATELY (NO DIMEN REVERSAL)
C
  90  IH = 1
      DO 100 J=1,NDIM
        IG = IH
        IH = IH + MEXA(J)
        CALL MFSORT(BUFA, IBEX, ICEX, IG, IH, MREAL)
 100  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  MFCOMP
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
      SUBROUTINE MFCOMP(MEXA, NDIM, ISGN, BUFA, B, C, M, NPAS, IPAS)
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
      PI = 4.0D0 * DATAN(1.0D0)
      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
C KPAS IS GLOBAL COMPUTING PASS NUMBER (JPAS-1 IS LOCAL)
C KPEFF IS EFFECTIVE GLOBAL PASS FOR MULTIDIMEN. FFT W GENERATION
C D IS USED FOR RECURSIVE MODIFICATION OF W PHASE FACTOR
C SPAN SEPARATES VALUES IN FFT KERNEL, STEP TO NEXT PAIR, SAME W
C NRPT COUNTS REPETITION OF W FOR MULTIDIMEN. FFT
C IMOD AND NCLR ARE USED TO COMPUTE W WITHOUT EXCEEDING SMALL INTEGER
C
      DO 80 JPAS=1,NPAS
        KPAS = IPAS + JPAS - 1
        KDIFF = M - MFSUM(MEXA,NDIM,KPAS)
        KPEFF = KPAS + KDIFF
        D = CEXP(CMPLX(0.,PIMOD*2.0**KPEFF))
        ITEM = C - JPAS
        SPAN = 2**ITEM
        STEP = 2*SPAN
        IF (B.LT.ITEM) ITEM = B
        NB = 2**ITEM
        IF (ITEM.GT.KDIFF) ITEM = KDIFF
        NRPT = 2**ITEM
        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
C BELOW, START OF ONE PASS THROUGH CORE, NOTING BLOCK BOUNDARIES
C
  20    DO 70 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
C NEW W COMPUTED DIRECTLY AT BEGINNING OF NEW BLOCK AREA
C
C BELOW, COMPUTATIONS WITHIN EACH BLOCK OF 2**B CMPLX ELMTS
C
  30      DO 60 I2=1,NB,NRPT
C
C BELOW, REPETITION OF SAME W DUE TO MULTIDIMEN FFT
C
            DO 50 I3=1,NRPT
              I4 = I1 + I2 + I3 - 2
C
C BELOW, STEPPING THOUGH INDICES HAVING SAME W IN ONE DIMEN FFT
C
              DO 40 J=I4,NC,STEP
                K = J + SPAN
                TEMP = (BUFA(J)-BUFA(K))*W
                BUFA(J) = BUFA(J) + BUFA(K)
                BUFA(K) = TEMP
  40          CONTINUE
  50        CONTINUE
C
C FFT 2-POINT KERNEL ARITHMETIC (ALGORITHM BIT-REVERSAL FOLLOWS COMPUT)
C
            W = W*D
C
C RECURSIVE MODIFICATION OF W WITHIN BLOCK BOUNDARIES
C
  60      CONTINUE
  70    CONTINUE
  80  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  MFSORT
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-----------------------------------------------------------------------
C
      SUBROUTINE MFSORT(BUFA, IBEX, ICEX, IG, IH, M)
C
      REAL BUFA(1)
      DATA LSET /1/, LPERM /2/, LREAD /1/, LWRIT /2/
C
      IDUM = MFINDX(LSET,IBEX,M,M,0)
C
C DUMMY CALL TO INITIALISE MFINDX (INITIALLY UNPERMUTED ARRAY)
C
      IF (IH-IG.LE.1) RETURN
      IF (IG.GE.IBEX) GO TO 50
      IF (IH.LE.ICEX) GO TO 60
C
C CHECK FOR SPECIAL CASES, REQUIRING SIMPLER TREATMENT
C
C BELOW, MIXED PERMUTATION OF BOTH ELEMENTS AND BLOCKS
C MOST EFFICIENT USE OF CORE STORE - TRIES TO DO ICEX-IBEX PASSES
C PER LOAD
C
      IPAS = 0
      NPAS = ICEX - IBEX
      MPAS = IH - IBEX
      IF (IBEX-IG.LT.MPAS) MPAS = IBEX - IG
C
C BELOW, FIRST VIRTUAL 'S' PERMUTATIONS
C
  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
C BYPASS UNNECESSARY CORE LOAD IF TRIVIAL CASES
C
      IDUM = MFINDX(LPERM,IBEX,IH,M,NPAS)
C
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
C
  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
C CARRY OUT IN-CORE, SYMMETRIC R PERMS. (ONE ONLY ON FIRST PASS)
C
      CALL MFLOAD(LWRIT, BUFA, IBEX, ICEX, IFLG)
C
C UNLOAD CORE AREA, WRITING BLOCKS BACK IN-PLACE TO MASS STORE
C
      GO TO 20
C
  30  IPAS = IPAS + NPAS
      MPAS = MPAS - NPAS
      GO TO 10
C
C END OF FIRST PART
C
C BELOW, FINAL 'R' PERM. OF BLOCKS IN MASS STORE, IF (IH-2*IBEX+IG).GT.1
C
  40  CALL MFREV(BUFA, 0, IH-2*IBEX+IG, M-IBEX, IBEX)
      RETURN
C
C BELOW, PERMUTATION OF BLOCKS ONLY REQUIRED
C
  50  CALL MFREV(BUFA, IG-IBEX, IH-IBEX, M-IBEX, IBEX)
      RETURN
C
C BELOW, PERMUTATION OF ELEMENTS IN CORE ONLY REQUIRED
C
  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
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: MFREV
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
      SUBROUTINE MFREV(BUFA, IG, IH, M, IBEX)
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
C NO PERMUTATION REQUIRED
C
      NB = 2**IBEX
      NB1 = NB + 1
      NG = 2**IG
      NGDB = NG*2
      NH = 2**IH
      NHHF = NH/2
C
C NG IS MOVEMENT GROUP SIZE, NH IS PERMUTATION REPLICATION SIZE
C
      NM = 2**M
      NPARS = NM - NH + 1
      NREV = NH/4
C
      DO 10 J=1,IHG
        IRA(J) = 0
        NREV = NREV/2
        NRA(J) = NREV
  10  CONTINUE
C
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
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)
C
  20  NOF = 3
      IF (IFOR.GE.IREV) NOF = 1
C
C SELECTS ONCE-ONLY SWAP PAIRS (EITHER 1 OR 3 PAIRS)
C
      DO 60 JOF=1,NOF
        IFOF = IFOFA(JOF)
        IROF = IROFA(JOF)
        DO 50 I1=1,NG
C
C REPETITION OVER GROUP OR SUPER ELEMENT OF NG ACTUAL ELEMENTS
C
          IN2F = IFOR + IFOF + I1
          IN2R = IREV + IROF + I1
          DO 40 I2=1,NPARS,NH
C
C REPETITION OF SAME PERMUTATION OVER ARRAY PARTS
C
            IN3F = IN2F + I2
            IN3R = IN2R + I2
            IF (IBEX.GE.0) GO TO 30
C
C BELOW, IN-CORE ELEMENT SORTING
C
            TEMP = BUFA(IN3R)
            BUFA(IN3R) = BUFA(IN3F)
            BUFA(IN3F) = TEMP
            GO TO 40
C
C BELOW, SORTING WHOLE BLOCKS IN MASS STORE
C
  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
  50    CONTINUE
  60  CONTINUE
C
C END OF INNER, REPETITION LOOPS
C
      IFOR = IFOR + NGDB
C
C INCREMENT FORWARD QUARTER-LENGTH INTEGER (ALREADY SCALED BY NG*2)
C
      IF (IFOR.GE.NHHF) RETURN
      IF (IREV.GE.NREV) GO TO 70
C
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
C
  70  DO 80 J=1,IHG
        IF (IRA(J).LT.NRA(J)) GO TO 90
  80  CONTINUE
C
C BELOW, SIMPLE INCREMENT OF REVERSE INTEGER, LOWER IN HIERARCHY
C
  90  IRA(J) = IRA(J) + NRA(J)
      IREV = IRA(J)
 100  IF (J.EQ.1) GO TO 20
      J = J - 1
      IRA(J) = IREV
      GO TO 100
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  MFLOAD
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-----------------------------------------------------------------------
C
      SUBROUTINE MFLOAD(LOAD, BUFA, IBEX, ICEX, IFLG)
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
C HOLDS CURRENT MFINDX VALUE FOR ENTRY 2 AND SUBRTN MFCOMP
C
      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
C READS BLOCK WITH NEXT VIRTUAL MFINDX INDEX
C
  10  CONTINUE
      RETURN
  20  IFLG = -1
      RETURN
C
C RESETS MFINDX TO START OF IN-PLACE BLOCK
C
  30  IDUM = MFINDX(LREST,0,0,0,0)
      DO 40 J=1,NCB
        K = (J-1)*NB
        JB = MFINDX(LINDX,0,0,0,0)
        CALL MFWRIT(BUFA(K+1), NB, JB)
C
C WRITES BLOCK WITH NEXT VIRTUAL MFINDX INDEX (REPEAT MFREAD SEQUENCE)
C
  40  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:  MFINDX
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
      FUNCTION MFINDX(LSPEC, B, H, M, N)
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
      COMMON /VARIAB/ JAY, JOFF, NRPT, IPERM, JAYH, JOFH, NRPTH, NHB,
     *    NMB, IHB, ISTEP, NPERM
      INTEGER B, H
C
      IF (LSPEC.EQ.1) GO TO 50
      IF (LSPEC.EQ.2) GO TO 60
      IF (LSPEC.EQ.3) GO TO 70
      IF (LSPEC.EQ.4) GO TO 90
C
      IF (ISTEP.NE.0) GO TO 10
C
C BELOW, PRECEDES FIRST MFINDX OF A PASS
C
      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
C
  10  MFINDX = JAY + JOFF
      IF (MFINDX.GT.NMB) GO TO 30
      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
C
      JOFF = JOFF + NHB
  20  NRPT = NHB
      JAY = 0
      RETURN
C
C 40 BELOW, END OF ONE PASS, PARS RESET, IPERM ALTERED IF INVERSE
C
  30  IF (NPERM.LT.0) IPERM = MOD(IPERM-NPERM,IHB)
      MFINDX = -1
  40  JOFF = 1
      ISTEP = 0
      GO TO 20
C
C LSPEC=1 (LSET) SETS IPERM=0 (UNPERMED), ENTERS B,H,M,N PARAMS
C
  50  IPERM = 0
C
C LSPEC=2 (LPERM) CHANGES THE B,H,M,N PARAMETERS (DUMMIES ELSEWHERE)
C
  60  IHB = H - B
      NPERM = N
      NHB = 2**IHB
      NMB = 2**(M-B)
      MFINDX = IPERM
      GO TO 40
C
C LSPEC=3 (LHOLD) HOLDS CURRENT MFINDX INDEXING PARAMETERS
C
  70  JAYH = JAY
      JOFH = JOFF
      NRPTH = NRPT
  80  MFINDX = IPERM
      RETURN
C
C LSPEC=4 (LREST) RESTORES PARAMETERS TO INDEX MFINDX AT LAST LHOLD
C
  90  JAY = JAYH
      JOFF = JOFH
      NRPT = NRPTH
      GO TO 80
      END
C
C-----------------------------------------------------------------------
C FUNCTION:  MFSUM
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-----------------------------------------------------------------------
C
      FUNCTION MFSUM(MEXA, NDIM, MLIM)
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
C
      NDIMH = NDIM/2
      DO 20 J=1,NDIMH
        K = NDIM + 1 - J
        MTEM = MEXA(J)
        MEXA(J) = MEXA(K)
        MEXA(K) = MTEM
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  MFRCMP
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-----------------------------------------------------------------------
C
      SUBROUTINE MFRCMP(MEXA, NDIM, ISGN, IDIR, IPAK, BUFA, B, M)
C
      COMPLEX ATEM, BTEM, TEMP, W, D, BUFA(1)
      INTEGER B, JAYA(4), KAYA(4), JWKA(4), KWKA(4), MEXA(1)
C
C JAYA,KAYA,JWKA,KWKA ALLOW UP TO 4 DIMENSIONS - INCREASE IF REQUIRED
C
      DATA LOWER /1/, LUPPR /2/, LCLR /-1/
      PI = 4.0D0 * DATAN(1.0D0)
C
      DO 10 IDIM=1,NDIM
        JAYA(IDIM) = 0
        KAYA(IDIM) = 0
  10  CONTINUE
C
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
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
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
C CONJUGATE-SYMMETRIC BASE INDICES COMPUTED FROM MULTIDIMEN. SET
C
  40  MEX1 = MEXA(1)
C
C 2**MEX1 IS NUMBER OF VALUES ADJACENT IN FIRST DIMENSION
C
      IFLG = -1
      IF (MEX1.LE.B) GO TO 160
C
C BELOW, FIRST DIMEN. GREATER THAN BLOCK SIZE, MULTIPLE BLOCKS
C
      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
C JB AND KB ARE BLOCK INDEX PAIRS CONTAINING CONJUGATE ELEMENTS
C
      J1 = 0
      K1 = 0
C
  50  NCNT = NB + 1
  60  CALL MFRLOD(LOWER, IOF, BUFA, NB, JB, JBMAX, JBOF, IDIR, NCNT)
C
C LOWER BLOCK LOADED
C
      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
C UPPER BLOCK LOADED
C
      IFLG = IFLG + 1
      K2 = K1 + IOF
      KB = KB + KBINC
C
C FIRST TIME, UPPER BLOCK STEPS HIGH, FOLLOWING STEPS SMALL NEGATIVE
C
      KBINC = -1
C
C J AND K INDEX CONJUGATE-SYMMETRIC PAIRS IN CORE
C
  80  J = J2 + J3
      K = K2 + K3
      JJ = J + NBDB
      KK = K + NBDB
      IF (IDIR.GT.0) GO TO 200
C
C BELOW, UNSCRAMBLING FOR REAL-TO-COMPLEX FFT
C
      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 170
      BUFA(K) = CONJG(BTEM)
      IF (IEXPND.GT.0) BUFA(KK) = CONJG(ATEM)
C
  90  J3 = J3 + 1
      K3 = KINC - J3
C
C IN-CORE INDEX PAIRS STEPPED IN OPPOSING DIRECTIONS
C
      IF (IDIM.NE.NDIM) GO TO 100
      IWFG = 1
      W = W*D
C
C RECURSIVE MODIFICATION OF W IF UNIDIMEN. TRANSFORM
C
 100  NCNT = NCNT - 1
      IF (NCNT.LE.0) GO TO 110
C
C ENTER RECURSION ROUTINE IF OPERATION COMPLETE IN CURRENT DIMEN
C
      IF (J3.EQ.1) GO TO 70
      IF (NCNT.GT.NEWBLK) GO TO 80
C
C END OF INNER LOOP (NOTE SPECIAL CASE WHEN J3=1 ABOVE)
C
C BELOW, MAY REQUIRE TO READ NEW BLOCKS
C
      NBLCNT = NBLCNT - 1
      IF (NBLCNT.GT.0) GO TO 50
      IF (JAY.EQ.KAY) GO TO 60
C
C JAY.EQ.KAY NEEDS SYMMETRICAL MIDDLE, OTHERWISE CURRENT DIMEN COMPLT
C
C BELOW, RECURSION TO COMPUTE MULTIDIMEN. CONJUGATE-SYMMETRY
C
 110  JAYA(IDIM) = 0
      KAYA(IDIM) = 0
      IDIM = IDIM + 1
      IF (IDIM.GT.NDIM) GO TO 140
      NUMD = 2**MEXA(IDIM)
      IF (NUMD.LE.1) GO TO 140
      IF (IDIM.NE.NDIM) GO TO 120
      IWFG = 1
      W = W*D
C
C RECURSIVE MODIFICATION OF W IF MULTIDIMEN. FFT
C
 120  IF (JAYA(IDIM).EQ.0) GO TO 130
      IF ((JWKA(IDIM)*NUMD+JAYA(IDIM)).EQ.(KWKA(IDIM)*NUMD+KAYA(IDIM)))
     *    GO TO 110
      IF (KAYA(IDIM).EQ.1) GO TO 110
C
 130  JAYA(IDIM) = JAYA(IDIM) + 1
      KAYA(IDIM) = NUMD - JAYA(IDIM)
C
C RECURSIVE STEPPING OF MULTIDIMEN. CONJUGATE-SYMMETRIC INDEX PAIRS
C
      GO TO 20
C
C BELOW, OPERATION COMPLETE, TIDY UP AND RETURN FROM SUBROUTINE
C
 140  DO 150 IAREA=1,2
        CALL MFRLOD(IAREA, IOF, BUFA, NB, LCLR, JBMAX, JBOF, IDIR, IFLG)
C
C DUMMY CALL TO MFRLOD TO WRITE ANY UNWRITTEN BLOCKS TO MASS STORE
C
 150  CONTINUE
      RETURN
C
C RETURN FROM SUBROUTINE
C
C BELOW, FIRST DIMEN. LESS THAN BLOCK SIZE, INDEX PAIRS ALL IN-CORE
C
 160  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
C JB AND KB ARE BLOCK INDEX PAIRS CONTAINING CONJUGATE ELEMENTS
C
      J1 = (JAY-(JB-1)*NBPD1)*NUMD1
      K1 = (KAY-(KB-1)*NBPD1)*NUMD1
      GO TO 60
C
C BELOW, UNSCRAMBLING WITH W0 (IWFG=0) MUST BE TREATED DIFFERENTLY
C
 170  IF (IEXPND.LT.0) GO TO 180
C
C BELOW, ARRAY EXPANSION (EITHER IPAK=+1 OR IPAK=0 AND STILL REDUNDANT)
C
      BUFA(K) = CONJG(ATEM)
      BUFA(KK) = CONJG(BTEM)
      GO TO 90
C
C BELOW, NO ARRAY EXPANSION (EITHER IPAK=-1 OR IPAK=0 NOT REDUNDANT)
C
 180  IF (J.EQ.K) GO TO 190
      BUFA(K) = CONJG(BTEM)
      GO TO 90
C
C BELOW, SPECIAL CASE IF IPAK=-1 AND ELEMENTS ARE SAME
C
 190  BUFA(J) = CMPLX(REAL(ATEM),REAL(BTEM))
      GO TO 90
C
C BELOW, SCRAMBLING FOR COMPLEX-TO-REAL FFT
C
 200  IF (IWFG.EQ.0) GO TO 220
      BTEM = CONJG(BUFA(K))
 210  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 BELOW, SCRAMBLING WITH W0 (IWFG=0) MUST BE TREATED DIFFERENTLY
C
 220  IF (IEXPND.LT.0) GO TO 230
      BTEM = BUFA(JJ)
      GO TO 210
C
C BELOW, NO REDUNDANCY (EITHER IPAK=-1 OR IPAK=0 OR 1 NOT REDUND)
C
 230  IF (J.EQ.K) GO TO 240
      BTEM = CONJG(BUFA(K))
      GO TO 210
C
C BELOW, SPECIAL CASE IF IPAK=-1 AND ELEMENTS ARE SAME
C
 240  BTEM = CMPLX(AIMAG(BUFA(J)),0.)
      BUFA(J) = CMPLX(REAL(BUFA(J)),0.)
      GO TO 210
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  MFRLOD
C LOADS, UNLOADS CORE STORE ARRAY BUFA, FOR REAL FFT UNSCRAMBLING ROUTN
C-----------------------------------------------------------------------
C
      SUBROUTINE MFRLOD(IAREA, IOF, BUFA, NB, JB, JBMAX, JBOF, IDIR,
     *    NCNT)
C
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
C JBAREA() HOLDS INDEX OF BLOCK LOADED IN AREA 1 OR 2 (FIRST TIME -1 BEL
C
      IF (JBAREA(1).GE.0) GO TO 20
      IEXIST = -1
      DO 10 I=1,NPARF
        JBPFA(I) = -1
  10  CONTINUE
C
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
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
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
C IOF IS BASE OFFSET OF CORE AREA WHERE BLOCK IS TO BE FOUND
C
  50  IOF = (IAREA-1)*NB + 1
      IOFDB = IOF + NBDB
      IF (JBAREA(IAREA).LT.0) GO TO 90
      CALL MFWRIT(BUFA(IOF), NBDB, JBAREA(IAREA))
C
C WRITE OUT BLOCK BEFORE READING NEW BLOCK
C
      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
C
      DO 60 I=1,NPARF
        IF (JBPFA(I).LT.0) GO TO 70
  60  CONTINUE
C
C NO ROOM IN LISTS, DEFAULTS TO ALL READ
C
      I = NPARF
      IEXIST = I
  70  JBPFA(I) = JBAREA(IAREA)
      NCPFA(I) = NCAREA(IAREA)
C
  80  CALL MFWRIT(BUFA(IOFDB), NBDB, JBOF+JBAREA(IAREA))
C
C SIMILARLY, WRITE OUT BLOCK PAIR IF EXPANDING
C
C BELOW, READ BLOCK NOTING BLOCK INDX (READ EXPANDED, IF PART FILLED)
C
  90  JBAREA(IAREA) = JB
      IF (JB.LT.0) GO TO 130
      CALL MFREAD(BUFA(IOF), NBDB, JB)
C
C READ REQUIRED BLOCK AND NOTE ACCESS COUNT
C
      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
C
      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
C IF BLOCK TO BE READ WAS ONLY PART FILLED, THEN IT EXISTS TO READ
C
 120  CALL MFREAD(BUFA(IOFDB), NBDB, JBOF+JB)
C
C READ EXPANDED BLOCK IF REQUIRED
C
 130  RETURN
C
C RETURN FROM SUBROUTINE
C
C BELOW, TRIVIAL CASE - BLOCK ALREADY LOADED
C IOF IS BASE OFFSET OF CORE AREA WHERE BLOCK IS TO BE FOUND
C
 140  IOF = (I-1)*NB + 1
      IF (I.NE.IAREA) GO TO 150
      NCAREA(I) = NCAREA(I) + NCHLD
C
C INCREASE ACCESS COUNT IF CURRENT IAREA MATCHES ORIGINAL IAREA
C
 150  RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:  MFPAR
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
      FUNCTION MFPAR(IRMF, ICOMP)
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
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 NMAX /4/
      FIXMAX = FLOAT(I1MACH(9))
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
        MEXA(I) = IFIX(ALOG(DIMA(I))/ALG2+0.5)
  10  CONTINUE
C
  20  M = 0
      DO 30 I=1,NMAX
        IF (I.GT.NDIM) MEXA(I) = 0
        M = M + MEXA(I)
        DIMA(I) = 2.**MEXA(I)
  30  CONTINUE
      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
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
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
C IPAK=-1 REQUIRES NO FILE EXPANSION
C
      IDIM = NDIM
      IF (IRMF.LT.0) IDIM = 1
      FADD = SIZE*2./DIMA(IDIM)
C
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
C FBLK IS MAXIMUM NUMBER OF 'RECDS', SIZE RBLK, POSSIBLE
C
      TBLK = FBLK
      IF (IRMF.GE.0) GO TO 60
C
C GENERALLY TBLK=FBLK, BUT FOR PACKED REAL NOT SO, BELOW
C
      FSIZ = SIZE
      TBLK = FSIZ/RBLK
C
  60  TDM1 = 1.
      RDM1 = FSIZ
      IF (NDIM.EQ.1) GO TO 70
C
C JOB COMPLETED IF NDIM=1
C
      RDM1 = DIMA(1)
      IF (IRMF.GE.0) RDM1 = RDM1*2.
      TDM1 = FSIZ/RDM1
C
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
C CONVERT VALUES IN /MFVAL/ TO INTEGERS IN /MFINT/ (-1 IF TOO LARGE)
C
      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
C FORCES ICEX.NGT.M+1 AND IBEX.NGT.ICEX-1, OR MFPAR=1 (IRMF=0)
C
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  DMPERM
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
      SUBROUTINE DMPERM(MEXA, NDIM, NSHFT, IREX, BUFA, IBEX, ICEX)
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
C FINDS M TOTAL AND REVERSES MEXA LIST
C
      CALL MFSORT(BUFA, IBEX, ICEX, IREX, M, M)
C
C INITIAL OVERALL BIT-REVERSAL M BITS ABOVE IREX BITS
C
      IF (NSHFT.LT.0) GO TO 10
C
C BELOW, REVERSAL OF TWO PARTS, TO FORM REQUIRED SHIFT
C
      IH = MFSUM(MEXA,NS,-1) + IREX
      CALL MFSORT(BUFA, IBEX, ICEX, IREX, IH, M)
C
C REVERSE LOWER PART OF MEXA LIST AND LOWER PART OF ARRAY BITS
C
      CALL MFSORT(BUFA, IBEX, ICEX, IH, M, M)
C
C SEPARATELY REVERSE UPPER PART OF ARRAY BITS
C
      IH = MFSUM(MEXA(NSHFT+1),NDIM-NS,-1)
C
C REVERSE UPPER PART OF MEXA LIST
C
      RETURN
C
C RETURN FROM SUBROUTINE AFTER CYCLIC SHIFTS
C
C BELOW, SEPARATELY REVERSE OVER EACH DIMENSION (DIMEN REVERSAL)
C
  10  IH = IREX
      DO 20 J=1,NDIM
        IG = IH
        IH = IH + MEXA(J)
        CALL MFSORT(BUFA, IBEX, ICEX, IG, IH, M)
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C TEST PROGRAM:  UNIVERSAL
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
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)
      DIMENSION XR(32)
      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 STO
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 = I1MACH(2)
      IPRINT = 1
C
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
C M SETS THE OVERALL ARRAY SIZE FOR AUTO IBEX,ICEX,MEXA,NDIM STEPPING
C FIXED VALUES CAN BE USED (SEE COMMENTS BELOW DO LOOP)
C
      IRMF = -1
C
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN = -1
      IDIR = -1
      IPAK = 1
C
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 SQUEEZED RESULT
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      IF (IPRINT.LE.0) WRITE (LP,9999) IPRINT, IRMF
9999  FORMAT (20H1MASS STORE FFT TEST/9H IPRINT =, I3, 13H (1=COPIOUS,0,
     *    22H=MAX DIFFS,-1=OVERALL)/9H   IRMF =, I3, 15H (0=CMFFT TEST,,
     *    14H-1=RMFFT TEST)/)
      DIFMG = 0.
C
C BELOW, COMPUTE ALL POSSIBLE MEXA FOR 1,2 AND 3 DIMENSIONS
C
      DO 200 NDIM1=1,3
        NDIM = NDIM1
        M2M = M - NDIM + 1
        IF (NDIM.EQ.1) M2M = 1
        DO 190 M2=1,M2M
          M3M = M2M - M2 + 1
          IF (NDIM.NE.3) M3M = 1
          DO 180 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
C REPLACE DO LOOP BY FIXED MEXA LIST, IF DESIRED
C
            IBEX = 2
            ICEX = 4
C
C DUMMY IBEX, ICEX SO NO ERROR IN MFPAR BELOW
C
            IERR = MFPAR(IRMF,0)
            IF (IERR.NE.0) GO TO 210
C
C CALL HELPER ROUTINE TO COMPUTE SIZES USED BY NAIVE DFT SUBRTN
C
            M = MFSUM(MEXA,NDIM,99)
            AR = RANMF(1)
C
C RESET RANDOM NUMBER GENERATOR AND COMPUTE M IN CASE NOT GIVEN
C
            DO 10 J=1,NSZE
              RANDA(J) = RANMF(-1)
              ANAIV(J) = CMPLX(RANDA(J),0.)
  10        CONTINUE
C
C LOAD RANDOM NUMBERS FOR FFT AND FOR NAIVE SUBROUTINE
C
            CALL NAIVE(ANAIV, BNAIV, NDMA(1), NDMA(2), NDMA(3), ISGN)
C
C NAIVE DFT SUBROUTINE CALLED TO COMPUTE 'SLOW' FOURIER TRANSFORM
C
            IF (IRMF.EQ.0) GO TO 20
C
C BELOW, SET LIMTS FOR STEPPING IBEX, ICEX, WHEN CALLING RMFFT
C
            IBEXL = 2
            ICEXL = 4
            ICEXM = M
            GO TO 30
C
C SETS DIFFERENT LIMITS FOR IBEX, ICEX FOR CMFFT BELOW (IRMF.EQ.0)
C
  20        IBEXL = 1
            ICEXL = 2
            ICEXM = M + 1
C
  30        IF (ICEXL.GT.ICEXM) GO TO 210
            DO 170 ICEX1=ICEXL,ICEXM
              ICEX = ICEX1
              IBEXM = ICEX - 2
              IF (IRMF.EQ.0) IBEXM = ICEX - 1
              DO 160 IBEX1=IBEXL,IBEXM
                IBEX = IBEX1
C
C IBEX AND ICEX COMPUTED; REPLACE DO LOOPS BY FIXED IBEX,ICEX IF REQ
C
                IERR = MFPAR(IRMF,0)
                IF (IERR.NE.0) GO TO 210
                NCNT = NRD1
C
C HELPER ROUTN, NTD1 TOTAL NUMB OF NRD1 REALS IN FIRST DIMENSION
C (NCNT IS NUMBER OF REAL ELMTS IN FIRST DIMENSION)
C
                SCAL = 1.
                IF (IRMF.EQ.0) GO TO 60
C
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
                DO 50 JB=1,NTD1
                  K = (JB-1)*NCNT
                  DO 40 I=1,NCNT
                    J = K + I
                    BUFA(I) = RANDA(J)
  40              CONTINUE
                  CALL MFWRIT(BUFA, NRD1, JB)
  50            CONTINUE
C
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
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
                GO TO 90
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FIRST DI
C
  60            NCNT = NCNT/2
                DO 80 JB=1,NTD1
                  K = (JB-1)*NCNT
                  DO 70 I=1,NCNT
                    J = K + I
                    CBUFA(I) = CMPLX(RANDA(J),0.)
  70              CONTINUE
                  CALL MFWRIT(BUFA, NRD1, JB)
  80            CONTINUE
C
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
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
  90            IF (IPRINT.LE.0) GO TO 100
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
                WRITE (LP,9999) IPRINT, IRMF
                WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *              (MEXA(J),J=1,NDIM)
9998            FORMAT (8H  NDIM =, I3, 8H  ISGN =, I3, 8H  IDIR =, I3,
     *              8H  IBEX =, I3, 8H  ICEX =, I3, 8H  IPAK =,
     *              I3/10H  MEXA() =, 3I5)
                WRITE (LP,9997)
9997            FORMAT (8H   INDEX, 9X, 3HFFT, 19X, 5HNAIVE, 18X,
     *              4HDIFF)
C
 100            IERR = MFPAR(-IRMF,0)
                IF (IERR.NE.0) GO TO 210
                NCNT = NRD1/2
C
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
C
                DIFM = 0.
                DO 120 JB=1,NTD1
                  K = (JB-1)*NCNT
                  CALL MFREAD(BUFA, NRD1, JB)
C
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 110 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 110
                    IF (IDIR.NE.0) WRITE (LP,9996) INDEX, CBUFA(I),
     *                  BNAIV(J), CDIF
                    IF (IDIR.EQ.0) WRITE (LP,9996) INDEX, CBUFA(I),
     *                  ANAIV(J), CDIF
9996                FORMAT (1X, I5, 2(1X, 2E11.3),1X,2E10.2)
 110              CONTINUE
 120            CONTINUE
C
C BELOW, PRINT INTERMEDIATE DIFFERENCES
C
                IF (IPRINT.GE.0) WRITE (LP,9995) DIFM, NDIM, ISGN,
     *              IDIR, IBEX, ICEX, IPAK, (MEXA(J),J=1,NDIM)
9995            FORMAT (10H MAX DIFF , E11.4/8H  NDIM =, I3, 8H  ISGN =,
     *              I3, 8H  IDIR =, I3, 8H  IBEX =, I3, 8H  ICEX =, I3,
     *              8H  IPAK =, I3/10H  MEXA() =, 3I5)
                IF (DIFM.GT.DIFMG) DIFMG = DIFM
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
C
                ISGN = -ISGN
                IDIR = -IDIR
                SCAL = 1./SIZE
                IF (IRMF.NE.0) CALL RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL,
     *              BUFA, IBEX, ICEX, IPAK)
C
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
                IF (IRMF.EQ.0) CALL CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL,
     *              BUFA, IBEX, ICEX)
C
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
                IF (IPRINT.LE.0) GO TO 130
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
                WRITE (LP,9999) IPRINT, IRMF
                WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *              (MEXA(J),J=1,NDIM)
                WRITE (LP,9994)
9994            FORMAT (8H   INDEX, 2X, 8HFFT/2**M, 4X, 5HINPUT, 7X,
     *              4HDIFF)
C
 130            IERR = MFPAR(IRMF,0)
                IF (IERR.NE.0) GO TO 210
                NCNT = NRD1
                IF (IRMF.EQ.0) NCNT = NCNT/2
C
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
C
                DIFM = 0.
                DO 150 JB=1,NTD1
                  K = (JB-1)*NCNT
                  CALL MFREAD(BUFA, NRD1, JB)
C
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 140 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,9996) INDEX, RM,
     *                  RANDA(J), DIF
 140              CONTINUE
 150            CONTINUE
C
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RANDOM SET)
C
                IF (IPRINT.GE.0) WRITE (LP,9995) DIFM, NDIM, ISGN,
     *              IDIR, IBEX, ICEX, IPAK, (MEXA(J),J=1,NDIM)
                IF (DIFM.GT.DIFMG) DIFMG = DIFM
                ISGN = -ISGN
                IDIR = -IDIR
C
C RESTORE ISGN AND IDIR FOR FORWARD TRANSFORM
C
 160          CONTINUE
 170        CONTINUE
 180      CONTINUE
 190    CONTINUE
 200  CONTINUE
C
C BELOW, PRINT OVERALL MAXIMUM DIFFERENCE
C
      WRITE (LP,9993) DIFMG, M, ISGN, IDIR, IPAK
9993  FORMAT (18H OVERALL MAX DIFF , E11.4/25H FOR ALL MEXA(1 TO 3 DIM),
     *    23H,IBEX,ICEX,MEXA FOR M =, I3/19H   ISGN,IDIR(+/-) =, 2I4,
     *    8H  IPAK =, I4)
      STOP
C
 210  WRITE (LP,9992) IERR
9992  FORMAT (25H IBEX FORCED TOO SMALL OR, 24H NOT ALL /MFINT/ CORRECT,
     *    6H, IERR, I4)
      STOP
      END
C
C-----------------------------------------------------------------------
C FUNCTION:  RANMF
C RANDOM NUMBER GENERATOR FOR MASS STORE FFT TEST
C-----------------------------------------------------------------------
C
      FUNCTION RANMF(J)
      IF (J.GE.0) GO TO 20
C
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
        K = MOD(5*K,MODULO)
  10  CONTINUE
      Z = FLOAT(K)/FLMOD
      RANMF = Z
      RETURN
C
  20  K = J
      RANMF = J
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  NAIVE
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-----------------------------------------------------------------------
C
      SUBROUTINE NAIVE(ANAIV, BNAIV, NJ, NK, NL, ISGN)
      COMPLEX TEMP, ANAIV(NJ,NK,NL), BNAIV(NL,NK,NJ)
C
      PI = 4.0D0 * DATAN(1.0D0)
      PI2 = PI*2.0
      IF (ISGN.LT.0) PI2 = -PI2
C
      DO 60 JB=1,NJ
        AJB = FLOAT(JB-1)/FLOAT(NJ)
        DO 50 KB=1,NK
          AKB = FLOAT(KB-1)/FLOAT(NK)
          DO 40 LB=1,NL
            ALB = FLOAT(LB-1)/FLOAT(NL)
C
            TEMP = (0.,0.)
            DO 30 JA=1,NJ
              AJA = FLOAT(JA-1)*AJB
              DO 20 KA=1,NK
                AKA = FLOAT(KA-1)*AKB
                DO 10 LA=1,NL
                  ALA = FLOAT(LA-1)*ALB
                  TEMP = TEMP + ANAIV(JA,KA,LA)*CEXP(CMPLX(0.,PI2*
     *                (AJA+AKA+ALA)))
  10            CONTINUE
  20          CONTINUE
  30        CONTINUE
C
            BNAIV(LB,KB,JB) = TEMP
  40      CONTINUE
  50    CONTINUE
  60  CONTINUE
C
      DO 90 JA=1,NJ
        DO 80 KA=1,NK
          DO 70 LA=1,NL
            ANAIV(JA,KA,LA) = BNAIV(LA,KA,JA)
  70      CONTINUE
  80    CONTINUE
  90  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  MFREAD
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-----------------------------------------------------------------------
C
      SUBROUTINE MFREAD(BUFA, NB, JB)
      COMMON /MASS/ RMAS(1024)
      REAL BUFA(NB)
C
      IOF = (JB-1)*NB
      DO 10 I=1,NB
        K = IOF + I
        BUFA(I) = RMAS(K)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: MFWRIT
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-----------------------------------------------------------------------
C
      SUBROUTINE MFWRIT(BUFA, NB, JB)
      COMMON /MASS/ RMAS(1024)
      REAL BUFA(NB)
C
      IOF = (JB-1)*NB
      DO 10 I=1,NB
        K = IOF + I
        RMAS(K) = BUFA(I)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C TEST PROGRAM:  MASTOM
C CONTROL DATA 6000 AND CYBER MASS STORE I/O FFT
C---------------------------------------------------------------------
C
      PROGRAM MASTOM(TAPE1,INPUT,OUTPUT,TAPE60=INPUT,TAPE5=OUTPUT)
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
C LUN IS LOGICAL UNIT FOR CYBER MASS STORE I/O
C
      LP = 5
      IPRINT = 0
C
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF = -1
C
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN = -1
      IDIR = -1
      IPAK = 1
C
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
C
      IF (IPRINT.LE.0) WRITE (LP,9999) IPRINT, IRMF
9999  FORMAT (31H1MASS STORE FFT TEST - IPRINT =, I3, 14H (1=COPIOUS,0=,
     *    11HMAX DIFFS),, 9H   IRMF =, I3, 25H (0=CMFFT TEST,1=RMFFT TE,
     *    3HST)/)
      DIFMG = 0.
C
      MEXA(1) = 8
      MEXA(2) = 6
      NDIM = 2
      IBEX = 9
      ICEX = 13
C
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 130
      CALL OPENMS(LUN, MINDX, NFBK+1, 0)
C
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
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 30
C
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 20 JB=1,NTBK
        DO 10 I=1,NCNT
          BUFA(I) = VALU
          VALU = VALU + VINC
  10    CONTINUE
        CALL MFWRIT(BUFA, NRBK, JB)
  20  CONTINUE
C
C LOAD RAMP FUNCTION 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
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
      GO TO 60
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FFT BLOC
C
  30  NCNT = NCNT/2
      DO 50 JB=1,NTBK
        DO 40 I=1,NCNT
          CBUFA(I) = CMPLX(VALU,0.)
          VALU = VALU + VINC
  40    CONTINUE
        CALL MFWRIT(BUFA, NRBK, JB)
  50  CONTINUE
C
C LOAD RAMP FUNCTION 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
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
  60  IF (IPRINT.LE.0) GO TO 90
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
9998  FORMAT (8H  NDIM =, I3, 8H  ISGN =, I3, 8H  IDIR =, I3, 7H  IBEX ,
     *    1H=, I3, 8H  ICEX =, I3, 8H  IPAK =, I3, 10H  MEXA() =, 3I5)
      WRITE (LP,9997)
9997  FORMAT (8H   INDEX, 12X, 3HFFT)
C
      IERR = MFPAR(-IRMF,0)
      IF (IERR.NE.0) GO TO 130
      NCNT = NRBK/2
C
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
C
      DO 80 JB=1,NTBK
        K = (JB-1)*NCNT
        CALL MFREAD(BUFA, NRBK, JB)
C
C READ 'RECDS' OF NRBK LENGTH, TOTALING 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 70 I=1,NCNT
          INDEX = J - 1
          IF (IPRINT.LE.0) GO TO 70
          WRITE (LP,9996) INDEX, CBUFA(I)
9996      FORMAT (1X, I5, 2X, 2E13.4)
  70    CONTINUE
  80  CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
C
  90  ISGN = -ISGN
      IDIR = -IDIR
      SCAL = 1./SIZE
      IF (IRMF.NE.0) CALL RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX, IPAK)
C
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF (IRMF.EQ.0) CALL CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX)
C
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF (IPRINT.LE.0) GO TO 100
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
      WRITE (LP,9995)
9995  FORMAT (8H   INDEX, 4X, 8HFFT/2**M, 6X, 5HINPUT, 10X, 4HDIFF)
C
 100  IERR = MFPAR(IRMF,0)
      IF (IERR.NE.0) GO TO 130
      NCNT = NRBK
      IF (IRMF.EQ.0) NCNT = NCNT/2
C
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
C
      DIFM = 0.
      VALU = 0.
      DO 120 JB=1,NTBK
        CALL MFREAD(BUFA, NRBK, JB)
C
C READ 'RECDS' OF NRBK LENGTH, TOTALING 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 110 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,9994) I, RM, VALU, DIF
9994      FORMAT (1X, I5, 3(2X, E13.4))
          VALU = VALU + VINC
 110    CONTINUE
 120  CONTINUE
C
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF (IPRINT.GE.0) WRITE (LP,9993) DIFM, NDIM, ISGN, IDIR, IBEX,
     *    ICEX, IPAK, (MEXA(J),J=1,NDIM)
9993  FORMAT (10H MAX DIFF , E11.4, 1X, 8H  NDIM =, I3, 8H  ISGN =, I3,
     *    8H  IDIR =, I3, 8H  IBEX =, I3, 8H  ICEX =, I3, 8H  IPAK =,
     *    I3, 10H  MEXA() =, 3I5)
      IF (DIFM.GT.DIFMG) DIFMG = DIFM
C
      STOP
C
 130  WRITE (LP,9992) IERR
9992  FORMAT (25H IBEX FORCED TOO SMALL OR, 24H NOT ALL /MFINT/ CORRECT,
     *    6H, IERR, I4)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINES: MFREAD, MFWRIT
C CONTROL DATA 6000 AND CYBER MASS STORE I/O ROUTINES FOR FFT
C-----------------------------------------------------------------------
C
      SUBROUTINE MFREAD(BUFA, NB, JB)
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
C WRITE BLOCK, INDEX JB, FROM BUFA TO MASS STORE, NB REAL VALUES
C
      CALL WRITMS(LUN, BUFA, NB, JB, -1, 0)
      RETURN
      END
C
C-----------------------------------------------------------------------
C TEST PROGRAM: LCMTOM
C CONTROL DATA 6000 AND CYBER EXTENDED CORE/LCM FFT
C-----------------------------------------------------------------------
C
      PROGRAM LCMTOM(INPUT,OUTPUT,TAPE60=INPUT,TAPE5=OUTPUT)
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
      LEVEL3, 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
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF = -1
C
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN = -1
      IDIR = -1
      IPAK = 1
C
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
C
      IF (IPRINT.LE.0) WRITE (LP,9999) IPRINT, IRMF
9999  FORMAT (31H1MASS STORE FFT TEST - IPRINT =, I3, 14H (1=COPIOUS,0=,
     *    11HMAX DIFFS),, 9H   IRMF =, I3, 25H (0=CMFFT TEST,1=RMFFT TE,
     *    3HST)/)
      DIFMG = 0.
C
      MEXA(1) = 8
      MEXA(2) = 6
      NDIM = 2
      IBEX = 7
      ICEX = 13
C
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 130
      NCNT = NRD1
C
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 30
C
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 20 JB=1,NTD1
        DO 10 I=1,NCNT
          BUFA(I) = VALU
          VALU = VALU + VINC
  10    CONTINUE
        CALL MFWRIT(BUFA, NRD1, JB)
  20  CONTINUE
C
C LOAD RAMP FUNCTION 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 STOR
C
      CALL RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA, IBEX, ICEX, IPAK)
C
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
      GO TO 60
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FFT BLOC
C
  30  NCNT = NCNT/2
      DO 50 JB=1,NTD1
        DO 40 I=1,NCNT
          CBUFA(I) = CMPLX(VALU,0.)
          VALU = VALU + VINC
  40    CONTINUE
        CALL MFWRIT(BUFA, NRD1, JB)
  50  CONTINUE
C
C LOAD RAMP FUNCTION 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 STOR
C
      CALL CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA, IBEX, ICEX)
C
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
  60  IF (IPRINT.LE.0) GO TO 90
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
9998  FORMAT (8H  NDIM =, I3, 8H  ISGN =, I3, 8H  IDIR =, I3, 7H  IBEX ,
     *    1H=, I3, 8H  ICEX =, I3, 8H  IPAK =, I3, 10H  MEXA() =, 3I5)
      WRITE (LP,9997)
9997  FORMAT (8H   INDEX, 12X, 3HFFT)
C
      IERR = MFPAR(-IRMF,0)
      IF (IERR.NE.0) GO TO 130
      NCNT = NRD1/2
C
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
C
      DO 80 JB=1,NTD1
        K = (JB-1)*NCNT
        CALL MFREAD(BUFA, NRD1, JB)
C
C READ 'RECDS' OF NRD1 LENGTH, TOTALING 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 STOR
C
        DO 70 I=1,NCNT
          INDEX = J - 1
          IF (IPRINT.LE.0) GO TO 70
          WRITE (LP,9996) INDEX, CBUFA(I)
9996      FORMAT (1X, I5, 2X, 2E13.4)
  70    CONTINUE
  80  CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
C
  90  ISGN = -ISGN
      IDIR = -IDIR
      SCAL = 1./SIZE
      IF (IRMF.NE.0) CALL RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX, IPAK)
C
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF (IRMF.EQ.0) CALL CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX)
C
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF (IPRINT.LE.0) GO TO 100
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
      WRITE (LP,9995)
9995  FORMAT (8H   INDEX, 4X, 8HFFT/2**M, 6X, 5HINPUT, 10X, 4HDIFF)
C
 100  IERR = MFPAR(IRMF,0)
      IF (IERR.NE.0) GO TO 130
      NCNT = NRD1
      IF (IRMF.EQ.0) NCNT = NCNT/2
C
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
C
      DIFM = 0.
      VALU = 0.
      DO 120 JB=1,NTD1
        CALL MFREAD(BUFA, NRD1, JB)
C
C READ 'RECDS' OF NRD1 LENGTH, TOTALING 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 STOR
C
        DO 110 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,9994) I, RM, VALU, DIF
9994      FORMAT (1X, I5, 3(2X, E13.4))
          VALU = VALU + VINC
 110    CONTINUE
 120  CONTINUE
C
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF (IPRINT.GE.0) WRITE (LP,9993) DIFM, NDIM, ISGN, IDIR, IBEX,
     *    ICEX, IPAK, (MEXA(J),J=1,NDIM)
9993  FORMAT (10H MAX DIFF , E11.4, 1X, 8H  NDIM =, I3, 8H  ISGN =, I3,
     *    8H  IDIR =, I3, 8H  IBEX =, I3, 8H  ICEX =, I3, 8H  IPAK =,
     *    I3, 10H  MEXA() =, 3I5)
      IF (DIFM.GT.DIFMG) DIFMG = DIFM
C
      STOP
C
 130  WRITE (LP,9992) IERR
9992  FORMAT (25H IBEX FORCED TOO SMALL OR, 24H NOT ALL /MFINT/ CORRECT,
     *    6H, IERR, I4)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINES: MFREAD, MFWRIT
C CONTROL DATA 6000 AND CYBER EXTENDED CORE STORE I/O ROUTINES FOR FFT
C-----------------------------------------------------------------------
C
      SUBROUTINE MFREAD(BUFA, NB, JB)
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
      LEVEL3, LBUFA
      COMMON /FFTCOM/ LBUFA(32768)
      REAL BUFA(NB)
C
      CALL MOVLEV(LBUFA((JB-1)*NB+1), BUFA, NB)
      RETURN
C
      ENTRY MFWRIT
C
C WRITE BLOCK, INDEX JB, FROM BUFA TO EXTENDED CORE, NB REAL VALUES
C
      CALL MOVLEV(BUFA, LBUFA((JB-1)*NB+1), NB)
      RETURN
      END
C
C-----------------------------------------------------------------------
C TEST PROGRAM: PDP 11 MASS STORE I/O FFT SAMPLE
C-----------------------------------------------------------------------
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
C LUN IS LOGICAL UNIT FOR PDP 11 MASS STORE I/O
C
      LP = 5
      IPRINT = 0
C
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF = -1
C
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN = -1
      IDIR = -1
      IPAK = 1
C
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
C
      IF (IPRINT.LE.0) WRITE (LP,9999) IPRINT, IRMF
9999  FORMAT (31H1MASS STORE FFT TEST - IPRINT =, I3, 14H (1=COPIOUS,0=,
     *    11HMAX DIFFS),, 9H   IRMF =, I3, 25H (0=CMFFT TEST,1=RMFFT TE,
     *    3HST)/)
      DIFMG = 0.
C
      MEXA(1) = 8
      MEXA(2) = 6
      NDIM = 2
      IBEX = 7
      ICEX = 11
C
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 130
      CALL ASSIGN(LUN,'FFTEST.DAT',0)
      NWD = NRBK*2
      DEFINE FILE LUN(NFBK,NWD,U,INDX)
C
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
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 30
C
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 20 JB=1,NTBK
        DO 10 I=1,NCNT
          BUFA(I) = VALU
          VALU = VALU + VINC
  10    CONTINUE
        CALL MFWRIT(BUFA, NRBK, JB)
  20  CONTINUE
C
C LOAD RAMP FUNCTION 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
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
      GO TO 60
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN FFT BLOC
C
  30  NCNT = NCNT/2
      DO 50 JB=1,NTBK
        DO 40 I=1,NCNT
          CBUFA(I) = CMPLX(VALU,0.)
          VALU = VALU + VINC
  40    CONTINUE
        CALL MFWRIT(BUFA, NRBK, JB)
  50  CONTINUE
C
C LOAD RAMP FUNCTION 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
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
  60  IF (IPRINT.LE.0) GO TO 90
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
9998  FORMAT (8H  NDIM =, I3, 8H  ISGN =, I3, 8H  IDIR =, I3, 7H  IBEX ,
     *    1H=, I3, 8H  ICEX =, I3, 8H  IPAK =, I3, 10H  MEXA() =, 3I5)
      WRITE (LP,9997)
9997  FORMAT (8H   INDEX, 12X, 3HFFT)
C
      IERR = MFPAR(-IRMF,0)
      IF (IERR.NE.0) GO TO 130
      IF (IERR.NE.0) GO TO 130
      NCNT = NRBK/2
C
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
C
      DO 80 JB=1,NTBK
        K = (JB-1)*NCNT
        CALL MFREAD(BUFA, NRBK, JB)
C
C READ 'RECDS' OF NRBK LENGTH, TOTALING 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 70 I=1,NCNT
          INDEX = J - 1
          IF (IPRINT.LE.0) GO TO 70
          WRITE (LP,9996) INDEX, CBUFA(I)
9996      FORMAT (1X, I5, 2X, 2E13.4)
  70    CONTINUE
  80  CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
C
  90  ISGN = -ISGN
      IDIR = -IDIR
      SCAL = 1./SIZE
      IF (IRMF.NE.0) CALL RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX, IPAK)
C
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF (IRMF.EQ.0) CALL CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX)
C
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF (IPRINT.LE.0) GO TO 100
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
      WRITE (LP,9995)
9995  FORMAT (8H   INDEX, 4X, 8HFFT/2**M, 6X, 5HINPUT, 10X, 4HDIFF)
C
 100  IERR = MFPAR(IRMF,0)
      IF (IERR.NE.0) GO TO 130
      NCNT = NRBK
      IF (IRMF.EQ.0) NCNT = NCNT/2
C
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
C
      DIFM = 0.
      VALU = 0.
      DO 120 JB=1,NTBK
        CALL MFREAD(BUFA, NRBK, JB)
C
C READ 'RECDS' OF NRBK LENGTH, TOTALING 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 110 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,9994) I, RM, VALU, DIF
9994      FORMAT (1X, I5, 3(2X, E13.4))
          VALU = VALU + VINC
 110    CONTINUE
 120  CONTINUE
C
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF (IPRINT.GE.0) WRITE (LP,9993) DIFM, NDIM, ISGN, IDIR, IBEX,
     *    ICEX, IPAK, (MEXA(J),J=1,NDIM)
9993  FORMAT (10H MAX DIFF , E11.4, 1X, 8H  NDIM =, I3, 8H  ISGN =, I3,
     *    8H  IDIR =, I3, 8H  IBEX =, I3, 8H  ICEX =, I3, 8H  IPAK =,
     *    I3, 10H  MEXA() =, 3I5)
      IF (DIFM.GT.DIFMG) DIFMG = DIFM
C
      STOP
C
 130  WRITE (LP,9992) IERR
9992  FORMAT (25H IBEX FORCED TOO SMALL OR, 24H NOT ALL /MFINT/ CORRECT,
     *    6H, IERR, I4)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  MFREAD
C PDP 11 FORTRAN DIRECT ACCESS READ ROUTINE FOR FFT
C-----------------------------------------------------------------------
C
      SUBROUTINE MFREAD(BUFA, NB, JB)
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
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: MFWRIT
C PDP 11 FORTRAN DIRECT ACCESS WRITE ROUTINE FOR FFT
C-----------------------------------------------------------------------
C
      SUBROUTINE MFWRIT(BUFA, NB, JB)
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
      END
C
C-----------------------------------------------------------------------
C TEST PROGRAM: PDP 11 FAST MACRO I/O FFT SAMPLE
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 EFFICIENT 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),
C CHANNEL (RT11) (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
C LUN IS LOGICAL UNIT FOR PDP 11 MASS STORE I/O
C
      LP = 5
      IPRINT = 0
C
C LP IS PRINTER LOGICAL UNIT, IPRINT=+1 COPIOUS, 0 MAX DIFFERENCES ONLY
C
      IRMF = -1
C
C IRMF=-1 REAL ROUTINE RMFFT TEST, 0 COMPLEX ROUTINE CMFFT TEST
C
      ISGN = -1
      IDIR = -1
      IPAK = 1
C
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
C
      IF (IPRINT.LE.0) WRITE (LP,9999) IPRINT, IRMF
9999  FORMAT (31H1MASS STORE FFT TEST - IPRINT =, I3, 14H (1=COPIOUS,0=,
     *    11HMAX DIFFS),, 9H   IRMF =, I3, 25H (0=CMFFT TEST,1=RMFFT TE,
     *    3HST)/)
      DIFMG = 0.
C
      MEXA(1) = 8
      MEXA(2) = 6
      NDIM = 2
      IBEX = 7
      ICEX = 11
C
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 130
C
C (NOTE HELPER ROUTN MFPAR RETURNS NFBK AS MAXIMUM FILE SIZE)
C
      STOP 'ASSIGN AND OPEN FILE HERE'
C
C DELETE ABOVE LINE AND INVOKE SYSTEM MACRO 'OPEN$' (RSX) OR '.OPEN (RT1
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
C PRESET IOERR IN /FFTCOM/; ZERO IF NO ERRORS IN I/O
C
      NCNT = NRD1
C
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 30
C
C SWITCH FOR COMPLEX ROUTINE CMFFT AT 50, REAL RMFFT BELOW
C
      DO 20 JB=1,NTD1
        DO 10 I=1,NCNT
          BUFA(I) = VALU
          VALU = VALU + VINC
  10    CONTINUE
        CALL MFWRIT(BUFA, NRD1, JB)
  20  CONTINUE
C
C LOAD RAMP FUNCTION 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 ST
C
      CALL RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA, IBEX, ICEX, IPAK)
C
C REAL MASS STORE ROUTINE RMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
      GO TO 60
C
C BELOW, USE COMPLEX ROUTINE CMFFT,  NCNT NUM OF CMPLX ELMTS IN 'FIRST'
C
  30  NCNT = NCNT/2
      DO 50 JB=1,NTD1
        DO 40 I=1,NCNT
          CBUFA(I) = CMPLX(VALU,0.)
          VALU = VALU + VINC
  40    CONTINUE
        CALL MFWRIT(BUFA, NRD1, JB)
  50  CONTINUE
C
C LOAD RAMP FUNCTION 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 ST
C
      CALL CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA, IBEX, ICEX)
C
C COMPLEX MASS STORE ROUTINE CMFFT TO TRANSFORM ARRAY TO COMPLEX RESULT
C
  60  IF (IPRINT.LE.0) GO TO 90
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
9998  FORMAT (8H  NDIM =, I3, 8H  ISGN =, I3, 8H  IDIR =, I3, 7H  IBEX ,
     *    1H=, I3, 8H  ICEX =, I3, 8H  IPAK =, I3, 10H  MEXA() =, 3I5)
      WRITE (LP,9997)
9997  FORMAT (8H   INDEX, 12X, 3HFFT)
C
      IERR = MFPAR(-IRMF,0)
      IF (IERR.NE.0) GO TO 130
      NCNT = NRD1/2
C
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
C
      DO 80 JB=1,NTD1
        K = (JB-1)*NCNT
        CALL MFREAD(BUFA, NRD1, JB)
C
C READ 'RECDS' OF NRD1 LENGTH, TOTALING 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 ST
C
        DO 70 I=1,NCNT
          INDEX = J - 1
          IF (IPRINT.LE.0) GO TO 70
          WRITE (LP,9996) INDEX, CBUFA(I)
9996      FORMAT (1X, I5, 2X, 2E13.4)
  70    CONTINUE
  80  CONTINUE
C
C BELOW, INVERT ISGN AND IDIR FOR INVERSE TRANSFORM (COMPLEX-TO-REAL)
C
  90  ISGN = -ISGN
      IDIR = -IDIR
      SCAL = 1./SIZE
      IF (IRMF.NE.0) CALL RMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX, IPAK)
C
C EITHER ROUTINE RMFFT INVERSE TRANSFORMS ARRAY TO PACKED REAL
C
      IF (IRMF.EQ.0) CALL CMFFT(MEXA, NDIM, ISGN, IDIR, SCAL, BUFA,
     *    IBEX, ICEX)
C
C OR ROUTINE CMFFT INVERSE TRANSFORMS ARRAY
C
      IF (IPRINT.LE.0) GO TO 100
C
C BELOW, PRINT HEADINGS FOR TEST OUTPUT
C
      WRITE (LP,9999) IPRINT, IRMF
      WRITE (LP,9998) NDIM, ISGN, IDIR, IBEX, ICEX, IPAK,
     *    (MEXA(J),J=1,NDIM)
      WRITE (LP,9995)
9995  FORMAT (8H   INDEX, 4X, 8HFFT/2**M, 6X, 5HINPUT, 10X, 4HDIFF)
C
 100  IERR = MFPAR(IRMF,0)
      IF (IERR.NE.0) GO TO 130
      NCNT = NRD1
      IF (IRMF.EQ.0) NCNT = NCNT/2
C
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
C
      DIFM = 0.
      VALU = 0.
      DO 120 JB=1,NTD1
        CALL MFREAD(BUFA, NRD1, JB)
C
C READ 'RECDS' OF NRD1 LENGTH, TOTALING 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 ST
C
        DO 110 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,9994) I, RM, VALU, DIF
9994      FORMAT (1X, I5, 3(2X, E13.4))
          VALU = VALU + VINC
 110    CONTINUE
 120  CONTINUE
C
C PRINT INVERSE RESULTS (SHOULD BE SAME AS INITIAL RAMP)
C
      IF (IPRINT.GE.0) WRITE (LP,9993) DIFM, NDIM, ISGN, IDIR, IBEX,
     *    ICEX, IPAK, (MEXA(J),J=1,NDIM)
9993  FORMAT (10H MAX DIFF , E11.4, 1X, 8H  NDIM =, I3, 8H  ISGN =, I3,
     *    8H  IDIR =, I3, 8H  IBEX =, I3, 8H  ICEX =, I3, 8H  IPAK =,
     *    I3, 10H  MEXA() =, 3I5)
      IF (DIFM.GT.DIFMG) DIFMG = DIFM
C
      STOP
C
 130  WRITE (LP,9992) IERR
9992  FORMAT (25H IBEX FORCED TOO SMALL OR, 24H NOT ALL /MFINT/ CORRECT,
     *    6H, IERR, I4)
      STOP
      END
;
;------------------------------------------------------------------
;  SUBROUTINE: MFREAD, MFWRIT
;  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
Return to Main Software Page