Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |