Home

Context Navigation

source: palm/trunk/SOURCE/singleton_mod.f90 @ 1884

Last change on this file since 1884 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 37.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |