Home

Context Navigation

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

Last change on this file since 2 was 1, checked in by raasch, 17 years ago
Initial repository layout and content
File size: 35.5 KB

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