Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 37.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |