Home

Context Navigation

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

Last change on this file since 4181 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 37.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |