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