← Back to branch summary

~ubuntu-branches/debian/sid/abinit/sid

« back to all changes in this revision

Viewing changes to src/lib01fftnew/fftstp.F90

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2007-09-14 13:05:00 UTC
  • Revision ID: james.westby@ubuntu.com-20070914130500-1kzh2mrgo6ir4b6i
Tags: upstream-5.3.4.dfsg
Import upstream version 5.3.4.dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib01fftnew/fftstp.F90
 
1
         subroutine fftstp(mm,n1dfft,m,nn,n,zin,zout,trig,after,now,before,isign)
 
2
         use defs_basis
 
3
!       RESTRICTIONS on USAGE
 
4
!  Copyright (C) 2002-2007 Stefan Goedecker, CEA Grenoble
 
5
!  This file is distributed under the terms of the
 
6
!  GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
 
7
 
 
8
! NEW FFTSTP
 
9
 
 
10
        implicit real(dp) (a-h,o-z)
 
11
!Arguments ------------------------------------
 
12
        integer after,before,atn,atb
 
13
        integer :: mm,n1dfft,m,nn,n,now,isign
 
14
        real(dp) :: zin,zout,trig
 
15
        dimension trig(2,2048),zin(2,mm,m),zout(2,nn,n)
 
16
!Local variables-------------------------------
 
17
 
 
18
! *************************************************************************
 
19
        atn=after*now
 
20
        atb=after*before
 
21
 
 
22
!         sqrt(.5d0)
 
23
        rt2i=half_sqrt2
 
24
        if (now.eq.2) then
 
25
        ia=1
 
26
        nin1=ia-after
 
27
        nout1=ia-atn
 
28
        do 2001,ib=1,before
 
29
        nin1=nin1+after
 
30
        nin2=nin1+atb
 
31
        nout1=nout1+atn
 
32
        nout2=nout1+after
 
33
        do 2001,j=1,n1dfft
 
34
        r1=zin(1,j,nin1)
 
35
        s1=zin(2,j,nin1)
 
36
        r2=zin(1,j,nin2)
 
37
        s2=zin(2,j,nin2)
 
38
        zout(1,j,nout1)= r2 + r1
 
39
        zout(2,j,nout1)= s2 + s1
 
40
        zout(1,j,nout2)= r1 - r2
 
41
        zout(2,j,nout2)= s1 - s2
 
42
2001        continue
 
43
        do 2000,ia=2,after
 
44
        ias=ia-1
 
45
        if (2*ias.eq.after) then
 
46
                if (isign.eq.1) then
 
47
                        nin1=ia-after
 
48
                        nout1=ia-atn
 
49
                        do 2010,ib=1,before
 
50
                        nin1=nin1+after
 
51
                        nin2=nin1+atb
 
52
                        nout1=nout1+atn
 
53
                        nout2=nout1+after
 
54
                        do 2010,j=1,n1dfft
 
55
                        r1=zin(1,j,nin1)
 
56
                        s1=zin(2,j,nin1)
 
57
                        r2=zin(2,j,nin2)
 
58
                        s2=zin(1,j,nin2)
 
59
                        zout(1,j,nout1)= r1 - r2
 
60
                        zout(2,j,nout1)= s2 + s1
 
61
                        zout(1,j,nout2)= r2 + r1
 
62
                        zout(2,j,nout2)= s1 - s2
 
63
2010                        continue
 
64
                else
 
65
                        nin1=ia-after
 
66
                        nout1=ia-atn
 
67
                        do 2020,ib=1,before
 
68
                        nin1=nin1+after
 
69
                        nin2=nin1+atb
 
70
                        nout1=nout1+atn
 
71
                        nout2=nout1+after
 
72
                        do 2020,j=1,n1dfft
 
73
                        r1=zin(1,j,nin1)
 
74
                        s1=zin(2,j,nin1)
 
75
                        r2=zin(2,j,nin2)
 
76
                        s2=zin(1,j,nin2)
 
77
                        zout(1,j,nout1)= r2 + r1
 
78
                        zout(2,j,nout1)= s1 - s2
 
79
                        zout(1,j,nout2)= r1 - r2
 
80
                        zout(2,j,nout2)= s2 + s1
 
81
2020                        continue
 
82
                end if
 
83
        else if (4*ias.eq.after) then
 
84
                if (isign.eq.1) then
 
85
                        nin1=ia-after
 
86
                        nout1=ia-atn
 
87
                        do 2030,ib=1,before
 
88
                        nin1=nin1+after
 
89
                        nin2=nin1+atb
 
90
                        nout1=nout1+atn
 
91
                        nout2=nout1+after
 
92
                        do 2030,j=1,n1dfft
 
93
                        r1=zin(1,j,nin1)
 
94
                        s1=zin(2,j,nin1)
 
95
                        r=zin(1,j,nin2)
 
96
                        s=zin(2,j,nin2)
 
97
                        r2=(r - s)*rt2i
 
98
                        s2=(r + s)*rt2i
 
99
                        zout(1,j,nout1)= r2 + r1
 
100
                        zout(2,j,nout1)= s2 + s1
 
101
                        zout(1,j,nout2)= r1 - r2
 
102
                        zout(2,j,nout2)= s1 - s2
 
103
2030                        continue
 
104
                else
 
105
                        nin1=ia-after
 
106
                        nout1=ia-atn
 
107
                        do 2040,ib=1,before
 
108
                        nin1=nin1+after
 
109
                        nin2=nin1+atb
 
110
                        nout1=nout1+atn
 
111
                        nout2=nout1+after
 
112
                        do 2040,j=1,n1dfft
 
113
                        r1=zin(1,j,nin1)
 
114
                        s1=zin(2,j,nin1)
 
115
                        r=zin(1,j,nin2)
 
116
                        s=zin(2,j,nin2)
 
117
                        r2=(r + s)*rt2i
 
118
                        s2=(s - r)*rt2i
 
119
                        zout(1,j,nout1)= r2 + r1
 
120
                        zout(2,j,nout1)= s2 + s1
 
121
                        zout(1,j,nout2)= r1 - r2
 
122
                        zout(2,j,nout2)= s1 - s2
 
123
2040                        continue
 
124
                end if
 
125
        else if (4*ias.eq.3*after) then
 
126
                if (isign.eq.1) then
 
127
                        nin1=ia-after
 
128
                        nout1=ia-atn
 
129
                        do 2050,ib=1,before
 
130
                        nin1=nin1+after
 
131
                        nin2=nin1+atb
 
132
                        nout1=nout1+atn
 
133
                        nout2=nout1+after
 
134
                        do 2050,j=1,n1dfft
 
135
                        r1=zin(1,j,nin1)
 
136
                        s1=zin(2,j,nin1)
 
137
                        r=zin(1,j,nin2)
 
138
                        s=zin(2,j,nin2)
 
139
                        r2=(r + s)*rt2i
 
140
                        s2=(r - s)*rt2i
 
141
                        zout(1,j,nout1)= r1 - r2
 
142
                        zout(2,j,nout1)= s2 + s1
 
143
                        zout(1,j,nout2)= r2 + r1
 
144
                        zout(2,j,nout2)= s1 - s2
 
145
2050                        continue
 
146
                else
 
147
                        nin1=ia-after
 
148
                        nout1=ia-atn
 
149
                        do 2060,ib=1,before
 
150
                        nin1=nin1+after
 
151
                        nin2=nin1+atb
 
152
                        nout1=nout1+atn
 
153
                        nout2=nout1+after
 
154
                        do 2060,j=1,n1dfft
 
155
                        r1=zin(1,j,nin1)
 
156
                        s1=zin(2,j,nin1)
 
157
                        r=zin(1,j,nin2)
 
158
                        s=zin(2,j,nin2)
 
159
                        r2=(s - r)*rt2i
 
160
                        s2=(r + s)*rt2i
 
161
                        zout(1,j,nout1)= r2 + r1
 
162
                        zout(2,j,nout1)= s1 - s2
 
163
                        zout(1,j,nout2)= r1 - r2
 
164
                        zout(2,j,nout2)= s2 + s1
 
165
2060                        continue
 
166
                end if
 
167
        else
 
168
                itrig=ias*before+1
 
169
                cr2=trig(1,itrig)
 
170
                ci2=trig(2,itrig)
 
171
                nin1=ia-after
 
172
                nout1=ia-atn
 
173
                do 2090,ib=1,before
 
174
                nin1=nin1+after
 
175
                nin2=nin1+atb
 
176
                nout1=nout1+atn
 
177
                nout2=nout1+after
 
178
                do 2090,j=1,n1dfft
 
179
                r1=zin(1,j,nin1)
 
180
                s1=zin(2,j,nin1)
 
181
                r=zin(1,j,nin2)
 
182
                s=zin(2,j,nin2)
 
183
                r2=r*cr2 - s*ci2
 
184
                s2=r*ci2 + s*cr2
 
185
                zout(1,j,nout1)= r2 + r1
 
186
                zout(2,j,nout1)= s2 + s1
 
187
                zout(1,j,nout2)= r1 - r2
 
188
                zout(2,j,nout2)= s1 - s2
 
189
2090                continue
 
190
        end if
 
191
2000        continue
 
192
        else if (now.eq.4) then
 
193
        if (isign.eq.1) then
 
194
                ia=1
 
195
                nin1=ia-after
 
196
                nout1=ia-atn
 
197
                do 4001,ib=1,before
 
198
                nin1=nin1+after
 
199
                nin2=nin1+atb
 
200
                nin3=nin2+atb
 
201
                nin4=nin3+atb
 
202
                nout1=nout1+atn
 
203
                nout2=nout1+after
 
204
                nout3=nout2+after
 
205
                nout4=nout3+after
 
206
                do 4001,j=1,n1dfft
 
207
                r1=zin(1,j,nin1)
 
208
                s1=zin(2,j,nin1)
 
209
                r2=zin(1,j,nin2)
 
210
                s2=zin(2,j,nin2)
 
211
                r3=zin(1,j,nin3)
 
212
                s3=zin(2,j,nin3)
 
213
                r4=zin(1,j,nin4)
 
214
                s4=zin(2,j,nin4)
 
215
                r=r1 + r3
 
216
                s=r2 + r4
 
217
                zout(1,j,nout1) = r + s
 
218
                zout(1,j,nout3) = r - s
 
219
                r=r1 - r3
 
220
                s=s2 - s4
 
221
                zout(1,j,nout2) = r - s
 
222
                zout(1,j,nout4) = r + s
 
223
                r=s1 + s3
 
224
                s=s2 + s4
 
225
                zout(2,j,nout1) = r + s
 
226
                zout(2,j,nout3) = r - s
 
227
                r=s1 - s3
 
228
                s=r2 - r4
 
229
                zout(2,j,nout2) = r + s
 
230
                zout(2,j,nout4) = r - s
 
231
4001                continue
 
232
                do 4000,ia=2,after
 
233
                ias=ia-1
 
234
                if (2*ias.eq.after) then
 
235
                        nin1=ia-after
 
236
                        nout1=ia-atn
 
237
                        do 4010,ib=1,before
 
238
                        nin1=nin1+after
 
239
                        nin2=nin1+atb
 
240
                        nin3=nin2+atb
 
241
                        nin4=nin3+atb
 
242
                        nout1=nout1+atn
 
243
                        nout2=nout1+after
 
244
                        nout3=nout2+after
 
245
                        nout4=nout3+after
 
246
                        do 4010,j=1,n1dfft
 
247
                        r1=zin(1,j,nin1)
 
248
                        s1=zin(2,j,nin1)
 
249
                        r=zin(1,j,nin2)
 
250
                        s=zin(2,j,nin2)
 
251
                        r2=(r-s)*rt2i
 
252
                        s2=(r+s)*rt2i
 
253
                        r3=zin(2,j,nin3)
 
254
                        s3=zin(1,j,nin3)
 
255
                        r=zin(1,j,nin4)
 
256
                        s=zin(2,j,nin4)
 
257
                        r4=(r + s)*rt2i
 
258
                        s4=(r - s)*rt2i
 
259
                        r=r1 - r3
 
260
                        s=r2 - r4
 
261
                        zout(1,j,nout1) = r + s
 
262
                        zout(1,j,nout3) = r - s
 
263
                        r=r1 + r3
 
264
                        s=s2 - s4
 
265
                        zout(1,j,nout2) = r - s
 
266
                        zout(1,j,nout4) = r + s
 
267
                        r=s1 + s3
 
268
                        s=s2 + s4
 
269
                        zout(2,j,nout1) = r + s
 
270
                        zout(2,j,nout3) = r - s
 
271
                        r=s1 - s3
 
272
                        s=r2 + r4
 
273
                        zout(2,j,nout2) = r + s
 
274
                        zout(2,j,nout4) = r - s
 
275
4010                        continue
 
276
                else
 
277
                        itt=ias*before
 
278
                        itrig=itt+1
 
279
                        cr2=trig(1,itrig)
 
280
                        ci2=trig(2,itrig)
 
281
                        itrig=itrig+itt
 
282
                        cr3=trig(1,itrig)
 
283
                        ci3=trig(2,itrig)
 
284
                        itrig=itrig+itt
 
285
                        cr4=trig(1,itrig)
 
286
                        ci4=trig(2,itrig)
 
287
                        nin1=ia-after
 
288
                        nout1=ia-atn
 
289
                        do 4020,ib=1,before
 
290
                        nin1=nin1+after
 
291
                        nin2=nin1+atb
 
292
                        nin3=nin2+atb
 
293
                        nin4=nin3+atb
 
294
                        nout1=nout1+atn
 
295
                        nout2=nout1+after
 
296
                        nout3=nout2+after
 
297
                        nout4=nout3+after
 
298
                        do 4020,j=1,n1dfft
 
299
                        r1=zin(1,j,nin1)
 
300
                        s1=zin(2,j,nin1)
 
301
                        r=zin(1,j,nin2)
 
302
                        s=zin(2,j,nin2)
 
303
                        r2=r*cr2 - s*ci2
 
304
                        s2=r*ci2 + s*cr2
 
305
                        r=zin(1,j,nin3)
 
306
                        s=zin(2,j,nin3)
 
307
                        r3=r*cr3 - s*ci3
 
308
                        s3=r*ci3 + s*cr3
 
309
                        r=zin(1,j,nin4)
 
310
                        s=zin(2,j,nin4)
 
311
                        r4=r*cr4 - s*ci4
 
312
                        s4=r*ci4 + s*cr4
 
313
                        r=r1 + r3
 
314
                        s=r2 + r4
 
315
                        zout(1,j,nout1) = r + s
 
316
                        zout(1,j,nout3) = r - s
 
317
                        r=r1 - r3
 
318
                        s=s2 - s4
 
319
                        zout(1,j,nout2) = r - s
 
320
                        zout(1,j,nout4) = r + s
 
321
                        r=s1 + s3
 
322
                        s=s2 + s4
 
323
                        zout(2,j,nout1) = r + s
 
324
                        zout(2,j,nout3) = r - s
 
325
                        r=s1 - s3
 
326
                        s=r2 - r4
 
327
                        zout(2,j,nout2) = r + s
 
328
                        zout(2,j,nout4) = r - s
 
329
4020                        continue
 
330
                end if
 
331
4000                continue
 
332
        else
 
333
                ia=1
 
334
                nin1=ia-after
 
335
                nout1=ia-atn
 
336
                do 4101,ib=1,before
 
337
                nin1=nin1+after
 
338
                nin2=nin1+atb
 
339
                nin3=nin2+atb
 
340
                nin4=nin3+atb
 
341
                nout1=nout1+atn
 
342
                nout2=nout1+after
 
343
                nout3=nout2+after
 
344
                nout4=nout3+after
 
345
                do 4101,j=1,n1dfft
 
346
                r1=zin(1,j,nin1)
 
347
                s1=zin(2,j,nin1)
 
348
                r2=zin(1,j,nin2)
 
349
                s2=zin(2,j,nin2)
 
350
                r3=zin(1,j,nin3)
 
351
                s3=zin(2,j,nin3)
 
352
                r4=zin(1,j,nin4)
 
353
                s4=zin(2,j,nin4)
 
354
                r=r1 + r3
 
355
                s=r2 + r4
 
356
                zout(1,j,nout1) = r + s
 
357
                zout(1,j,nout3) = r - s
 
358
                r=r1 - r3
 
359
                s=s2 - s4
 
360
                zout(1,j,nout2) = r + s
 
361
                zout(1,j,nout4) = r - s
 
362
                r=s1 + s3
 
363
                s=s2 + s4
 
364
                zout(2,j,nout1) = r + s
 
365
                zout(2,j,nout3) = r - s
 
366
                r=s1 - s3
 
367
                s=r2 - r4
 
368
                zout(2,j,nout2) = r - s
 
369
                zout(2,j,nout4) = r + s
 
370
4101                continue
 
371
                do 4100,ia=2,after
 
372
                ias=ia-1
 
373
                if (2*ias.eq.after) then
 
374
                        nin1=ia-after
 
375
                        nout1=ia-atn
 
376
                        do 4110,ib=1,before
 
377
                        nin1=nin1+after
 
378
                        nin2=nin1+atb
 
379
                        nin3=nin2+atb
 
380
                        nin4=nin3+atb
 
381
                        nout1=nout1+atn
 
382
                        nout2=nout1+after
 
383
                        nout3=nout2+after
 
384
                        nout4=nout3+after
 
385
                        do 4110,j=1,n1dfft
 
386
                        r1=zin(1,j,nin1)
 
387
                        s1=zin(2,j,nin1)
 
388
                        r=zin(1,j,nin2)
 
389
                        s=zin(2,j,nin2)
 
390
                        r2=(r + s)*rt2i
 
391
                        s2=(s - r)*rt2i
 
392
                        r3=zin(2,j,nin3)
 
393
                        s3=zin(1,j,nin3)
 
394
                        r=zin(1,j,nin4)
 
395
                        s=zin(2,j,nin4)
 
396
                        r4=(s - r)*rt2i
 
397
                        s4=(r + s)*rt2i
 
398
                        r=r1 + r3
 
399
                        s=r2 + r4
 
400
                        zout(1,j,nout1) = r + s
 
401
                        zout(1,j,nout3) = r - s
 
402
                        r=r1 - r3
 
403
                        s=s2 + s4
 
404
                        zout(1,j,nout2) = r + s
 
405
                        zout(1,j,nout4) = r - s
 
406
                        r=s1 - s3
 
407
                        s=s2 - s4
 
408
                        zout(2,j,nout1) = r + s
 
409
                        zout(2,j,nout3) = r - s
 
410
                        r=s1 + s3
 
411
                        s=r2 - r4
 
412
                        zout(2,j,nout2) = r - s
 
413
                        zout(2,j,nout4) = r + s
 
414
4110                        continue
 
415
                else
 
416
                        itt=ias*before
 
417
                        itrig=itt+1
 
418
                        cr2=trig(1,itrig)
 
419
                        ci2=trig(2,itrig)
 
420
                        itrig=itrig+itt
 
421
                        cr3=trig(1,itrig)
 
422
                        ci3=trig(2,itrig)
 
423
                        itrig=itrig+itt
 
424
                        cr4=trig(1,itrig)
 
425
                        ci4=trig(2,itrig)
 
426
                        nin1=ia-after
 
427
                        nout1=ia-atn
 
428
                        do 4120,ib=1,before
 
429
                        nin1=nin1+after
 
430
                        nin2=nin1+atb
 
431
                        nin3=nin2+atb
 
432
                        nin4=nin3+atb
 
433
                        nout1=nout1+atn
 
434
                        nout2=nout1+after
 
435
                        nout3=nout2+after
 
436
                        nout4=nout3+after
 
437
                        do 4120,j=1,n1dfft
 
438
                        r1=zin(1,j,nin1)
 
439
                        s1=zin(2,j,nin1)
 
440
                        r=zin(1,j,nin2)
 
441
                        s=zin(2,j,nin2)
 
442
                        r2=r*cr2 - s*ci2
 
443
                        s2=r*ci2 + s*cr2
 
444
                        r=zin(1,j,nin3)
 
445
                        s=zin(2,j,nin3)
 
446
                        r3=r*cr3 - s*ci3
 
447
                        s3=r*ci3 + s*cr3
 
448
                        r=zin(1,j,nin4)
 
449
                        s=zin(2,j,nin4)
 
450
                        r4=r*cr4 - s*ci4
 
451
                        s4=r*ci4 + s*cr4
 
452
                        r=r1 + r3
 
453
                        s=r2 + r4
 
454
                        zout(1,j,nout1) = r + s
 
455
                        zout(1,j,nout3) = r - s
 
456
                        r=r1 - r3
 
457
                        s=s2 - s4
 
458
                        zout(1,j,nout2) = r + s
 
459
                        zout(1,j,nout4) = r - s
 
460
                        r=s1 + s3
 
461
                        s=s2 + s4
 
462
                        zout(2,j,nout1) = r + s
 
463
                        zout(2,j,nout3) = r - s
 
464
                        r=s1 - s3
 
465
                        s=r2 - r4
 
466
                        zout(2,j,nout2) = r - s
 
467
                        zout(2,j,nout4) = r + s
 
468
4120                        continue
 
469
                end if
 
470
4100                continue
 
471
        end if
 
472
        else if (now.eq.8) then
 
473
        if (isign.eq.-1) then
 
474
                ia=1
 
475
                        nin1=ia-after
 
476
                        nout1=ia-atn
 
477
                        do 8120,ib=1,before
 
478
                        nin1=nin1+after
 
479
                        nin2=nin1+atb
 
480
                        nin3=nin2+atb
 
481
                        nin4=nin3+atb
 
482
                        nin5=nin4+atb
 
483
                        nin6=nin5+atb
 
484
                        nin7=nin6+atb
 
485
                        nin8=nin7+atb
 
486
                        nout1=nout1+atn
 
487
                        nout2=nout1+after
 
488
                        nout3=nout2+after
 
489
                        nout4=nout3+after
 
490
                        nout5=nout4+after
 
491
                        nout6=nout5+after
 
492
                        nout7=nout6+after
 
493
                        nout8=nout7+after
 
494
                        do 8120,j=1,n1dfft
 
495
                        r1=zin(1,j,nin1)
 
496
                        s1=zin(2,j,nin1)
 
497
                        r2=zin(1,j,nin2)
 
498
                        s2=zin(2,j,nin2)
 
499
                        r3=zin(1,j,nin3)
 
500
                        s3=zin(2,j,nin3)
 
501
                        r4=zin(1,j,nin4)
 
502
                        s4=zin(2,j,nin4)
 
503
                        r5=zin(1,j,nin5)
 
504
                        s5=zin(2,j,nin5)
 
505
                        r6=zin(1,j,nin6)
 
506
                        s6=zin(2,j,nin6)
 
507
                        r7=zin(1,j,nin7)
 
508
                        s7=zin(2,j,nin7)
 
509
                        r8=zin(1,j,nin8)
 
510
                        s8=zin(2,j,nin8)
 
511
                        r=r1 + r5
 
512
                        s=r3 + r7
 
513
                        ap=r + s
 
514
                        am=r - s
 
515
                        r=r2 + r6
 
516
                        s=r4 + r8
 
517
                        bp=r + s
 
518
                        bm=r - s
 
519
                        r=s1 + s5
 
520
                        s=s3 + s7
 
521
                        cp=r + s
 
522
                        cm=r - s
 
523
                        r=s2 + s6
 
524
                        s=s4 + s8
 
525
                        dpp=r + s
 
526
                        dm=r - s
 
527
                        zout(1,j,nout1) = ap + bp
 
528
                        zout(2,j,nout1) = cp + dpp
 
529
                        zout(1,j,nout5) = ap - bp
 
530
                        zout(2,j,nout5) = cp - dpp
 
531
                        zout(1,j,nout3) = am + dm
 
532
                        zout(2,j,nout3) = cm - bm
 
533
                        zout(1,j,nout7) = am - dm
 
534
                        zout(2,j,nout7) = cm + bm
 
535
                        r=r1 - r5
 
536
                        s=s3 - s7
 
537
                        ap=r + s
 
538
                        am=r - s
 
539
                        r=s1 - s5
 
540
                        s=r3 - r7
 
541
                        bp=r + s
 
542
                        bm=r - s
 
543
                        r=s4 - s8
 
544
                        s=r2 - r6
 
545
                        cp=r + s
 
546
                        cm=r - s
 
547
                        r=s2 - s6
 
548
                        s=r4 - r8
 
549
                        dpp=r + s
 
550
                        dm=r - s
 
551
                        r = ( cp + dm)*rt2i
 
552
                        s = ( dm - cp)*rt2i
 
553
                        cp= ( cm + dpp)*rt2i
 
554
                        dpp = ( cm - dpp)*rt2i
 
555
                        zout(1,j,nout2) = ap + r
 
556
                        zout(2,j,nout2) = bm + s
 
557
                        zout(1,j,nout6) = ap - r
 
558
                        zout(2,j,nout6) = bm - s
 
559
                        zout(1,j,nout4) = am + cp
 
560
                        zout(2,j,nout4) = bp + dpp
 
561
                        zout(1,j,nout8) = am - cp
 
562
                        zout(2,j,nout8) = bp - dpp
 
563
8120                        continue
 
564
                do 8000,ia=2,after
 
565
                ias=ia-1
 
566
                        itt=ias*before
 
567
                        itrig=itt+1
 
568
                        cr2=trig(1,itrig)
 
569
                        ci2=trig(2,itrig)
 
570
                        itrig=itrig+itt
 
571
                        cr3=trig(1,itrig)
 
572
                        ci3=trig(2,itrig)
 
573
                        itrig=itrig+itt
 
574
                        cr4=trig(1,itrig)
 
575
                        ci4=trig(2,itrig)
 
576
                        itrig=itrig+itt
 
577
                        cr5=trig(1,itrig)
 
578
                        ci5=trig(2,itrig)
 
579
                        itrig=itrig+itt
 
580
                        cr6=trig(1,itrig)
 
581
                        ci6=trig(2,itrig)
 
582
                        itrig=itrig+itt
 
583
                        cr7=trig(1,itrig)
 
584
                        ci7=trig(2,itrig)
 
585
                        itrig=itrig+itt
 
586
                        cr8=trig(1,itrig)
 
587
                        ci8=trig(2,itrig)
 
588
                        nin1=ia-after
 
589
                        nout1=ia-atn
 
590
                        do 8020,ib=1,before
 
591
                        nin1=nin1+after
 
592
                        nin2=nin1+atb
 
593
                        nin3=nin2+atb
 
594
                        nin4=nin3+atb
 
595
                        nin5=nin4+atb
 
596
                        nin6=nin5+atb
 
597
                        nin7=nin6+atb
 
598
                        nin8=nin7+atb
 
599
                        nout1=nout1+atn
 
600
                        nout2=nout1+after
 
601
                        nout3=nout2+after
 
602
                        nout4=nout3+after
 
603
                        nout5=nout4+after
 
604
                        nout6=nout5+after
 
605
                        nout7=nout6+after
 
606
                        nout8=nout7+after
 
607
                        do 8020,j=1,n1dfft
 
608
                        r1=zin(1,j,nin1)
 
609
                        s1=zin(2,j,nin1)
 
610
                        r=zin(1,j,nin2)
 
611
                        s=zin(2,j,nin2)
 
612
                        r2=r*cr2 - s*ci2
 
613
                        s2=r*ci2 + s*cr2
 
614
                        r=zin(1,j,nin3)
 
615
                        s=zin(2,j,nin3)
 
616
                        r3=r*cr3 - s*ci3
 
617
                        s3=r*ci3 + s*cr3
 
618
                        r=zin(1,j,nin4)
 
619
                        s=zin(2,j,nin4)
 
620
                        r4=r*cr4 - s*ci4
 
621
                        s4=r*ci4 + s*cr4
 
622
                        r=zin(1,j,nin5)
 
623
                        s=zin(2,j,nin5)
 
624
                        r5=r*cr5 - s*ci5
 
625
                        s5=r*ci5 + s*cr5
 
626
                        r=zin(1,j,nin6)
 
627
                        s=zin(2,j,nin6)
 
628
                        r6=r*cr6 - s*ci6
 
629
                        s6=r*ci6 + s*cr6
 
630
                        r=zin(1,j,nin7)
 
631
                        s=zin(2,j,nin7)
 
632
                        r7=r*cr7 - s*ci7
 
633
                        s7=r*ci7 + s*cr7
 
634
                        r=zin(1,j,nin8)
 
635
                        s=zin(2,j,nin8)
 
636
                        r8=r*cr8 - s*ci8
 
637
                        s8=r*ci8 + s*cr8
 
638
                        r=r1 + r5
 
639
                        s=r3 + r7
 
640
                        ap=r + s
 
641
                        am=r - s
 
642
                        r=r2 + r6
 
643
                        s=r4 + r8
 
644
                        bp=r + s
 
645
                        bm=r - s
 
646
                        r=s1 + s5
 
647
                        s=s3 + s7
 
648
                        cp=r + s
 
649
                        cm=r - s
 
650
                        r=s2 + s6
 
651
                        s=s4 + s8
 
652
                        dpp=r + s
 
653
                        dm=r - s
 
654
                        zout(1,j,nout1) = ap + bp
 
655
                        zout(2,j,nout1) = cp + dpp
 
656
                        zout(1,j,nout5) = ap - bp
 
657
                        zout(2,j,nout5) = cp - dpp
 
658
                        zout(1,j,nout3) = am + dm
 
659
                        zout(2,j,nout3) = cm - bm
 
660
                        zout(1,j,nout7) = am - dm
 
661
                        zout(2,j,nout7) = cm + bm
 
662
                        r=r1 - r5
 
663
                        s=s3 - s7
 
664
                        ap=r + s
 
665
                        am=r - s
 
666
                        r=s1 - s5
 
667
                        s=r3 - r7
 
668
                        bp=r + s
 
669
                        bm=r - s
 
670
                        r=s4 - s8
 
671
                        s=r2 - r6
 
672
                        cp=r + s
 
673
                        cm=r - s
 
674
                        r=s2 - s6
 
675
                        s=r4 - r8
 
676
                        dpp=r + s
 
677
                        dm=r - s
 
678
                        r = ( cp + dm)*rt2i
 
679
                        s = ( dm - cp)*rt2i
 
680
                        cp= ( cm + dpp)*rt2i
 
681
                        dpp = ( cm - dpp)*rt2i
 
682
                        zout(1,j,nout2) = ap + r
 
683
                        zout(2,j,nout2) = bm + s
 
684
                        zout(1,j,nout6) = ap - r
 
685
                        zout(2,j,nout6) = bm - s
 
686
                        zout(1,j,nout4) = am + cp
 
687
                        zout(2,j,nout4) = bp + dpp
 
688
                        zout(1,j,nout8) = am - cp
 
689
                        zout(2,j,nout8) = bp - dpp
 
690
 
 
691
8020                        continue
 
692
8000                continue
 
693
 
 
694
        else
 
695
                ia=1
 
696
                        nin1=ia-after
 
697
                        nout1=ia-atn
 
698
                        do 8121,ib=1,before
 
699
                        nin1=nin1+after
 
700
                        nin2=nin1+atb
 
701
                        nin3=nin2+atb
 
702
                        nin4=nin3+atb
 
703
                        nin5=nin4+atb
 
704
                        nin6=nin5+atb
 
705
                        nin7=nin6+atb
 
706
                        nin8=nin7+atb
 
707
                        nout1=nout1+atn
 
708
                        nout2=nout1+after
 
709
                        nout3=nout2+after
 
710
                        nout4=nout3+after
 
711
                        nout5=nout4+after
 
712
                        nout6=nout5+after
 
713
                        nout7=nout6+after
 
714
                        nout8=nout7+after
 
715
                        do 8121,j=1,n1dfft
 
716
                        r1=zin(1,j,nin1)
 
717
                        s1=zin(2,j,nin1)
 
718
                        r2=zin(1,j,nin2)
 
719
                        s2=zin(2,j,nin2)
 
720
                        r3=zin(1,j,nin3)
 
721
                        s3=zin(2,j,nin3)
 
722
                        r4=zin(1,j,nin4)
 
723
                        s4=zin(2,j,nin4)
 
724
                        r5=zin(1,j,nin5)
 
725
                        s5=zin(2,j,nin5)
 
726
                        r6=zin(1,j,nin6)
 
727
                        s6=zin(2,j,nin6)
 
728
                        r7=zin(1,j,nin7)
 
729
                        s7=zin(2,j,nin7)
 
730
                        r8=zin(1,j,nin8)
 
731
                        s8=zin(2,j,nin8)
 
732
                        r=r1 + r5
 
733
                        s=r3 + r7
 
734
                        ap=r + s
 
735
                        am=r - s
 
736
                        r=r2 + r6
 
737
                        s=r4 + r8
 
738
                        bp=r + s
 
739
                        bm=r - s
 
740
                        r=s1 + s5
 
741
                        s=s3 + s7
 
742
                        cp=r + s
 
743
                        cm=r - s
 
744
                        r=s2 + s6
 
745
                        s=s4 + s8
 
746
                        dpp=r + s
 
747
                        dm=r - s
 
748
                        zout(1,j,nout1) = ap + bp
 
749
                        zout(2,j,nout1) = cp + dpp
 
750
                        zout(1,j,nout5) = ap - bp
 
751
                        zout(2,j,nout5) = cp - dpp
 
752
                        zout(1,j,nout3) = am - dm
 
753
                        zout(2,j,nout3) = cm + bm
 
754
                        zout(1,j,nout7) = am + dm
 
755
                        zout(2,j,nout7) = cm - bm
 
756
                        r= r1 - r5
 
757
                        s=-s3 + s7
 
758
                        ap=r + s
 
759
                        am=r - s
 
760
                        r=s1 - s5
 
761
                        s=r7 - r3
 
762
                        bp=r + s
 
763
                        bm=r - s
 
764
                        r=-s4 + s8
 
765
                        s= r2 - r6
 
766
                        cp=r + s
 
767
                        cm=r - s
 
768
                        r=-s2 + s6
 
769
                        s= r4 - r8
 
770
                        dpp=r + s
 
771
                        dm=r - s
 
772
                        r = ( cp + dm)*rt2i
 
773
                        s = ( cp - dm)*rt2i
 
774
                        cp= ( cm + dpp)*rt2i
 
775
                        dpp= ( dpp - cm)*rt2i
 
776
                        zout(1,j,nout2) = ap + r
 
777
                        zout(2,j,nout2) = bm + s
 
778
                        zout(1,j,nout6) = ap - r
 
779
                        zout(2,j,nout6) = bm - s
 
780
                        zout(1,j,nout4) = am + cp
 
781
                        zout(2,j,nout4) = bp + dpp
 
782
                        zout(1,j,nout8) = am - cp
 
783
                        zout(2,j,nout8) = bp - dpp
 
784
8121                        continue
 
785
 
 
786
                do 8001,ia=2,after
 
787
                ias=ia-1
 
788
                        itt=ias*before
 
789
                        itrig=itt+1
 
790
                        cr2=trig(1,itrig)
 
791
                        ci2=trig(2,itrig)
 
792
                        itrig=itrig+itt
 
793
                        cr3=trig(1,itrig)
 
794
                        ci3=trig(2,itrig)
 
795
                        itrig=itrig+itt
 
796
                        cr4=trig(1,itrig)
 
797
                        ci4=trig(2,itrig)
 
798
                        itrig=itrig+itt
 
799
                        cr5=trig(1,itrig)
 
800
                        ci5=trig(2,itrig)
 
801
                        itrig=itrig+itt
 
802
                        cr6=trig(1,itrig)
 
803
                        ci6=trig(2,itrig)
 
804
                        itrig=itrig+itt
 
805
                        cr7=trig(1,itrig)
 
806
                        ci7=trig(2,itrig)
 
807
                        itrig=itrig+itt
 
808
                        cr8=trig(1,itrig)
 
809
                        ci8=trig(2,itrig)
 
810
                        nin1=ia-after
 
811
                        nout1=ia-atn
 
812
                        do 8021,ib=1,before
 
813
                        nin1=nin1+after
 
814
                        nin2=nin1+atb
 
815
                        nin3=nin2+atb
 
816
                        nin4=nin3+atb
 
817
                        nin5=nin4+atb
 
818
                        nin6=nin5+atb
 
819
                        nin7=nin6+atb
 
820
                        nin8=nin7+atb
 
821
                        nout1=nout1+atn
 
822
                        nout2=nout1+after
 
823
                        nout3=nout2+after
 
824
                        nout4=nout3+after
 
825
                        nout5=nout4+after
 
826
                        nout6=nout5+after
 
827
                        nout7=nout6+after
 
828
                        nout8=nout7+after
 
829
                        do 8021,j=1,n1dfft
 
830
                        r1=zin(1,j,nin1)
 
831
                        s1=zin(2,j,nin1)
 
832
                        r=zin(1,j,nin2)
 
833
                        s=zin(2,j,nin2)
 
834
                        r2=r*cr2 - s*ci2
 
835
                        s2=r*ci2 + s*cr2
 
836
                        r=zin(1,j,nin3)
 
837
                        s=zin(2,j,nin3)
 
838
                        r3=r*cr3 - s*ci3
 
839
                        s3=r*ci3 + s*cr3
 
840
                        r=zin(1,j,nin4)
 
841
                        s=zin(2,j,nin4)
 
842
                        r4=r*cr4 - s*ci4
 
843
                        s4=r*ci4 + s*cr4
 
844
                        r=zin(1,j,nin5)
 
845
                        s=zin(2,j,nin5)
 
846
                        r5=r*cr5 - s*ci5
 
847
                        s5=r*ci5 + s*cr5
 
848
                        r=zin(1,j,nin6)
 
849
                        s=zin(2,j,nin6)
 
850
                        r6=r*cr6 - s*ci6
 
851
                        s6=r*ci6 + s*cr6
 
852
                        r=zin(1,j,nin7)
 
853
                        s=zin(2,j,nin7)
 
854
                        r7=r*cr7 - s*ci7
 
855
                        s7=r*ci7 + s*cr7
 
856
                        r=zin(1,j,nin8)
 
857
                        s=zin(2,j,nin8)
 
858
                        r8=r*cr8 - s*ci8
 
859
                        s8=r*ci8 + s*cr8
 
860
                        r=r1 + r5
 
861
                        s=r3 + r7
 
862
                        ap=r + s
 
863
                        am=r - s
 
864
                        r=r2 + r6
 
865
                        s=r4 + r8
 
866
                        bp=r + s
 
867
                        bm=r - s
 
868
                        r=s1 + s5
 
869
                        s=s3 + s7
 
870
                        cp=r + s
 
871
                        cm=r - s
 
872
                        r=s2 + s6
 
873
                        s=s4 + s8
 
874
                        dpp=r + s
 
875
                        dm=r - s
 
876
                        zout(1,j,nout1) = ap + bp
 
877
                        zout(2,j,nout1) = cp + dpp
 
878
                        zout(1,j,nout5) = ap - bp
 
879
                        zout(2,j,nout5) = cp - dpp
 
880
                        zout(1,j,nout3) = am - dm
 
881
                        zout(2,j,nout3) = cm + bm
 
882
                        zout(1,j,nout7) = am + dm
 
883
                        zout(2,j,nout7) = cm - bm
 
884
                        r= r1 - r5
 
885
                        s=-s3 + s7
 
886
                        ap=r + s
 
887
                        am=r - s
 
888
                        r=s1 - s5
 
889
                        s=r7 - r3
 
890
                        bp=r + s
 
891
                        bm=r - s
 
892
                        r=-s4 + s8
 
893
                        s= r2 - r6
 
894
                        cp=r + s
 
895
                        cm=r - s
 
896
                        r=-s2 + s6
 
897
                        s= r4 - r8
 
898
                        dpp=r + s
 
899
                        dm=r - s
 
900
                        r = ( cp + dm)*rt2i
 
901
                        s = ( cp - dm)*rt2i
 
902
                        cp= ( cm + dpp)*rt2i
 
903
                        dpp= ( dpp - cm)*rt2i
 
904
                        zout(1,j,nout2) = ap + r
 
905
                        zout(2,j,nout2) = bm + s
 
906
                        zout(1,j,nout6) = ap - r
 
907
                        zout(2,j,nout6) = bm - s
 
908
                        zout(1,j,nout4) = am + cp
 
909
                        zout(2,j,nout4) = bp + dpp
 
910
                        zout(1,j,nout8) = am - cp
 
911
                        zout(2,j,nout8) = bp - dpp
 
912
8021                        continue
 
913
8001                continue
 
914
 
 
915
        end if
 
916
        else if (now.eq.3) then
 
917
!         .5d0*sqrt(3.d0)
 
918
        bb=isign*0.8660254037844387d0
 
919
        ia=1
 
920
        nin1=ia-after
 
921
        nout1=ia-atn
 
922
        do 3001,ib=1,before
 
923
        nin1=nin1+after
 
924
        nin2=nin1+atb
 
925
        nin3=nin2+atb
 
926
        nout1=nout1+atn
 
927
        nout2=nout1+after
 
928
        nout3=nout2+after
 
929
        do 3001,j=1,n1dfft
 
930
        r1=zin(1,j,nin1)
 
931
        s1=zin(2,j,nin1)
 
932
        r2=zin(1,j,nin2)
 
933
        s2=zin(2,j,nin2)
 
934
        r3=zin(1,j,nin3)
 
935
        s3=zin(2,j,nin3)
 
936
        r=r2 + r3
 
937
        s=s2 + s3
 
938
        zout(1,j,nout1) = r + r1
 
939
        zout(2,j,nout1) = s + s1
 
940
        r1=r1 - .5d0*r
 
941
        s1=s1 - .5d0*s
 
942
        r2=bb*(r2-r3)
 
943
        s2=bb*(s2-s3)
 
944
        zout(1,j,nout2) = r1 - s2
 
945
        zout(2,j,nout2) = s1 + r2
 
946
        zout(1,j,nout3) = r1 + s2
 
947
        zout(2,j,nout3) = s1 - r2
 
948
3001        continue
 
949
        do 3000,ia=2,after
 
950
        ias=ia-1
 
951
        if (4*ias.eq.3*after) then
 
952
        if (isign.eq.1) then
 
953
                nin1=ia-after
 
954
                nout1=ia-atn
 
955
                do 3010,ib=1,before
 
956
                nin1=nin1+after
 
957
                nin2=nin1+atb
 
958
                nin3=nin2+atb
 
959
                nout1=nout1+atn
 
960
                nout2=nout1+after
 
961
                nout3=nout2+after
 
962
                do 3010,j=1,n1dfft
 
963
                r1=zin(1,j,nin1)
 
964
                s1=zin(2,j,nin1)
 
965
                r2=zin(2,j,nin2)
 
966
                s2=zin(1,j,nin2)
 
967
                r3=zin(1,j,nin3)
 
968
                s3=zin(2,j,nin3)
 
969
                r=r3 + r2
 
970
                s=s2 - s3
 
971
                zout(1,j,nout1) = r1 - r
 
972
                zout(2,j,nout1) = s + s1
 
973
                r1=r1 + .5d0*r
 
974
                s1=s1 - .5d0*s
 
975
                r2=bb*(r2-r3)
 
976
                s2=bb*(s2+s3)
 
977
                zout(1,j,nout2) = r1 - s2
 
978
                zout(2,j,nout2) = s1 - r2
 
979
                zout(1,j,nout3) = r1 + s2
 
980
                zout(2,j,nout3) = s1 + r2
 
981
3010                continue
 
982
        else
 
983
                nin1=ia-after
 
984
                nout1=ia-atn
 
985
                do 3020,ib=1,before
 
986
                nin1=nin1+after
 
987
                nin2=nin1+atb
 
988
                nin3=nin2+atb
 
989
                nout1=nout1+atn
 
990
                nout2=nout1+after
 
991
                nout3=nout2+after
 
992
                do 3020,j=1,n1dfft
 
993
                r1=zin(1,j,nin1)
 
994
                s1=zin(2,j,nin1)
 
995
                r2=zin(2,j,nin2)
 
996
                s2=zin(1,j,nin2)
 
997
                r3=zin(1,j,nin3)
 
998
                s3=zin(2,j,nin3)
 
999
                r=r2 - r3
 
1000
                s=s2 + s3
 
1001
                zout(1,j,nout1) = r + r1
 
1002
                zout(2,j,nout1) = s1 - s
 
1003
                r1=r1 - .5d0*r
 
1004
                s1=s1 + .5d0*s
 
1005
                r2=bb*(r2+r3)
 
1006
                s2=bb*(s2-s3)
 
1007
                zout(1,j,nout2) = r1 + s2
 
1008
                zout(2,j,nout2) = s1 + r2
 
1009
                zout(1,j,nout3) = r1 - s2
 
1010
                zout(2,j,nout3) = s1 - r2
 
1011
3020                continue
 
1012
        end if
 
1013
        else if (8*ias.eq.3*after) then
 
1014
        if (isign.eq.1) then
 
1015
                nin1=ia-after
 
1016
                nout1=ia-atn
 
1017
                do 3030,ib=1,before
 
1018
                nin1=nin1+after
 
1019
                nin2=nin1+atb
 
1020
                nin3=nin2+atb
 
1021
                nout1=nout1+atn
 
1022
                nout2=nout1+after
 
1023
                nout3=nout2+after
 
1024
                do 3030,j=1,n1dfft
 
1025
                r1=zin(1,j,nin1)
 
1026
                s1=zin(2,j,nin1)
 
1027
                r=zin(1,j,nin2)
 
1028
                s=zin(2,j,nin2)
 
1029
                r2=(r - s)*rt2i
 
1030
                s2=(r + s)*rt2i
 
1031
                r3=zin(2,j,nin3)
 
1032
                s3=zin(1,j,nin3)
 
1033
                r=r2 - r3
 
1034
                s=s2 + s3
 
1035
                zout(1,j,nout1) = r + r1
 
1036
                zout(2,j,nout1) = s + s1
 
1037
                r1=r1 - .5d0*r
 
1038
                s1=s1 - .5d0*s
 
1039
                r2=bb*(r2+r3)
 
1040
                s2=bb*(s2-s3)
 
1041
                zout(1,j,nout2) = r1 - s2
 
1042
                zout(2,j,nout2) = s1 + r2
 
1043
                zout(1,j,nout3) = r1 + s2
 
1044
                zout(2,j,nout3) = s1 - r2
 
1045
3030                continue
 
1046
        else
 
1047
                nin1=ia-after
 
1048
                nout1=ia-atn
 
1049
                do 3040,ib=1,before
 
1050
                nin1=nin1+after
 
1051
                nin2=nin1+atb
 
1052
                nin3=nin2+atb
 
1053
                nout1=nout1+atn
 
1054
                nout2=nout1+after
 
1055
                nout3=nout2+after
 
1056
                do 3040,j=1,n1dfft
 
1057
                r1=zin(1,j,nin1)
 
1058
                s1=zin(2,j,nin1)
 
1059
                r=zin(1,j,nin2)
 
1060
                s=zin(2,j,nin2)
 
1061
                r2=(r + s)*rt2i
 
1062
                s2=(s - r)*rt2i
 
1063
                r3=zin(2,j,nin3)
 
1064
                s3=zin(1,j,nin3)
 
1065
                r=r2 + r3
 
1066
                s=s2 - s3
 
1067
                zout(1,j,nout1) = r + r1
 
1068
                zout(2,j,nout1) = s + s1
 
1069
                r1=r1 - .5d0*r
 
1070
                s1=s1 - .5d0*s
 
1071
                r2=bb*(r2-r3)
 
1072
                s2=bb*(s2+s3)
 
1073
                zout(1,j,nout2) = r1 - s2
 
1074
                zout(2,j,nout2) = s1 + r2
 
1075
                zout(1,j,nout3) = r1 + s2
 
1076
                zout(2,j,nout3) = s1 - r2
 
1077
3040                continue
 
1078
        end if
 
1079
        else
 
1080
        itt=ias*before
 
1081
        itrig=itt+1
 
1082
        cr2=trig(1,itrig)
 
1083
        ci2=trig(2,itrig)
 
1084
        itrig=itrig+itt
 
1085
        cr3=trig(1,itrig)
 
1086
        ci3=trig(2,itrig)
 
1087
        nin1=ia-after
 
1088
        nout1=ia-atn
 
1089
        do 3090,ib=1,before
 
1090
        nin1=nin1+after
 
1091
        nin2=nin1+atb
 
1092
        nin3=nin2+atb
 
1093
        nout1=nout1+atn
 
1094
        nout2=nout1+after
 
1095
        nout3=nout2+after
 
1096
        do 3090,j=1,n1dfft
 
1097
        r1=zin(1,j,nin1)
 
1098
        s1=zin(2,j,nin1)
 
1099
        r=zin(1,j,nin2)
 
1100
        s=zin(2,j,nin2)
 
1101
        r2=r*cr2 - s*ci2
 
1102
        s2=r*ci2 + s*cr2
 
1103
        r=zin(1,j,nin3)
 
1104
        s=zin(2,j,nin3)
 
1105
        r3=r*cr3 - s*ci3
 
1106
        s3=r*ci3 + s*cr3
 
1107
        r=r2 + r3
 
1108
        s=s2 + s3
 
1109
        zout(1,j,nout1) = r + r1
 
1110
        zout(2,j,nout1) = s + s1
 
1111
        r1=r1 - .5d0*r
 
1112
        s1=s1 - .5d0*s
 
1113
        r2=bb*(r2-r3)
 
1114
        s2=bb*(s2-s3)
 
1115
        zout(1,j,nout2) = r1 - s2
 
1116
        zout(2,j,nout2) = s1 + r2
 
1117
        zout(1,j,nout3) = r1 + s2
 
1118
        zout(2,j,nout3) = s1 - r2
 
1119
3090        continue
 
1120
        end if
 
1121
3000        continue
 
1122
        else if (now.eq.5) then
 
1123
!         cos(2.d0*pi/5.d0)
 
1124
        cos2=0.3090169943749474d0
 
1125
!         cos(4.d0*pi/5.d0)
 
1126
        cos4=-0.8090169943749474d0
 
1127
!        sin(2.d0*pi/5.d0)
 
1128
        sin2=isign*0.9510565162951536d0
 
1129
!         sin(4.d0*pi/5.d0)
 
1130
        sin4=isign*0.5877852522924731d0
 
1131
        ia=1
 
1132
        nin1=ia-after
 
1133
        nout1=ia-atn
 
1134
        do 5001,ib=1,before
 
1135
        nin1=nin1+after
 
1136
        nin2=nin1+atb
 
1137
        nin3=nin2+atb
 
1138
        nin4=nin3+atb
 
1139
        nin5=nin4+atb
 
1140
        nout1=nout1+atn
 
1141
        nout2=nout1+after
 
1142
        nout3=nout2+after
 
1143
        nout4=nout3+after
 
1144
        nout5=nout4+after
 
1145
        do 5001,j=1,n1dfft
 
1146
        r1=zin(1,j,nin1)
 
1147
        s1=zin(2,j,nin1)
 
1148
        r2=zin(1,j,nin2)
 
1149
        s2=zin(2,j,nin2)
 
1150
        r3=zin(1,j,nin3)
 
1151
        s3=zin(2,j,nin3)
 
1152
        r4=zin(1,j,nin4)
 
1153
        s4=zin(2,j,nin4)
 
1154
        r5=zin(1,j,nin5)
 
1155
        s5=zin(2,j,nin5)
 
1156
        r25 = r2 + r5
 
1157
        r34 = r3 + r4
 
1158
        s25 = s2 - s5
 
1159
        s34 = s3 - s4
 
1160
        zout(1,j,nout1) = r1 + r25 + r34
 
1161
        r = r1 + cos2*r25 + cos4*r34
 
1162
        s = sin2*s25 + sin4*s34
 
1163
        zout(1,j,nout2) = r - s
 
1164
        zout(1,j,nout5) = r + s
 
1165
        r = r1 + cos4*r25 + cos2*r34
 
1166
        s = sin4*s25 - sin2*s34
 
1167
        zout(1,j,nout3) = r - s
 
1168
        zout(1,j,nout4) = r + s
 
1169
        r25 = r2 - r5
 
1170
        r34 = r3 - r4
 
1171
        s25 = s2 + s5
 
1172
        s34 = s3 + s4
 
1173
        zout(2,j,nout1) = s1 + s25 + s34
 
1174
        r = s1 + cos2*s25 + cos4*s34
 
1175
        s = sin2*r25 + sin4*r34
 
1176
        zout(2,j,nout2) = r + s
 
1177
        zout(2,j,nout5) = r - s
 
1178
        r = s1 + cos4*s25 + cos2*s34
 
1179
        s = sin4*r25 - sin2*r34
 
1180
        zout(2,j,nout3) = r + s
 
1181
        zout(2,j,nout4) = r - s
 
1182
5001        continue
 
1183
        do 5000,ia=2,after
 
1184
        ias=ia-1
 
1185
        if (8*ias.eq.5*after) then
 
1186
                if (isign.eq.1) then
 
1187
                        nin1=ia-after
 
1188
                        nout1=ia-atn
 
1189
                        do 5010,ib=1,before
 
1190
                        nin1=nin1+after
 
1191
                        nin2=nin1+atb
 
1192
                        nin3=nin2+atb
 
1193
                        nin4=nin3+atb
 
1194
                        nin5=nin4+atb
 
1195
                        nout1=nout1+atn
 
1196
                        nout2=nout1+after
 
1197
                        nout3=nout2+after
 
1198
                        nout4=nout3+after
 
1199
                        nout5=nout4+after
 
1200
                        do 5010,j=1,n1dfft
 
1201
                        r1=zin(1,j,nin1)
 
1202
                        s1=zin(2,j,nin1)
 
1203
                        r=zin(1,j,nin2)
 
1204
                        s=zin(2,j,nin2)
 
1205
                        r2=(r - s)*rt2i
 
1206
                        s2=(r + s)*rt2i
 
1207
                        r3=zin(2,j,nin3)
 
1208
                        s3=zin(1,j,nin3)
 
1209
                        r=zin(1,j,nin4)
 
1210
                        s=zin(2,j,nin4)
 
1211
                        r4=(r + s)*rt2i
 
1212
                        s4=(r - s)*rt2i
 
1213
                        r5=zin(1,j,nin5)
 
1214
                        s5=zin(2,j,nin5)
 
1215
                        r25 = r2 - r5
 
1216
                        r34 = r3 + r4
 
1217
                        s25 = s2 + s5
 
1218
                        s34 = s3 - s4
 
1219
                        zout(1,j,nout1) = r1 + r25 - r34
 
1220
                        r = r1 + cos2*r25 - cos4*r34
 
1221
                        s = sin2*s25 + sin4*s34
 
1222
                        zout(1,j,nout2) = r - s
 
1223
                        zout(1,j,nout5) = r + s
 
1224
                        r = r1 + cos4*r25 - cos2*r34
 
1225
                        s = sin4*s25 - sin2*s34
 
1226
                        zout(1,j,nout3) = r - s
 
1227
                        zout(1,j,nout4) = r + s
 
1228
                        r25 = r2 + r5
 
1229
                        r34 = r4 - r3
 
1230
                        s25 = s2 - s5
 
1231
                        s34 = s3 + s4
 
1232
                        zout(2,j,nout1) = s1 + s25 + s34
 
1233
                        r = s1 + cos2*s25 + cos4*s34
 
1234
                        s = sin2*r25 + sin4*r34
 
1235
                        zout(2,j,nout2) = r + s
 
1236
                        zout(2,j,nout5) = r - s
 
1237
                        r = s1 + cos4*s25 + cos2*s34
 
1238
                        s = sin4*r25 - sin2*r34
 
1239
                        zout(2,j,nout3) = r + s
 
1240
                        zout(2,j,nout4) = r - s
 
1241
5010                        continue
 
1242
                else
 
1243
                        nin1=ia-after
 
1244
                        nout1=ia-atn
 
1245
                        do 5020,ib=1,before
 
1246
                        nin1=nin1+after
 
1247
                        nin2=nin1+atb
 
1248
                        nin3=nin2+atb
 
1249
                        nin4=nin3+atb
 
1250
                        nin5=nin4+atb
 
1251
                        nout1=nout1+atn
 
1252
                        nout2=nout1+after
 
1253
                        nout3=nout2+after
 
1254
                        nout4=nout3+after
 
1255
                        nout5=nout4+after
 
1256
                        do 5020,j=1,n1dfft
 
1257
                        r1=zin(1,j,nin1)
 
1258
                        s1=zin(2,j,nin1)
 
1259
                        r=zin(1,j,nin2)
 
1260
                        s=zin(2,j,nin2)
 
1261
                        r2=(r + s)*rt2i
 
1262
                        s2=(s - r)*rt2i
 
1263
                        r3=zin(2,j,nin3)
 
1264
                        s3=zin(1,j,nin3)
 
1265
                        r=zin(1,j,nin4)
 
1266
                        s=zin(2,j,nin4)
 
1267
                        r4=(s - r)*rt2i
 
1268
                        s4=(r + s)*rt2i
 
1269
                        r5=zin(1,j,nin5)
 
1270
                        s5=zin(2,j,nin5)
 
1271
                        r25 = r2 - r5
 
1272
                        r34 = r3 + r4
 
1273
                        s25 = s2 + s5
 
1274
                        s34 = s4 - s3
 
1275
                        zout(1,j,nout1) = r1 + r25 + r34
 
1276
                        r = r1 + cos2*r25 + cos4*r34
 
1277
                        s = sin2*s25 + sin4*s34
 
1278
                        zout(1,j,nout2) = r - s
 
1279
                        zout(1,j,nout5) = r + s
 
1280
                        r = r1 + cos4*r25 + cos2*r34
 
1281
                        s = sin4*s25 - sin2*s34
 
1282
                        zout(1,j,nout3) = r - s
 
1283
                        zout(1,j,nout4) = r + s
 
1284
                        r25 = r2 + r5
 
1285
                        r34 = r3 - r4
 
1286
                        s25 = s2 - s5
 
1287
                        s34 = s3 + s4
 
1288
                        zout(2,j,nout1) = s1 + s25 - s34
 
1289
                        r = s1 + cos2*s25 - cos4*s34
 
1290
                        s = sin2*r25 + sin4*r34
 
1291
                        zout(2,j,nout2) = r + s
 
1292
                        zout(2,j,nout5) = r - s
 
1293
                        r = s1 + cos4*s25 - cos2*s34
 
1294
                        s = sin4*r25 - sin2*r34
 
1295
                        zout(2,j,nout3) = r + s
 
1296
                        zout(2,j,nout4) = r - s
 
1297
5020                        continue
 
1298
                end if
 
1299
        else
 
1300
                ias=ia-1
 
1301
                itt=ias*before
 
1302
                itrig=itt+1
 
1303
                cr2=trig(1,itrig)
 
1304
                ci2=trig(2,itrig)
 
1305
                itrig=itrig+itt
 
1306
                cr3=trig(1,itrig)
 
1307
                ci3=trig(2,itrig)
 
1308
                itrig=itrig+itt
 
1309
                cr4=trig(1,itrig)
 
1310
                ci4=trig(2,itrig)
 
1311
                itrig=itrig+itt
 
1312
                cr5=trig(1,itrig)
 
1313
                ci5=trig(2,itrig)
 
1314
                nin1=ia-after
 
1315
                nout1=ia-atn
 
1316
                do 5100,ib=1,before
 
1317
                nin1=nin1+after
 
1318
                nin2=nin1+atb
 
1319
                nin3=nin2+atb
 
1320
                nin4=nin3+atb
 
1321
                nin5=nin4+atb
 
1322
                nout1=nout1+atn
 
1323
                nout2=nout1+after
 
1324
                nout3=nout2+after
 
1325
                nout4=nout3+after
 
1326
                nout5=nout4+after
 
1327
                do 5100,j=1,n1dfft
 
1328
                r1=zin(1,j,nin1)
 
1329
                s1=zin(2,j,nin1)
 
1330
                r=zin(1,j,nin2)
 
1331
                s=zin(2,j,nin2)
 
1332
                r2=r*cr2 - s*ci2
 
1333
                s2=r*ci2 + s*cr2
 
1334
                r=zin(1,j,nin3)
 
1335
                s=zin(2,j,nin3)
 
1336
                r3=r*cr3 - s*ci3
 
1337
                s3=r*ci3 + s*cr3
 
1338
                r=zin(1,j,nin4)
 
1339
                s=zin(2,j,nin4)
 
1340
                r4=r*cr4 - s*ci4
 
1341
                s4=r*ci4 + s*cr4
 
1342
                r=zin(1,j,nin5)
 
1343
                s=zin(2,j,nin5)
 
1344
                r5=r*cr5 - s*ci5
 
1345
                s5=r*ci5 + s*cr5
 
1346
                r25 = r2 + r5
 
1347
                r34 = r3 + r4
 
1348
                s25 = s2 - s5
 
1349
                s34 = s3 - s4
 
1350
                zout(1,j,nout1) = r1 + r25 + r34
 
1351
                r = r1 + cos2*r25 + cos4*r34
 
1352
                s = sin2*s25 + sin4*s34
 
1353
                zout(1,j,nout2) = r - s
 
1354
                zout(1,j,nout5) = r + s
 
1355
                r = r1 + cos4*r25 + cos2*r34
 
1356
                s = sin4*s25 - sin2*s34
 
1357
                zout(1,j,nout3) = r - s
 
1358
                zout(1,j,nout4) = r + s
 
1359
                r25 = r2 - r5
 
1360
                r34 = r3 - r4
 
1361
                s25 = s2 + s5
 
1362
                s34 = s3 + s4
 
1363
                zout(2,j,nout1) = s1 + s25 + s34
 
1364
                r = s1 + cos2*s25 + cos4*s34
 
1365
                s = sin2*r25 + sin4*r34
 
1366
                zout(2,j,nout2) = r + s
 
1367
                zout(2,j,nout5) = r - s
 
1368
                r = s1 + cos4*s25 + cos2*s34
 
1369
                s = sin4*r25 - sin2*r34
 
1370
                zout(2,j,nout3) = r + s
 
1371
                zout(2,j,nout4) = r - s
 
1372
5100                continue
 
1373
        end if
 
1374
5000        continue
 
1375
       else if (now.eq.6) then
 
1376
!         .5d0*sqrt(3.d0)
 
1377
        bb=isign*0.8660254037844387d0
 
1378
 
 
1379
        ia=1
 
1380
        nin1=ia-after
 
1381
        nout1=ia-atn
 
1382
        do 6120,ib=1,before
 
1383
        nin1=nin1+after
 
1384
        nin2=nin1+atb
 
1385
        nin3=nin2+atb
 
1386
        nin4=nin3+atb
 
1387
        nin5=nin4+atb
 
1388
        nin6=nin5+atb
 
1389
        nout1=nout1+atn
 
1390
        nout2=nout1+after
 
1391
        nout3=nout2+after
 
1392
        nout4=nout3+after
 
1393
        nout5=nout4+after
 
1394
        nout6=nout5+after
 
1395
        do 6120,j=1,n1dfft
 
1396
        r2=zin(1,j,nin3)
 
1397
        s2=zin(2,j,nin3)
 
1398
        r3=zin(1,j,nin5)
 
1399
        s3=zin(2,j,nin5)
 
1400
        r=r2 + r3
 
1401
        s=s2 + s3
 
1402
        r1=zin(1,j,nin1)
 
1403
        s1=zin(2,j,nin1)
 
1404
        ur1 = r + r1
 
1405
        ui1 = s + s1
 
1406
        r1=r1 - .5d0*r
 
1407
        s1=s1 - .5d0*s
 
1408
        r=r2-r3
 
1409
        s=s2-s3
 
1410
        ur2 = r1 - s*bb
 
1411
        ui2 = s1 + r*bb
 
1412
        ur3 = r1 + s*bb
 
1413
        ui3 = s1 - r*bb
 
1414
 
 
1415
        r2=zin(1,j,nin6)
 
1416
        s2=zin(2,j,nin6)
 
1417
        r3=zin(1,j,nin2)
 
1418
        s3=zin(2,j,nin2)
 
1419
        r=r2 + r3
 
1420
        s=s2 + s3
 
1421
        r1=zin(1,j,nin4)
 
1422
        s1=zin(2,j,nin4)
 
1423
        vr1 = r + r1
 
1424
        vi1 = s + s1
 
1425
        r1=r1 - .5d0*r
 
1426
        s1=s1 - .5d0*s
 
1427
        r=r2-r3
 
1428
        s=s2-s3
 
1429
        vr2 = r1 - s*bb
 
1430
        vi2 = s1 + r*bb
 
1431
        vr3 = r1 + s*bb
 
1432
        vi3 = s1 - r*bb
 
1433
 
 
1434
        zout(1,j,nout1)=ur1+vr1
 
1435
        zout(2,j,nout1)=ui1+vi1
 
1436
        zout(1,j,nout5)=ur2+vr2
 
1437
        zout(2,j,nout5)=ui2+vi2
 
1438
        zout(1,j,nout3)=ur3+vr3
 
1439
        zout(2,j,nout3)=ui3+vi3
 
1440
        zout(1,j,nout4)=ur1-vr1
 
1441
        zout(2,j,nout4)=ui1-vi1
 
1442
        zout(1,j,nout2)=ur2-vr2
 
1443
        zout(2,j,nout2)=ui2-vi2
 
1444
        zout(1,j,nout6)=ur3-vr3
 
1445
        zout(2,j,nout6)=ui3-vi3
 
1446
6120        continue
 
1447
 
 
1448
        else
 
1449
        stop 'error fftstp'
 
1450
        end if
 
1451
 
 
1452
        return
 
1453
end subroutine fftstp