Home

Context Navigation

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

Last change on this file since 1139 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 35.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |