Home

Context Navigation

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

Last change on this file since 4372 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 37.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |