Home

Context Navigation

source: palm/trunk/SOURCE/singleton.f90 @ 1463

Last change on this file since 1463 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 34.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |