Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |