Home

Context Navigation

source: palm/trunk/SOURCE/singleton.f90 @ 1320

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
Property svn:keywords set to `Id`
File size: 34.8 KB

Line
1	MODULE singleton
2
3	!-----------------------------------------------------------------------------
4	! Current revisions:
5	! -----------------
6	! kind-parameters added to all INTEGER and REAL declaration statements,
7	! kinds are defined in new module kinds,
8	! revision history before 2012 removed,
9	!
10	! Former revisions:
11	! -----------------
12	! Revision 1.1 2002/05/02 18:56:59 raasch
13	! Initial revision
14	!
15	!
16	! Description:
17	! ------------
18	! Multivariate Fast Fourier Transform
19	!
20	! Fortran 90 Implementation of Singleton's mixed-radix algorithm,
21	! RC Singleton, Stanford Research Institute, Sept. 1968.
22	!
23	! Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and
24	! John Beale.
25	!
26	! Fourier transforms can be computed either in place, using assumed size
27	! arguments, or by generic function, using assumed shape arguments.
28	!
29	!
30	! Public:
31	!
32	! fftkind kind parameter of complex arguments
33	! and function results.
34	!
35	! fft(array, dim, inv, stat) generic transform function
36	! COMPLEX(fftkind), DIMENSION(:,...,:), INTENT(IN) :: array
37	! INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim
38	! LOGICAL, INTENT(IN), OPTIONAL:: inv
39	! INTEGER, INTENT(OUT), OPTIONAL:: stat
40	!
41	! fftn(array, shape, dim, inv, stat) in place transform subroutine
42	! COMPLEX(fftkind), DIMENSION(*), INTENT(INOUT) :: array
43	! INTEGER, DIMENSION(:), INTENT(IN) :: shape
44	! INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim
45	! LOGICAL, INTENT(IN), OPTIONAL:: inv
46	! INTEGER, INTENT(OUT), OPTIONAL:: stat
47	!
48	!
49	! Formal Parameters:
50	!
51	! array The complex array to be transformed. array can be of arbitrary
52	! rank (i.e. up to seven).
53	!
54	! shape With subroutine fftn, the shape of the array to be transformed
55	! has to be passed separately, since fftradix - the internal trans-
56	! formation routine - will treat array always as one dimensional.
57	! The product of elements in shape must be the number of
58	! elements in array.
59	! Although passing array with assumed shape would have been nicer,
60	! I prefered assumed size in order to prevent the compiler from
61	! using a copy-in-copy-out mechanism. That would generally be
62	! necessary with fftn passing array to fftradix and with fftn
63	! being prepared for accepting non consecutive array sections.
64	! Using assumed size, it's up to the user to pass an array argu-
65	! ment, that can be addressed as continous one dimensional array
66	! without copying. Otherwise, transformation will not really be
67	! performed in place.
68	! On the other hand, since the rank of array and the size of
69	! shape needn't match, fftn is appropriate for handling more than
70	! seven dimensions.
71	! As far as function fft is concerned all this doesn't matter,
72	! because the argument will be copied anyway. Thus no extra
73	! shape argument is needed for fft.
74	!
75	! Optional Parameters:
76	!
77	! dim One dimensional integer array, containing the dimensions to be
78	! transformed. Default is (/1,...,N/) with N being the rank of
79	! array, i.e. complete transform. dim can restrict transformation
80	! to a subset of available dimensions. Its size must not exceed the
81	! rank of array or the size of shape respectivly.
82	!
83	! inv If .true., inverse transformation will be performed. Default is
84	! .false., i.e. forward transformation.
85	!
86	! stat If present, a system dependent nonzero status value will be
87	! returned in stat, if allocation of temporary storage failed.
88	!
89	!
90	! Scaling:
91	!
92	! Transformation results will always be scaled by the square root of the
93	! product of sizes of each dimension in dim. (See examples below)
94	!
95	!
96	! Examples:
97	!
98	! Let A be a LMN three dimensional complex array. Then
99	!
100	! result = fft(A)
101	!
102	! will produce a three dimensional transform, scaled by sqrt(LMN), while
103	!
104	! call fftn(A, SHAPE(A))
105	!
106	! will do the same in place.
107	!
108	! result = fft(A, dim=(/1,3/))
109	!
110	! will transform with respect to the first and the third dimension, scaled
111	! by sqrt(L*N).
112	!
113	! result = fft(fft(A), inv=.true.)
114	!
115	! should (approximately) reproduce A.
116	! With B having the same shape as A
117	!
118	! result = fft(fft(A) * CONJG(fft(B)), inv=.true.)
119	!
120	! will correlate A and B.
121	!
122	!
123	! Remarks:
124	!
125	! Following changes have been introduced with respect to fftn.c:
126	! - complex arguments and results are of type complex, rather than
127	! real an imaginary part separately.
128	! - increment parameter (magnitude of isign) has been dropped,
129	! inc is always one, direction of transform is given by inv.
130	! - maxf and maxp have been dropped. The amount of temporary storage
131	! needed is determined by the fftradix routine. Both fftn and fft
132	! can handle any size of array. (Maybe they take a lot of time and
133	! memory, but they will do it)
134	!
135	! Redesigning fftradix in a way, that it handles assumed shape arrays
136	! would have been desirable. However, I found it rather hard to do this
137	! in an efficient way. Problems were:
138	! - to prevent stride multiplications when indexing arrays. At least our
139	! compiler was not clever enough to discover that in fact additions
140	! would do the job as well. On the other hand, I haven't been clever
141	! enough to find an implementation using array operations.
142	! - fftradix is rather large and different versions would be necessaray
143	! for each possible rank of array.
144	! Consequently, in place transformation still needs the argument stored
145	! in a consecutive bunch of memory and can't be performed on array
146	! sections like A(100:199:-3, 50:1020). Calling fftn with such sections
147	! will most probably imply copy-in-copy-out. However, the function fft
148	! works with everything it gets and should be convenient to use.
149	!
150	! Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de>
151	! Restructured fftradix for better optimization. M. Steffens, 4 June 1997
152	!-----------------------------------------------------------------------------
153
154	USE kinds
155
156	IMPLICIT NONE
157
158	PRIVATE
159	PUBLIC:: fft, fftn
160
161	REAL(wp), PARAMETER:: sin60 = 0.86602540378443865_wp
162	REAL(wp), PARAMETER:: cos72 = 0.30901699437494742_wp
163	REAL(wp), PARAMETER:: sin72 = 0.95105651629515357_wp
164	REAL(wp), PARAMETER:: pi = 3.14159265358979323_wp
165
166	INTERFACE fft
167	MODULE PROCEDURE fft1d
168	MODULE PROCEDURE fft2d
169	MODULE PROCEDURE fft3d
170	MODULE PROCEDURE fft4d
171	MODULE PROCEDURE fft5d
172	MODULE PROCEDURE fft6d
173	MODULE PROCEDURE fft7d
174	END INTERFACE
175
176
177	CONTAINS
178
179
180	FUNCTION fft1d(array, dim, inv, stat) RESULT(ft)
181	!
182	!-- Formal parameters
183	COMPLEX(wp), DIMENSION(:), INTENT(IN) :: array
184	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
185	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
186	LOGICAL, INTENT(IN), OPTIONAL:: inv
187	!
188	!-- Function result
189	COMPLEX(wp), DIMENSION(SIZE(array, 1)):: ft
190
191	INTEGER(iwp):: ishape(1)
192
193	!
194	!-- Intrinsics used
195	INTRINSIC SIZE, SHAPE
196
197	ft = array
198	ishape = SHAPE( array )
199	CALL fftn(ft, ishape, inv = inv, stat = stat)
200
201	END FUNCTION fft1d
202
203
204	FUNCTION fft2d(array, dim, inv, stat) RESULT(ft)
205	!
206	!-- Formal parameters
207	COMPLEX(wp), DIMENSION(:,:), INTENT(IN) :: array
208	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
209	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
210	LOGICAL, INTENT(IN), OPTIONAL:: inv
211	!
212	!-- Function result
213	COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)):: ft
214
215	INTEGER(iwp) :: ishape(2)
216	!
217	!-- Intrinsics used
218	INTRINSIC SIZE, SHAPE
219
220	ft = array
221	ishape = SHAPE( array )
222	CALL fftn(ft, ishape, dim, inv, stat)
223
224	END FUNCTION fft2d
225
226
227	FUNCTION fft3d(array, dim, inv, stat) RESULT(ft)
228	!
229	!-- Formal parameters
230	COMPLEX(wp), DIMENSION(:,:,:), INTENT(IN) :: array
231	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
232	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
233	LOGICAL, INTENT(IN), OPTIONAL:: inv
234	!
235	!-- Function result
236	COMPLEX(wp), &
237	DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)):: ft
238
239	INTEGER(iwp) :: ishape(3)
240
241	!
242	!-- Intrinsics used
243	INTRINSIC SIZE, SHAPE
244
245	ft = array
246	ishape = SHAPE( array)
247	CALL fftn(ft, ishape, dim, inv, stat)
248
249	END FUNCTION fft3d
250
251
252	FUNCTION fft4d(array, dim, inv, stat) RESULT(ft)
253	!
254	!-- Formal parameters
255	COMPLEX(wp), DIMENSION(:,:,:,:), INTENT(IN) :: array
256	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
257	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
258	LOGICAL, INTENT(IN), OPTIONAL:: inv
259	!
260	!-- Function result
261	COMPLEX(wp), DIMENSION( &
262	SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)):: ft
263
264	INTEGER(iwp) :: ishape(4)
265	!
266	!-- Intrinsics used
267	INTRINSIC SIZE, SHAPE
268
269	ft = array
270	ishape = SHAPE( array )
271	CALL fftn(ft, ishape, dim, inv, stat)
272
273	END FUNCTION fft4d
274
275
276	FUNCTION fft5d(array, dim, inv, stat) RESULT(ft)
277	!
278	!-- Formal parameters
279	COMPLEX(wp), DIMENSION(:,:,:,:,:), INTENT(IN) :: array
280	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
281	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
282	LOGICAL, INTENT(IN), OPTIONAL:: inv
283	!
284	!-- Function result
285	COMPLEX(wp), DIMENSION( &
286	SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
287	SIZE(array, 5)):: ft
288
289	INTEGER(iwp) :: ishape(5)
290
291	!
292	!-- Intrinsics used
293	INTRINSIC SIZE, SHAPE
294
295	ft = array
296	ishape = SHAPE( array )
297	CALL fftn(ft, ishape, dim, inv, stat)
298
299	END FUNCTION fft5d
300
301
302	FUNCTION fft6d(array, dim, inv, stat) RESULT(ft)
303	!
304	!-- Formal parameters
305	COMPLEX(wp), DIMENSION(:,:,:,:,:,:), INTENT(IN) :: array
306	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
307	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
308	LOGICAL, INTENT(IN), OPTIONAL:: inv
309	!
310	!-- Function result
311	COMPLEX(wp), DIMENSION( &
312	SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
313	SIZE(array, 5), SIZE(array, 6)):: ft
314
315	INTEGER(iwp) :: ishape(6)
316
317	!
318	!-- Intrinsics used
319	INTRINSIC SIZE, SHAPE
320
321	ft = array
322	ishape = SHAPE( array )
323	CALL fftn(ft, ishape, dim, inv, stat)
324
325	END FUNCTION fft6d
326
327
328	FUNCTION fft7d(array, dim, inv, stat) RESULT(ft)
329	!
330	!-- Formal parameters
331	COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: array
332	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
333	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
334	LOGICAL, INTENT(IN), OPTIONAL:: inv
335	!
336	!-- Function result
337	COMPLEX(wp), DIMENSION( &
338	SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
339	SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)):: ft
340
341	INTEGER(iwp) :: ishape(7)
342
343	!
344	!-- Intrinsics used
345	INTRINSIC SIZE, SHAPE
346
347	ft = array
348	ishape = SHAPE( array )
349	CALL fftn(ft, ishape, dim, inv, stat)
350
351	END FUNCTION fft7d
352
353
354	SUBROUTINE fftn(array, shape, dim, inv, stat)
355	!
356	!-- Formal parameters
357	COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array
358	INTEGER(iwp), DIMENSION(:), INTENT(IN) :: shape
359	INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim
360	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
361	LOGICAL, INTENT(IN), OPTIONAL:: inv
362	!
363	!-- Local arrays
364	INTEGER(iwp), DIMENSION(SIZE(shape)):: d
365	!
366	!-- Local scalars
367	LOGICAL :: inverse
368	INTEGER(iwp) :: i, ndim, ntotal
369	REAL(wp):: scale
370	!
371	!-- Intrinsics used
372	INTRINSIC PRESENT, MIN, PRODUCT, SIZE, SQRT
373
374	!
375	!-- Optional parameter settings
376	IF (PRESENT(inv)) THEN
377	inverse = inv
378	ELSE
379	inverse = .FALSE.
380	END IF
381	IF (PRESENT(dim)) THEN
382	ndim = MIN(SIZE(dim), SIZE(d))
383	d(1:ndim) = DIM(1:ndim)
384	ELSE
385	ndim = SIZE(d)
386	d = (/(i, i = 1, SIZE(d))/)
387	END IF
388
389	ntotal = PRODUCT(shape)
390	scale = SQRT(1.0_wp / PRODUCT(shape(d(1:ndim))))
391	DO i = 1, ntotal
392	array(i) = CMPLX(REAL(array(i)) * scale, AIMAG(array(i)) * scale, &
393	KIND=wp)
394	END DO
395
396	DO i = 1, ndim
397	CALL fftradix(array, ntotal, shape(d(i)), PRODUCT(shape(1:d(i))), &
398	inverse, stat)
399	IF (PRESENT(stat)) THEN
400	IF (stat /=0) RETURN
401	END IF
402	END DO
403
404	END SUBROUTINE fftn
405
406
407	SUBROUTINE fftradix(array, ntotal, npass, nspan, inv, stat)
408	!
409	!-- Formal parameters
410	COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array
411	INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan
412	INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat
413	LOGICAL, INTENT(IN) :: inv
414	!
415	!-- Local arrays
416	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: ctmp
417	INTEGER(iwp), DIMENSION(BIT_SIZE(0)) :: factor
418	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: perm
419	REAL(wp), DIMENSION(:), ALLOCATABLE :: sine, cosine
420	!
421	!-- Local scalars
422	INTEGER(iwp) :: maxfactor, nfactor, nsquare, nperm
423	!
424	!-- Intrinsics used
425	INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, &
426	CMPLX, REAL, AIMAG
427
428	IF (npass <= 1) RETURN
429
430	CALL factorize(npass, factor, nfactor, nsquare)
431
432	maxfactor = MAXVAL(factor(:nfactor))
433	IF (nfactor - ISHFT(nsquare, 1) > 0) THEN
434	nperm = MAX(nfactor + 1, PRODUCT(factor(nsquare+1: nfactor-nsquare)) - 1)
435	ELSE
436	nperm = nfactor + 1
437	END IF
438
439	IF (PRESENT(stat)) THEN
440	ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT=stat)
441	IF (stat /= 0) RETURN
442	CALL transform(array, ntotal, npass, nspan, &
443	factor, nfactor, ctmp, sine, cosine, inv)
444	DEALLOCATE(sine, cosine, STAT=stat)
445	IF (stat /= 0) RETURN
446	ALLOCATE(perm(nperm), STAT=stat)
447	IF (stat /= 0) RETURN
448	CALL permute(array, ntotal, npass, nspan, &
449	factor, nfactor, nsquare, maxfactor, &
450	ctmp, perm)
451	DEALLOCATE(perm, ctmp, STAT=stat)
452	IF (stat /= 0) RETURN
453	ELSE
454	ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor))
455	CALL transform(array, ntotal, npass, nspan, &
456	factor, nfactor, ctmp, sine, cosine, inv)
457	DEALLOCATE(sine, cosine)
458	ALLOCATE(perm(nperm))
459	CALL permute(array, ntotal, npass, nspan, &
460	factor, nfactor, nsquare, maxfactor, &
461	ctmp, perm)
462	DEALLOCATE(perm, ctmp)
463	END IF
464
465
466	CONTAINS
467
468
469	SUBROUTINE factorize(npass, factor, nfactor, nsquare)
470	!
471	!-- Formal parameters
472	INTEGER(iwp), INTENT(IN) :: npass
473	INTEGER(iwp), DIMENSION(*), INTENT(OUT):: factor
474	INTEGER(iwp), INTENT(OUT):: nfactor, nsquare
475	!
476	!-- Local scalars
477	INTEGER(iwp):: j, jj, k
478
479	nfactor = 0
480	k = npass
481	DO WHILE (MOD(k, 16) == 0)
482	nfactor = nfactor + 1
483	factor(nfactor) = 4
484	k = k / 16
485	END DO
486	j = 3
487	jj = 9
488	DO
489	DO WHILE (MOD(k, jj) == 0)
490	nfactor = nfactor + 1
491	factor(nfactor) = j
492	k = k / jj
493	END DO
494	j = j + 2
495	jj = j * j
496	IF (jj > k) EXIT
497	END DO
498	IF (k <= 4) THEN
499	nsquare = nfactor
500	factor(nfactor + 1) = k
501	IF (k /= 1) nfactor = nfactor + 1
502	ELSE
503	IF (k - ISHFT(k / 4, 2) == 0) THEN
504	nfactor = nfactor + 1
505	factor(nfactor) = 2
506	k = k / 4
507	END IF
508	nsquare = nfactor
509	j = 2
510	DO
511	IF (MOD(k, j) == 0) THEN
512	nfactor = nfactor + 1
513	factor(nfactor) = j
514	k = k / j
515	END IF
516	j = ISHFT((j + 1) / 2, 1) + 1
517	IF (j > k) EXIT
518	END DO
519	END IF
520	IF (nsquare > 0) THEN
521	j = nsquare
522	DO
523	nfactor = nfactor + 1
524	factor(nfactor) = factor(j)
525	j = j - 1
526	IF (j==0) EXIT
527	END DO
528	END IF
529
530	END SUBROUTINE factorize
531
532
533	SUBROUTINE transform(array, ntotal, npass, nspan, &
534	factor, nfactor, ctmp, sine, cosine, inv) !-- compute fourier transform
535	!
536	!-- Formal parameters
537	COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array
538	COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp
539	INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan
540	INTEGER(iwp), DIMENSION(*), INTENT(IN) :: factor
541	INTEGER(iwp), INTENT(IN) :: nfactor
542	LOGICAL, INTENT(IN) :: inv
543	REAL(wp), DIMENSION(*), INTENT(OUT) :: sine, cosine
544	!
545	!-- Local scalars
546	INTEGER(iwp):: ii, ispan
547	INTEGER(iwp):: j, jc, jf, jj
548	INTEGER(iwp):: k, kk, kspan, k1, k2, k3, k4
549	INTEGER(iwp):: nn, nt
550	REAL(wp) :: s60, c72, s72, pi2, radf
551	REAL(wp) :: c1, s1, c2, s2, c3, s3, cd, sd, ak
552	COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm
553
554	c72 = cos72
555	IF (inv) THEN
556	s72 = sin72
557	s60 = sin60
558	pi2 = pi
559	ELSE
560	s72 = -sin72
561	s60 = -sin60
562	pi2 = -pi
563	END IF
564
565	nt = ntotal
566	nn = nt - 1
567	kspan = nspan
568	jc = nspan / npass
569	radf = pi2 * jc
570	pi2 = pi2 * 2.0_wp !-- use 2 PI from here on
571
572	ii = 0
573	jf = 0
574	DO
575	sd = radf / kspan
576	cd = SIN(sd)
577	cd = 2.0_wp * cd * cd
578	sd = SIN(sd + sd)
579	kk = 1
580	ii = ii + 1
581
582	SELECT CASE (factor(ii))
583	CASE (2)
584	!
585	!-- Transform for factor of 2 (including rotation factor)
586	kspan = kspan / 2
587	k1 = kspan + 2
588	DO
589	DO
590	k2 = kk + kspan
591	ck = array(k2)
592	array(k2) = array(kk)-ck
593	array(kk) = array(kk) + ck
594	kk = k2 + kspan
595	IF (kk > nn) EXIT
596	END DO
597	kk = kk - nn
598	IF (kk > jc) EXIT
599	END DO
600	IF (kk > kspan) RETURN
601	DO
602	c1 = 1.0_wp - cd
603	s1 = sd
604	DO
605	DO
606	DO
607	k2 = kk + kspan
608	ck = array(kk) - array(k2)
609	array(kk) = array(kk) + array(k2)
610	array(k2) = ck * CMPLX(c1, s1, KIND=wp)
611	kk = k2 + kspan
612	IF (kk >= nt) EXIT
613	END DO
614	k2 = kk - nt
615	c1 = -c1
616	kk = k1 - k2
617	IF (kk <= k2) EXIT
618	END DO
619	ak = c1 - (cd * c1 + sd * s1)
620	s1 = sd * c1 - cd * s1 + s1
621	c1 = 2.0_wp - (ak * ak + s1 * s1)
622	s1 = s1 * c1
623	c1 = c1 * ak
624	kk = kk + jc
625	IF (kk >= k2) EXIT
626	END DO
627	k1 = k1 + 1 + 1
628	kk = (k1 - kspan) / 2 + jc
629	IF (kk > jc + jc) EXIT
630	END DO
631
632	CASE (4) !-- transform for factor of 4
633	ispan = kspan
634	kspan = kspan / 4
635
636	DO
637	c1 = 1.0_wp
638	s1 = 0.0_wp
639	DO
640	DO
641	k1 = kk + kspan
642	k2 = k1 + kspan
643	k3 = k2 + kspan
644	ckp = array(kk) + array(k2)
645	ckm = array(kk) - array(k2)
646	cjp = array(k1) + array(k3)
647	cjm = array(k1) - array(k3)
648	array(kk) = ckp + cjp
649	cjp = ckp - cjp
650	IF (inv) THEN
651	ckp = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp)
652	ckm = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp)
653	ELSE
654	ckp = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp)
655	ckm = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp)
656	END IF
657	!
658	!-- Avoid useless multiplies
659	IF (s1 == 0.0_wp) THEN
660	array(k1) = ckp
661	array(k2) = cjp
662	array(k3) = ckm
663	ELSE
664	array(k1) = ckp * CMPLX(c1, s1, KIND=wp)
665	array(k2) = cjp * CMPLX(c2, s2, KIND=wp)
666	array(k3) = ckm * CMPLX(c3, s3, KIND=wp)
667	END IF
668	kk = k3 + kspan
669	IF (kk > nt) EXIT
670	END DO
671
672	c2 = c1 - (cd * c1 + sd * s1)
673	s1 = sd * c1 - cd * s1 + s1
674	c1 = 2.0_wp - (c2 * c2 + s1 * s1)
675	s1 = s1 * c1
676	c1 = c1 * c2
677	!
678	!-- Values of c2, c3, s2, s3 that will get used next time
679	c2 = c1 * c1 - s1 * s1
680	s2 = 2.0_wp * c1 * s1
681	c3 = c2 * c1 - s2 * s1
682	s3 = c2 * s1 + s2 * c1
683	kk = kk - nt + jc
684	IF (kk > kspan) EXIT
685	END DO
686	kk = kk - kspan + 1
687	IF (kk > jc) EXIT
688	END DO
689	IF (kspan == jc) RETURN
690
691	CASE default
692	!
693	!-- Transform for odd factors
694	k = factor(ii)
695	ispan = kspan
696	kspan = kspan / k
697
698	SELECT CASE (k)
699	CASE (3) !-- transform for factor of 3 (optional code)
700	DO
701	DO
702	k1 = kk + kspan
703	k2 = k1 + kspan
704	ck = array(kk)
705	cj = array(k1) + array(k2)
706	array(kk) = ck + cj
707	ck = ck - CMPLX( &
708	0.5_wp * REAL (cj), &
709	0.5_wp * AIMAG(cj), &
710	KIND=wp)
711	cj = CMPLX( &
712	(REAL (array(k1)) - REAL (array(k2))) * s60, &
713	(AIMAG(array(k1)) - AIMAG(array(k2))) * s60, &
714	KIND=wp)
715	array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
716	array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
717	kk = k2 + kspan
718	IF (kk >= nn) EXIT
719	END DO
720	kk = kk - nn
721	IF (kk > kspan) EXIT
722	END DO
723
724	CASE (5) !-- transform for factor of 5 (optional code)
725	c2 = c72 * c72 - s72 * s72
726	s2 = 2.0_wp * c72 * s72
727	DO
728	DO
729	k1 = kk + kspan
730	k2 = k1 + kspan
731	k3 = k2 + kspan
732	k4 = k3 + kspan
733	ckp = array(k1) + array(k4)
734	ckm = array(k1) - array(k4)
735	cjp = array(k2) + array(k3)
736	cjm = array(k2) - array(k3)
737	cc = array(kk)
738	array(kk) = cc + ckp + cjp
739	ck = CMPLX(REAL(ckp) * c72, AIMAG(ckp) * c72, &
740	KIND=wp) + &
741	CMPLX(REAL(cjp) * c2, AIMAG(cjp) * c2, &
742	KIND=wp) + cc
743	cj = CMPLX(REAL(ckm) * s72, AIMAG(ckm) * s72, &
744	KIND=wp) + &
745	CMPLX(REAL(cjm) * s2, AIMAG(cjm) * s2, &
746	KIND=wp)
747	array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
748	array(k4) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
749	ck = CMPLX(REAL(ckp) * c2, AIMAG(ckp) * c2, &
750	KIND=wp) + &
751	CMPLX(REAL(cjp) * c72, AIMAG(cjp) * c72, &
752	KIND=wp) + cc
753	cj = CMPLX(REAL(ckm) * s2, AIMAG(ckm) * s2, &
754	KIND=wp) - &
755	CMPLX(REAL(cjm) * s72, AIMAG(cjm) * s72, &
756	KIND=wp)
757	array(k2) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
758	array(k3) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
759	kk = k4 + kspan
760	IF (kk >= nn) EXIT
761	END DO
762	kk = kk - nn
763	IF (kk > kspan) EXIT
764	END DO
765
766	CASE default
767	IF (k /= jf) THEN
768	jf = k
769	s1 = pi2 / k
770	c1 = COS(s1)
771	s1 = SIN(s1)
772	cosine (jf) = 1.0_wp
773	sine (jf) = 0.0_wp
774	j = 1
775	DO
776	cosine (j) = cosine (k) * c1 + sine (k) * s1
777	sine (j) = cosine (k) * s1 - sine (k) * c1
778	k = k-1
779	cosine (k) = cosine (j)
780	sine (k) = -sine (j)
781	j = j + 1
782	IF (j >= k) EXIT
783	END DO
784	END IF
785	DO
786	DO
787	k1 = kk
788	k2 = kk + ispan
789	cc = array(kk)
790	ck = cc
791	j = 1
792	k1 = k1 + kspan
793	DO
794	k2 = k2 - kspan
795	j = j + 1
796	ctmp(j) = array(k1) + array(k2)
797	ck = ck + ctmp(j)
798	j = j + 1
799	ctmp(j) = array(k1) - array(k2)
800	k1 = k1 + kspan
801	IF (k1 >= k2) EXIT
802	END DO
803	array(kk) = ck
804	k1 = kk
805	k2 = kk + ispan
806	j = 1
807	DO
808	k1 = k1 + kspan
809	k2 = k2 - kspan
810	jj = j
811	ck = cc
812	cj = (0.0_wp, 0.0_wp)
813	k = 1
814	DO
815	k = k + 1
816	ck = ck + CMPLX( &
817	REAL (ctmp(k)) * cosine(jj), &
818	AIMAG(ctmp(k)) * cosine(jj), KIND=wp)
819	k = k + 1
820	cj = cj + CMPLX( &
821	REAL (ctmp(k)) * sine(jj), &
822	AIMAG(ctmp(k)) * sine(jj), KIND=wp)
823	jj = jj + j
824	IF (jj > jf) jj = jj - jf
825	IF (k >= jf) EXIT
826	END DO
827	k = jf - j
828	array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), &
829	KIND=wp)
830	array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), &
831	KIND=wp)
832	j = j + 1
833	IF (j >= k) EXIT
834	END DO
835	kk = kk + ispan
836	IF (kk > nn) EXIT
837	END DO
838	kk = kk - nn
839	IF (kk > kspan) EXIT
840	END DO
841
842	END SELECT
843	!
844	!-- Multiply by rotation factor (except for factors of 2 and 4)
845	IF (ii == nfactor) RETURN
846	kk = jc + 1
847	DO
848	c2 = 1.0_wp - cd
849	s1 = sd
850	DO
851	c1 = c2
852	s2 = s1
853	kk = kk + kspan
854	DO
855	DO
856	array(kk) = CMPLX(c2, s2, KIND=wp) * array(kk)
857	kk = kk + ispan
858	IF (kk > nt) EXIT
859	END DO
860	ak = s1 * s2
861	s2 = s1 * c2 + c1 * s2
862	c2 = c1 * c2 - ak
863	kk = kk - nt + kspan
864	IF (kk > ispan) EXIT
865	END DO
866	c2 = c1 - (cd * c1 + sd * s1)
867	s1 = s1 + sd * c1 - cd * s1
868	c1 = 2.0_wp - (c2 * c2 + s1 * s1)
869	s1 = s1 * c1
870	c2 = c2 * c1
871	kk = kk - ispan + jc
872	IF (kk > kspan) EXIT
873	END DO
874	kk = kk - kspan + jc + 1
875	IF (kk > jc + jc) EXIT
876	END DO
877
878	END SELECT
879	END DO
880	END SUBROUTINE transform
881
882
883	SUBROUTINE permute(array, ntotal, npass, nspan, &
884	factor, nfactor, nsquare, maxfactor, &
885	ctmp, perm)
886	!
887	!-- Formal parameters
888	COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array
889	COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp
890	INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan
891	INTEGER(iwp), DIMENSION(*), INTENT(IN OUT):: factor
892	INTEGER(iwp), INTENT(IN) :: nfactor, nsquare
893	INTEGER(iwp), INTENT(IN) :: maxfactor
894	INTEGER(iwp), DIMENSION(*), INTENT(OUT) :: perm
895	!
896	!-- Local scalars
897	COMPLEX(wp) :: ck
898	INTEGER(iwp):: ii, ispan
899	INTEGER(iwp):: j, jc, jj
900	INTEGER(iwp):: k, kk, kspan, kt, k1, k2, k3
901	INTEGER(iwp):: nn, nt
902	!
903	!-- Permute the results to normal order---done in two stages
904	!-- Permutation for square factors of n
905
906	nt = ntotal
907	nn = nt - 1
908	kt = nsquare
909	kspan = nspan
910	jc = nspan / npass
911
912	perm (1) = nspan
913	IF (kt > 0) THEN
914	k = kt + kt + 1
915	IF (nfactor < k) k = k - 1
916	j = 1
917	perm (k + 1) = jc
918	DO
919	perm (j + 1) = perm (j) / factor(j)
920	perm (k) = perm (k + 1) * factor(j)
921	j = j + 1
922	k = k - 1
923	IF (j >= k) EXIT
924	END DO
925	k3 = perm (k + 1)
926	kspan = perm (2)
927	kk = jc + 1
928	k2 = kspan + 1
929	j = 1
930
931	IF (npass /= ntotal) THEN
932	permute_multi: DO
933	DO
934	DO
935	k = kk + jc
936	DO
937	!
938	!-- Swap array(kk) <> array(k2)
939	ck = array(kk)
940	array(kk) = array(k2)
941	array(k2) = ck
942	kk = kk + 1
943	k2 = k2 + 1
944	IF (kk >= k) EXIT
945	END DO
946	kk = kk + nspan - jc
947	k2 = k2 + nspan - jc
948	IF (kk >= nt) EXIT
949	END DO
950	kk = kk - nt + jc
951	k2 = k2 - nt + kspan
952	IF (k2 >= nspan) EXIT
953	END DO
954	DO
955	DO
956	k2 = k2 - perm (j)
957	j = j + 1
958	k2 = perm (j + 1) + k2
959	IF (k2 <= perm (j)) EXIT
960	END DO
961	j = 1
962	DO
963	IF (kk < k2) CYCLE permute_multi
964	kk = kk + jc
965	k2 = k2 + kspan
966	IF (k2 >= nspan) EXIT
967	END DO
968	IF (kk >= nspan) EXIT
969	END DO
970	EXIT
971	END DO permute_multi
972	ELSE
973	permute_single: DO
974	DO
975	!
976	!-- Swap array(kk) <> array(k2)
977	ck = array(kk)
978	array(kk) = array(k2)
979	array(k2) = ck
980	kk = kk + 1
981	k2 = k2 + kspan
982	IF (k2 >= nspan) EXIT
983	END DO
984	DO
985	DO
986	k2 = k2 - perm (j)
987	j = j + 1
988	k2 = perm (j + 1) + k2
989	IF (k2 <= perm (j)) EXIT
990	END DO
991	j = 1
992	DO
993	IF (kk < k2) CYCLE permute_single
994	kk = kk + 1
995	k2 = k2 + kspan
996	IF (k2 >= nspan) EXIT
997	END DO
998	IF (kk >= nspan) EXIT
999	END DO
1000	EXIT
1001	END DO permute_single
1002	END IF
1003	jc = k3
1004	END IF
1005
1006	IF (ISHFT(kt, 1) + 1 >= nfactor) RETURN
1007
1008	ispan = perm (kt + 1)
1009	!
1010	!-- Permutation for square-free factors of n
1011	j = nfactor - kt
1012	factor(j + 1) = 1
1013	DO
1014	factor(j) = factor(j) * factor(j+1)
1015	j = j - 1
1016	IF (j == kt) EXIT
1017	END DO
1018	kt = kt + 1
1019	nn = factor(kt) - 1
1020	j = 0
1021	jj = 0
1022	DO
1023	k = kt + 1
1024	k2 = factor(kt)
1025	kk = factor(k)
1026	j = j + 1
1027	IF (j > nn) EXIT !-- exit infinite loop
1028	jj = jj + kk
1029	DO WHILE (jj >= k2)
1030	jj = jj - k2
1031	k2 = kk
1032	k = k + 1
1033	kk = factor(k)
1034	jj = jj + kk
1035	END DO
1036	perm (j) = jj
1037	END DO
1038	!
1039	!-- Determine the permutation cycles of length greater than 1
1040	j = 0
1041	DO
1042	DO
1043	j = j + 1
1044	kk = perm(j)
1045	IF (kk >= 0) EXIT
1046	END DO
1047	IF (kk /= j) THEN
1048	DO
1049	k = kk
1050	kk = perm (k)
1051	perm (k) = -kk
1052	IF (kk == j) EXIT
1053	END DO
1054	k3 = kk
1055	ELSE
1056	perm (j) = -j
1057	IF (j == nn) EXIT !-- exit infinite loop
1058	END IF
1059	END DO
1060	!
1061	!-- Reorder a and b, following the permutation cycles
1062	DO
1063	j = k3 + 1
1064	nt = nt - ispan
1065	ii = nt - 1 + 1
1066	IF (nt < 0) EXIT !-- exit infinite loop
1067	DO
1068	DO
1069	j = j-1
1070	IF (perm(j) >= 0) EXIT
1071	END DO
1072	jj = jc
1073	DO
1074	kspan = jj
1075	IF (jj > maxfactor) kspan = maxfactor
1076	jj = jj - kspan
1077	k = perm(j)
1078	kk = jc * k + ii + jj
1079	k1 = kk + kspan
1080	k2 = 0
1081	DO
1082	k2 = k2 + 1
1083	ctmp(k2) = array(k1)
1084	k1 = k1 - 1
1085	IF (k1 == kk) EXIT
1086	END DO
1087	DO
1088	k1 = kk + kspan
1089	k2 = k1 - jc * (k + perm(k))
1090	k = -perm(k)
1091	DO
1092	array(k1) = array(k2)
1093	k1 = k1 - 1
1094	k2 = k2 - 1
1095	IF (k1 == kk) EXIT
1096	END DO
1097	kk = k2
1098	IF (k == j) EXIT
1099	END DO
1100	k1 = kk + kspan
1101	k2 = 0
1102	DO
1103	k2 = k2 + 1
1104	array(k1) = ctmp(k2)
1105	k1 = k1 - 1
1106	IF (k1 == kk) EXIT
1107	END DO
1108	IF (jj == 0) EXIT
1109	END DO
1110	IF (j == 1) EXIT
1111	END DO
1112	END DO
1113
1114	END SUBROUTINE permute
1115
1116	END SUBROUTINE fftradix
1117
1118	END MODULE singleton

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |