fft551.inc
References to this file elsewhere.
1 !DECK FFT551
2 ! SUBROUTINE 'FFT551' - MULTIPLE FAST COSINE TRANSFORM
3 !
4 ! AUTHOR: CLIVE TEMPERTON, MAY 1988
5 ! [ALL-FORTRAN VERSION: C.T., OCTOBER 1995]
6 !
7 ! COSINE TRANSFORM OF LENGTH N IS CONVERTED TO
8 ! REAL PERIODIC TRANSFORM OF LENGTH N BY PRE- AND POST-
9 ! PROCESSING. REAL PERIODIC TRANSFORM IS PERFORMED BY
10 ! PRUNING REDUNDANT OPERATIONS FROM COMPLEX TRANSFORM.
11 !
12 ! SEE FOR EXAMPLE PAUL SWARZTRAUBER, "SYMMETRIC FFT'S",
13 ! MATH. COMP. 47 (1986), 323-346.
14 !
15 ! A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA
16 ! WORK IS AN AREA OF SIZE (N+1)*MIN(LOT,64)
17 ! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
18 ! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N
19 ! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
20 ! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
21 ! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
22 ! N+1 IS THE LENGTH OF THE DATA VECTORS
23 ! (WHICH INCLUDE NONZERO VALUES AT BOTH ENDPOINTS)
24 ! LOT IS THE NUMBER OF DATA VECTORS
25 ! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
26 ! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
27 !
28 ! ORDERING OF COEFFICIENTS: Z(0) , Z(1) , Z(2) , ... , Z(N)
29 !
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=0,...,N)(E(J)*Z(J)*COS(I*J*PI/N))
41 ! WHERE E(J)=0.5 FOR J=0,N --- ELSE E(J)=1
42 !
43 ! ISIGN=-1: Z(J)=(2/N)*SUM(I=0,...,N)(E(I)*X(I)*COS(I*J*PI/N))
44 !
45 ! N.B. FFT551 has an unusual definition of the FFTs,
46 ! such that the the coeff of wave0 is NOT the mean.
47 !
48 !---------------------------------------------------------------------
49 Subroutine FFT551 & ! in
50 ( ISIGN, & ! in
51 INC, & ! in
52 JUMP, & ! in
53 LOT, & ! in
54 N, & ! in
55 IFAX, & ! in
56 TRIGS, & ! in
57 A, & ! inout
58 IDIM ) ! in
59
60 ! Code Description: ORIGINAL CODE F77 IS HARDLY TOUCHED !!!
61
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
71 Integer , intent (in) :: IFAX(10) ! previously prepared list of
72 ! factors of N
73
74 Real , intent (in) :: TRIGS(3*N) ! previously prepared list of
75 ! trigonometric function values
76 Real , intent (inout) :: A( INC*(N+1) + JUMP*(LOT-1) ) ! data array ! vectors (which include zeros
77 ! at the endpoints)
78 Integer , intent (in) :: IDIM ! dimension workspace
79
80 Real :: WORK(IDIM) ! size (n+1)*min(lot,VectorLength)
81 Integer :: NFAX,NX,NH
82 Integer :: NBLOX,NVEX,NB
83 Integer :: K, IC, J, LA, IGO, JA,JB,IA,IB
84 Integer :: IFAC,IERR,ISTART
85
86 Real :: CO,S, t1,t2,si,scale, vectorlength
87
88 CHARACTER (LEN=*), PARAMETER :: RoutineName = "Var_FFT551"
89
90 VectorLength = LOT
91 NFAX=IFAX(1)
92 NX=N+1
93 NH=N/2
94 NBLOX=1+(LOT-1)/VectorLength
95 NVEX=LOT-(NBLOX-1)*VectorLength
96 ISTART=1
97 !
98 DO 200 NB=1,NBLOX
99 !
100 ! PREPROCESSING
101 ! -------------
102 IA=ISTART
103 IB=IA+NH*INC
104 IC=IA+N*INC
105 JA=1
106 JB=NH+1
107 IF (MOD(NFAX,2).EQ.1) THEN
108 !DIR$ IVDEP
109 DO 105 J=1,NVEX
110 T1=0.5*(A(IA)+A(IC))
111 T2=0.5*(A(IA)-A(IC))
112 A(IA)=T1
113 A(IC)=T2
114 IA=IA+JUMP
115 IC=IC+JUMP
116 105 CONTINUE
117 ELSE
118 !DIR$ IVDEP
119 DO 110 J=1,NVEX
120 WORK(JA)=0.5*(A(IA)+A(IC))
121 WORK(JB)=A(IB)
122 A(IC)=0.5*(A(IA)-A(IC))
123 IA=IA+JUMP
124 IB=IB+JUMP
125 IC=IC+JUMP
126 JA=JA+NX
127 JB=JB+NX
128 110 CONTINUE
129 ENDIF
130 !
131 DO 130 K=1,NH-1
132 JA=K+1
133 JB=N+1-K
134 IA=ISTART+K*INC
135 IB=ISTART+(JB-1)*INC
136 IC=ISTART+N*INC
137 SI=TRIGS(2*N+K)
138 CO=TRIGS(2*N+NH-K)
139 IF (MOD(NFAX,2).EQ.1) THEN
140 !DIR$ IVDEP
141 DO 115 J=1,NVEX
142 T1 = 0.5*(A(IA)+A(IB)) - SI*(A(IA)-A(IB))
143 T2 = 0.5*(A(IA)+A(IB)) + SI*(A(IA)-A(IB))
144 A(IC) = A(IC) + CO*(A(IA)-A(IB))
145 A(IA) = T1
146 A(IB) = T2
147 IA=IA+JUMP
148 IB=IB+JUMP
149 IC=IC+JUMP
150 115 CONTINUE
151 ELSE
152 !DIR$ IVDEP
153 DO 120 J=1,NVEX
154 WORK(JA) = 0.5*(A(IA)+A(IB)) - SI*(A(IA)-A(IB))
155 WORK(JB) = 0.5*(A(IA)+A(IB)) + SI*(A(IA)-A(IB))
156 A(IC) = A(IC) + CO*(A(IA)-A(IB))
157 IA=IA+JUMP
158 IB=IB+JUMP
159 IC=IC+JUMP
160 JA=JA+NX
161 JB=JB+NX
162 120 CONTINUE
163 ENDIF
164 130 CONTINUE
165 !
166 ! PERIODIC FOURIER ANALYSIS
167 ! -------------------------
168 IA=ISTART
169 LA=N
170 IGO=1-2*MOD(NFAX,2)
171 !
172 DO 140 K=1,NFAX
173 IFAC=IFAX(NFAX+2-K)
174 LA=LA/IFAC
175 IERR=-1
176 IF (IGO.EQ.+1) THEN
177 CALL qpassm(WORK,WORK(IFAC*LA+1),A(IA),A(IA+LA*INC), &
178 TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
179 ELSE IF (IGO.EQ.-1) THEN
180 CALL qpassm(A(IA),A(IA+IFAC*LA*INC),WORK,WORK(LA+1), &
181 TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
182 ENDIF
183 IF (IERR.NE.0) GO TO 500
184 IGO=-IGO
185 140 CONTINUE
186 !
187 ! POSTPROCESSING
188 ! --------------
189 SCALE=2.0
190 IF (ISIGN.EQ.+1) SCALE = FLOAT(N)
191 S=1.0
192 IF (ISIGN.EQ.-1) S = 2.0/FLOAT(N)
193 JA=ISTART
194 JB=JA+N*INC
195 IA=1
196 IB=N
197 !DIR$ IVDEP
198 DO 150 J=1,NVEX
199 A(JA)=SCALE*WORK(IA)
200 A(JA+INC)=S*A(JB)
201 A(JB)=SCALE*WORK(IB)
202 IA=IA+NX
203 IB=IB+NX
204 JA=JA+JUMP
205 JB=JB+JUMP
206 150 CONTINUE
207 !
208 DO 170 K=2,N-2,2
209 JA=ISTART+K*INC
210 IA=K
211 !DIR$ IVDEP
212 DO 160 J=1,NVEX
213 A(JA)=SCALE*WORK(IA)
214 A(JA+INC)=-SCALE*WORK(IA+1)+A(JA-INC)
215 IA=IA+NX
216 JA=JA+JUMP
217 160 CONTINUE
218 170 CONTINUE
219 !
220 ISTART=ISTART+NVEX*JUMP
221 NVEX=VectorLength
222 200 CONTINUE
223 GO TO 570
224 !
225 ! ERROR MESSAGES
226 ! --------------
227 500 CONTINUE
228 GO TO (510,530,550) IERR
229 510 CONTINUE
230 WRITE(UNIT=0,FMT='(A,I4,A)') &
231 'VECTOR LENGTH =',NVEX,', GREATER THAN VectorLength'
232 GO TO 570
233 530 CONTINUE
234 WRITE(UNIT=0,FMT='(A,I3,A)') &
235 'FACTOR =',IFAC,', NOT CATERED FOR'
236 GO TO 570
237 550 CONTINUE
238 WRITE(UNIT=0,FMT='(A,I3,A)') &
239 'FACTOR =',IFAC,', ONLY CATERED FOR IF LA*IFAC=N'
240 570 CONTINUE
241
242 RETURN
243 END SUBROUTINE FFT551
244
245