Home

Context Navigation

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

Last change on this file since 4051 was 3761, checked in by raasch, 6 years ago
unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings
Property svn:keywords set to `Id`
File size: 38.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |