Home

Context Navigation

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

Last change on this file since 1850 was 1850, checked in by maronga, 9 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 37.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |