Home

Context Navigation

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

Last change on this file since 4756 was 4591, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 39.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |