Home

Context Navigation

singleton_mod.f90 @ 4868

Last change on this file since 4868 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 39.1 KB

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

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

Context Navigation

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

Download in other formats: