Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |