Home

Context Navigation

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

Last change on this file since 2605 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 38.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |