-
Notifications
You must be signed in to change notification settings - Fork 0
/
mfftis.f
90 lines (88 loc) · 2.63 KB
/
mfftis.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
SUBROUTINE MFFTIS(C,IMSX,IVSX,IESX,NMX,NVX,NEX,TABLES,IERR)
* PURPOSE:
* THE SAME AS MFFTIM. IT IS A VARIANT OF MFFTIM, OPTIMIZED FOR
* MAXIMUM PERFORMANCE ON MATRICES WHOSE FIRST DIMENSION IS
* LESS THAN 64.
* IT REQUIRES THAT TABLES HAS BEEN PREPARED BY A CALL TO MFFTP
* WITH ID .NE. 0 .
* IMPLEMENTATION:
* THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
* "BUTTERFLY" TRANSFORMATION MFFT?9; PARAMETERS OF THE "BUTTERFLY"
* ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
* ARGUMENTS:
* INPUT:
* C : ARRAY TO BE TRANSFORMED.
* IMSX,IVSX,IESX,NMX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE
* OF C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED
* ON OUTPUT
* TABLES : ARRAY PREPARED BY MFFTP. IT IS NOT CHANGED ON OUTPUT.
* IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
* IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
* OUTPUT:
* C : TRANSFORM OF THE ORIGINAL ARRAY; "BIT REVERSED" ORDER
* IERR : ERROR CODE : =0 : SUCCESSFUL
* : =3 : 'TABLES' NOT CORRECTLY INITIALIZED
COMPLEX C(*)
INTEGER TABLES(-14:*)
INTEGER IDERR,FACERR,TBERR
PARAMETER (IDERR=1,FACERR=2,TBERR=3)
COMMON /MFFTPA/ IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
$ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
SAVE IBASE
* LOADING THE COMMON BLOCK : CONSTANTS
IMS=IMSX
IVS=IVSX
IES=1
NM=NMX
NV=NVX
NE=NEX
IVLIM=(NV-1)*IVS
LSTEP=IMS
* LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
LX=1
MX=NM
IBASE=(2+IMS)*NM
IFAC=TABLES(-1)
IF(IFAC.GT.3)THEN
IERR=TBERR
RETURN
ENDIF
*.. RADIX 2 LOOP
200 CONTINUE
DO 210 IM=1,TABLES(-14)
MX=MX/2
NUSTEP=LX*LSTEP
ILIM=NUSTEP-1
MSTEP=NUSTEP*2
MLIM=NM*LSTEP-MSTEP
CALL MFFTA9(C,TABLES(2*IBASE))
LX=LX+LX
IBASE=IBASE+NUSTEP
210 CONTINUE
IF(IFAC.EQ.1)RETURN
*.. RADIX 3 LOOP
300 CONTINUE
DO 310 IM=1,TABLES(-13)
MX=MX/3
NUSTEP=LX*LSTEP
ILIM=NUSTEP-1
MSTEP=NUSTEP*3
MLIM=NM*LSTEP-MSTEP
CALL MFFTB9(C,TABLES(2*IBASE))
LX=LX*3
IBASE=IBASE+NUSTEP*2
310 CONTINUE
IF(IFAC.EQ.2)RETURN
*.. RADIX 5 LOOP
500 CONTINUE
DO 510 IM=1,TABLES(-12)
MX=MX/5
NUSTEP=LX*LSTEP
ILIM=NUSTEP-1
MSTEP=NUSTEP*5
MLIM=NM*LSTEP-MSTEP
CALL MFFTC9(C,TABLES(2*IBASE))
LX=LX*5
IBASE=IBASE+NUSTEP*4
510 CONTINUE
END