Home

Context Navigation

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

Last change on this file since 2343 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 38.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |