Home

Context Navigation

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

Last change on this file since 1876 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 37.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |