fft661.inc
References to this file elsewhere.
1 Subroutine FFT661 & ! in
2 ( ISIGN, & ! in
3 INC, & ! in
4 JUMP, & ! in
5 LOT, & ! in
6 N, & ! in
7 IFAX, & ! in
8 TRIGS, & ! in
9 A, & ! inout
10 DIM ) ! in
11 !
12 !
13 ! Description:
14 ! MULTIPLE FAST SINE TRANSFORM
15 ! (Originally called FFT661, then Var_SinTrans)
16 ! author: clive temperton, may 1988
17 ! (slightly modified for all-fortran version)
18 !
19 ! Sine transform of length n is converted to
20 ! Real periodic transform of length n by pre- and post-
21 ! processing. Real periodic transform is performed by
22 ! pruning redundant operations from complex transform.
23 !
24 ! see for example paul swarztrauber, "symmetric fft's",
25 ! math. comp. 47 (1986), 323-346.
26 !
27 ! Method:
28 !
29 ! ordering of coefficients: z(0) , z(1) , z(2) , ... , z(n)
30 ! ordering of data: x(0) , x(1) , x(2) , ... , x(n)
31 !
32 ! vectorization is achieved on cray by doing the transforms
33 ! in parallel
34 !
35 ! N must be composed of factors 2,3 & 5 and must be even
36 !
37 ! definition of transforms:
38 ! -------------------------
39 !
40 ! isign=+1: x(i)=sum(j=1,...,n-1)(z(j)*sin(i*j*pi/n))
41 !
42 ! isign=-1: z(j)=(2/n)*sum(i=1,...,n-1)(x(i)*sin(i*j*pi/n))
43 !
44 ! Current Code Owner: Andrew Lorenc
45 !
46 ! History:
47 ! Version Date Comment
48 ! ------- ---- -------
49 ! 0.1 14/12/93 Original code. Phil Andrews
50 ! 0.2 16/09/94 Small Modifications for the
51 ! incorporation in the VAR project. HB
52 ! 1.1 21/04/95 placed under control. JB
53 ! 1.2 01/06/95 Tracing added. JB
54 !
55 ! Code Description:
56 ! NB BECAUSE OF THE TRICKY NESTED LOOPS
57 ! ORIGINAL CODE F77 IS HARDLY TOUCHED !!!
58
59 Implicit none
60
61 ! Subroutine arguments
62 Integer , intent (in) :: ISIGN ! Switch forward (-1) or inverse (+1)
63 Integer , intent (in) :: INC ! increment within each data
64 ! vector (e.g. INC=1 for
65 ! consecutively stored data)
66 Integer , intent (in) :: Jump ! increment between start of
67 ! data vectors
68 Integer , intent (in) :: LOT ! Number of data vectors
69 Integer , intent (in) :: N ! N+1 is the length of the data
70 ! vectors (which include zeros
71 ! at the endpoints)
72 Integer , intent (in) :: DIM ! dimension workspace
73 Integer , intent (in) :: IFAX(10) ! previously prepared list of
74 ! factors of N
75
76 Real , intent (in) :: TRIGS(3*N) ! previously prepared list of
77 ! trigonometric function values
78 Real , intent (inout) :: A( INC*(N+1) + JUMP*(LOT-1) ) ! data array
79
80 ! No descriptions given
81 Integer :: NFAX,NX,NH
82 Integer :: NBLOX,NVEX,NB
83 Integer :: K,JA,JB,IA,IB,IGO,LA,J
84 Integer :: IFAC,IERR,ISTART
85
86 Real :: SI,T1,T2,SCALE, vectorlength
87 Real :: WORK(DIM) ! size (n+1)*min(lot,VectorLength)
88
89 VectorLength = LOT
90 NFAX=IFAX(1)
91 NX=N+1
92 NH=N/2
93 NBLOX=1+(LOT-1)/VectorLength
94 NVEX=LOT-(NBLOX-1)*VectorLength
95 ISTART=1
96 !
97 DO 200 NB=1,NBLOX
98 !
99 ! PREPROCESSING
100 ! -------------
101 DO 120 K=1,NH-1
102 JA=K+1
103 JB=N+1-K
104 IA=ISTART+K*INC
105 IB=ISTART+(JB-1)*INC
106 SI=TRIGS(2*N+K)
107 IF (MOD(NFAX,2).EQ.0) THEN
108 !DIR$ IVDEP
109 DO 110 J=1,NVEX
110 WORK(JA) = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB))
111 WORK(JB) = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB))
112 IA=IA+JUMP
113 IB=IB+JUMP
114 JA=JA+NX
115 JB=JB+NX
116 110 CONTINUE
117 ELSE
118 !DIR$ IVDEP
119 DO 115 J=1,NVEX
120 T1 = SI*(A(IA)+A(IB)) + 0.5*(A(IA)-A(IB))
121 T2 = SI*(A(IA)+A(IB)) - 0.5*(A(IA)-A(IB))
122 A(IA) = T1
123 A(IB) = T2
124 IA=IA+JUMP
125 IB=IB+JUMP
126 115 CONTINUE
127 ENDIF
128 120 CONTINUE
129
130 JA=1
131 JB=NH+1
132 IA=ISTART
133 IB=ISTART+NH*INC
134 IF (MOD(NFAX,2).EQ.0) THEN
135 !DIR$ IVDEP
136 DO 130 J=1,NVEX
137 WORK(JA)=0.0
138 WORK(JB)=2.0*A(IB)
139 IB=IB+JUMP
140 JA=JA+NX
141 JB=JB+NX
142 130 CONTINUE
143 IGO = +1
144 ELSE
145 !DIR$ IVDEP
146 DO 135 J=1,NVEX
147 A(IA)=0.0
148 A(IB)=2.0*A(IB)
149 IA=IA+JUMP
150 IB=IB+JUMP
151 135 CONTINUE
152 IGO = -1
153 ENDIF
154 !
155 ! PERIODIC FOURIER ANALYSIS
156 ! -------------------------
157 IA=ISTART
158 LA=N
159 !
160 DO 140 K=1,NFAX
161 IFAC=IFAX(NFAX+2-K)
162 LA=LA/IFAC
163 IERR=-1
164 IF (IGO.EQ.+1) THEN
165 CALL qpassm(WORK,WORK(IFAC*LA+1),A(IA),A(LA*INC+IA), &
166 TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
167 ELSE IF (IGO.EQ.-1) THEN
168 CALL qpassm(A(IA),A(IFAC*LA*INC+IA),WORK,WORK(LA+1), &
169 TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
170 ENDIF
171 IF (IERR.NE.0) GO TO 500
172 IGO=-IGO
173 140 CONTINUE
174 !
175 ! POSTPROCESSING
176 ! --------------
177 SCALE=2.0
178 IF (ISIGN.EQ.+1) SCALE = FLOAT(N)
179 JA=ISTART
180 JB=JA+N*INC
181 IA=1
182 !DIR$ IVDEP
183 DO 150 J=1,NVEX
184 A(JA)=0.0
185 A(JA+INC)=0.5*SCALE*WORK(IA)
186 A(JB)=0.0
187 IA=IA+NX
188 JA=JA+JUMP
189 JB=JB+JUMP
190 150 CONTINUE
191 !
192 DO 170 K=2,N-2,2
193 JA=ISTART+K*INC
194 IA=K
195 !DIR$ IVDEP
196 DO 160 J=1,NVEX
197 A(JA)=-SCALE*WORK(IA+1)
198 A(JA+INC)=SCALE*WORK(IA)+A(JA-INC)
199 IA=IA+NX
200 JA=JA+JUMP
201 160 CONTINUE
202 170 CONTINUE
203 !
204 ISTART=ISTART+NVEX*JUMP
205 NVEX=VectorLength
206 200 CONTINUE
207
208 Go To 570
209 !
210 ! ERROR MESSAGES
211 ! --------------
212 500 CONTINUE
213 GO TO (510,530,550) IERR
214 510 CONTINUE
215 WRITE(UNIT=0,FMT='(A,I5,A)') 'NVEX=', NVEX ,' GREATER THAN VectorLength'
216 GO TO 570
217 530 CONTINUE
218 WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, 'NOT CATERED FOR'
219 GO TO 570
220 550 CONTINUE
221 WRITE(UNIT=0,FMT='(A,I5,A)') 'IFAC=', IFAC, ' ONLY CATERED FOR IF LA*IFAC=N'
222 570 CONTINUE
223
224
225 RETURN
226
227
228 End subroutine FFT661
229
230
231