Home

Context Navigation

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

Last change on this file since 231 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 35.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |