Home

Context Navigation

source: palm/trunk/SOURCE/temperton_fft_mod.f90 @ 3945

Last change on this file since 3945 was 3725, checked in by raasch, 6 years ago
modifications to avoid compiler warnings about unused variables, temperton-fft: GOTO statements replaced, file re-formatted corresponding to coding standards, ssh-calls for compilations on remote systems modified to avoid output of login messages on specific systems changed again (palmbuild, reverted as before r3549), error messages for failed restarts extended (palmrun)
Property svn:keywords set to `Id`
File size: 71.2 KB

Rev	Line
[1850]	1	!> @file temperton_fft_mod.f90
[1]	2	!------------------------------------------------------------------------------!
[1682]	3	!
[258]	4	! Current revisions:
[1]	5	! -----------------
[1343]	6	!
[3049]	7	!
[1321]	8	! Former revisions:
	9	! -----------------
	10	! $Id: temperton_fft_mod.f90 3725 2019-02-07 10:11:02Z raasch $
[3725]	11	! GOTO statements replaced, file re-formatted corresponding to coding standards
	12	!
	13	! 3241 2018-09-12 15:02:00Z raasch
[3241]	14	! unused variables removed
	15	!
	16	! 3049 2018-05-29 13:52:36Z Giersch
[3049]	17	! Error messages revised
	18	!
	19	! 3045 2018-05-28 07:55:41Z Giersch
[3045]	20	! Error message revised
	21	!
	22	! 2300 2017-06-29 13:31:14Z raasch
[2300]	23	! NEC related CPP directives removed
	24	!
	25	! 1851 2016-04-08 13:32:50Z maronga
[1321]	26	!
[1851]	27	! 1850 2016-04-08 13:29:27Z maronga
	28	! Module renamed
	29	!
[1683]	30	! 1682 2015-10-07 23:56:08Z knoop
	31	! Code annotations made doxygen readable
	32	!
[1343]	33	! 1342 2014-03-26 17:04:47Z kanani
	34	! REAL constants defined as wp-kind
	35	!
[1323]	36	! 1322 2014-03-20 16:38:49Z raasch
	37	! REAL constants defined as wp-kind
	38	!
[1321]	39	! 1320 2014-03-20 08:40:49Z raasch
[1320]	40	! ONLY-attribute added to USE-statements,
	41	! kind-parameters added to all INTEGER and REAL declaration statements,
	42	! kinds are defined in new module kinds,
	43	! old module precision_kind is removed,
	44	! revision history before 2012 removed,
	45	! comment fields (!:) to be used for variable explanations added to
	46	! all variable declaration statements
[1]	47	!
	48	! Revision 1.1 2003/03/12 16:41:59 raasch
	49	! Initial revision
	50	!
	51	!
	52	! Description:
	53	! ------------
[1682]	54	!> Fast Fourier transformation developed by Clive Temperton, ECMWF.
[1]	55	!------------------------------------------------------------------------------!
[1682]	56	MODULE temperton_fft
[1]	57
[1320]	58	USE kinds
	59
[1]	60	IMPLICIT NONE
	61
	62	PRIVATE
	63
	64	PUBLIC set99, fft991cy
	65
[2300]	66	INTEGER(iwp), PARAMETER :: nfft = 32 !< maximum length of calls to *fft
	67	INTEGER(iwp), PARAMETER :: nout = 6 !< standard output stream
	68
[1]	69	CONTAINS
	70
[1682]	71	!------------------------------------------------------------------------------!
	72	! Description:
	73	! ------------
[3725]	74	!> Calls Fortran-versions of fft's.
[1682]	75	!>
	76	!> Method:
	77	!>
	78	!> Subroutine 'fft991cy' - multiple fast real periodic transform
	79	!> supersedes previous routine 'fft991cy'.
	80	!>
	81	!> Real transform of length n performed by removing redundant
	82	!> operations from complex transform of length n.
	83	!>
	84	!> a is the array containing input & output data.
	85	!> work is an area of size (n+1)*min(lot,nfft).
	86	!> trigs is a previously prepared list of trig function values.
	87	!> ifax is a previously prepared list of factors of n.
	88	!> inc is the increment within each data 'vector'
	89	!> (e.g. inc=1 for consecutively stored data).
	90	!> jump is the increment between the start of each data vector.
	91	!> n is the length of the data vectors.
	92	!> lot is the number of data vectors.
	93	!> isign = +1 for transform from spectral to gridpoint
	94	!> = -1 for transform from gridpoint to spectral.
	95	!>
	96	!> ordering of coefficients:
	97	!> a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
	98	!> where b(0)=b(n/2)=0; (n+2) locations required.
	99	!>
	100	!> ordering of data:
	101	!> x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
	102	!>
	103	!> Vectorization is achieved on cray by doing the transforms
	104	!> in parallel.
	105	!>
	106	!> n must be composed of factors 2,3 & 5 but does not have to be even.
	107	!>
	108	!> definition of transforms:
	109	!>
	110	!> isign=+1: x(j)=sum(k=0,.,n-1)(c(k)exp(2ijk*pi/n))
	111	!> where c(k)=a(k)+ib(k) and c(n-k)=a(k)-ib(k)
	112	!>
	113	!> isign=-1: a(k)=(1/n)sum(j=0,.,n-1)(x(j)cos(2jk*pi/n))
	114	!> b(k)=-(1/n)sum(j=0,.,n-1)(x(j)sin(2jk*pi/n))
	115	!>
	116	!> calls fortran-versions of fft's !!!
	117	!> dimension a(n),work(n),trigs(n),ifax(1)
	118	!------------------------------------------------------------------------------!
[3725]	119	SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign )
[1]	120
[1320]	121	USE kinds
[1]	122
	123	IMPLICIT NONE
	124
[1682]	125	INTEGER(iwp) :: inc !<
	126	INTEGER(iwp) :: isign !<
	127	INTEGER(iwp) :: jump !<
	128	INTEGER(iwp) :: lot !<
	129	INTEGER(iwp) :: n !<
[1]	130
[1682]	131	REAL(wp) :: a(*) !<
	132	REAL(wp) :: trigs(*) !<
	133	REAL(wp) :: work(*) !<
	134	INTEGER(iwp) :: ifax(*) !<
[1]	135
[1682]	136	INTEGER(iwp) :: i !<
	137	INTEGER(iwp) :: ia !<
	138	INTEGER(iwp) :: ibase !<
	139	INTEGER(iwp) :: ierr !<
	140	INTEGER(iwp) :: ifac !<
	141	INTEGER(iwp) :: igo !<
	142	INTEGER(iwp) :: ii !<
	143	INTEGER(iwp) :: istart !<
	144	INTEGER(iwp) :: ix !<
	145	INTEGER(iwp) :: iz !<
	146	INTEGER(iwp) :: j !<
	147	INTEGER(iwp) :: jbase !<
	148	INTEGER(iwp) :: jj !<
	149	INTEGER(iwp) :: k !<
	150	INTEGER(iwp) :: la !<
	151	INTEGER(iwp) :: nb !<
	152	INTEGER(iwp) :: nblox !<
	153	INTEGER(iwp) :: nfax !<
	154	INTEGER(iwp) :: nvex !<
	155	INTEGER(iwp) :: nx !<
[1]	156
	157
[3725]	158	IF ( ifax(10) /= n ) CALL set99( trigs, ifax, n )
	159	nfax = ifax(1)
	160	nx = n + 1
	161	IF ( MOD(n,2) == 1 ) nx = n
	162	nblox = 1 + (lot-1) / nfft
	163	nvex = lot - (nblox-1) * nfft
[1]	164
	165
[3725]	166	IF ( isign == 1 ) THEN
	167	!
	168	!-- Backward fft: spectral to gridpoint transform
[1]	169
[3725]	170	istart = 1
	171	DO nb = 1, nblox
[1]	172
	173	ia = istart
[3725]	174	i = istart
	175	!DIR$ IVDEP
	176	!CDIR NODEP
	177	!OCL NOVREC
	178	DO j = 1, nvex
	179	a(i+inc) = 0.5_wp * a(i)
	180	i = i + jump
	181	ENDDO
	182	IF ( MOD(n,2) /= 1 ) THEN
	183	i = istart + n * inc
	184	DO j = 1, nvex
	185	a(i) = 0.5_wp * a(i)
	186	i = i + jump
	187	ENDDO
	188	ENDIF
	189	ia = istart + inc
	190	la = 1
	191	igo = + 1
[1]	192
[3725]	193	DO k = 1, nfax
[1]	194
[3725]	195	ifac = ifax(k+1)
	196	ierr = -1
	197	IF ( igo /= -1 ) THEN
	198	CALL rpassm( a(ia), a(ia+lainc), work(1), work(ifacla+1), trigs, inc, 1, jump, &
	199	nx, nvex, n, ifac, la, ierr )
	200	ELSE
	201	CALL rpassm( work(1), work(la+1), a(ia), a(ia+ifaclainc), trigs, 1, inc, nx, &
	202	jump, nvex, n, ifac, la, ierr )
	203	ENDIF
	204	!
	205	!-- Following messages shouldn't appear in PALM applications
	206	IF ( ierr /= 0 ) THEN
	207	SELECT CASE (ierr)
	208	CASE (:-1)
	209	WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
	210	CASE (0)
	211	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
	212	CASE (1:)
	213	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
	214	END SELECT
	215	RETURN
	216	ENDIF
[1]	217
[3725]	218	la = ifac * la
	219	igo = -igo
	220	ia = istart
[1]	221
[3725]	222	ENDDO
	223	!
	224	!-- If necessary, copy results back to a
	225	IF ( MOD(nfax,2) /= 0 ) THEN
	226	ibase = 1
	227	jbase = ia
	228	DO jj = 1, nvex
	229	i = ibase
	230	j = jbase
	231	DO ii = 1, n
	232	a(j) = work(i)
	233	i = i + 1
	234	j = j + inc
	235	ENDDO
	236	ibase = ibase + nx
	237	jbase = jbase + jump
	238	ENDDO
	239	ENDIF
	240	!
	241	!-- Fill in zeros at end
	242	ix = istart + n*inc
	243	!DIR$ IVDEP
	244	!CDIR NODEP
	245	!OCL NOVREC
	246	DO j = 1, nvex
	247	a(ix) = 0.0_wp
	248	a(ix+inc) = 0.0_wp
	249	ix = ix + jump
	250	ENDDO
	251	istart = istart + nvex*jump
	252	nvex = nfft
[1]	253
[3725]	254	ENDDO
[1]	255
[3725]	256	ELSEIF ( isign == -1 ) THEN
	257	!
	258	!-- Forward fft: gridpoint to spectral transform
	259	istart = 1
	260	DO nb = 1, nblox
	261	ia = istart
	262	la = n
	263	igo = + 1
[1]	264
[3725]	265	DO k = 1, nfax
[1]	266
[3725]	267	ifac = ifax(nfax+2-k)
	268	la = la / ifac
	269	ierr = -1
	270	IF ( igo /= -1 ) THEN
	271	CALL qpassm( a(ia), a(ia+ifaclainc), work(1), work(la+1), trigs, inc, 1, jump, &
	272	nx, nvex, n, ifac, la, ierr )
	273	ELSE
	274	CALL qpassm( work(1), work(ifacla+1), a(ia), a(ia+lainc), trigs, 1, inc, nx, &
	275	jump, nvex, n, ifac, la, ierr )
	276	ENDIF
	277	!
	278	!-- Following messages shouldn't appear in PALM applications
	279	IF ( ierr /= 0 ) THEN
	280	SELECT CASE (ierr)
	281	CASE (:-1)
	282	WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
	283	CASE (0)
	284	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
	285	CASE (1:)
	286	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
	287	END SELECT
	288	RETURN
	289	ENDIF
[1]	290
[3725]	291	igo = -igo
	292	ia = istart + inc
[1]	293
[3725]	294	ENDDO
	295	!
	296	!-- If necessary, copy results back to a
	297	IF ( MOD(nfax,2) /= 0 ) THEN
	298	ibase = 1
	299	jbase = ia
	300	DO jj = 1, nvex
	301	i = ibase
	302	j = jbase
	303	DO ii = 1, n
	304	a(j) = work(i)
	305	i = i + 1
	306	j = j + inc
	307	ENDDO
	308	ibase = ibase + nx
	309	jbase = jbase + jump
	310	ENDDO
	311	ENDIF
	312	!
	313	!-- Shift a(0) and fill in zero imag parts
	314	ix = istart
	315	!DIR$ IVDEP
	316	!CDIR NODEP
	317	!OCL NOVREC
	318	DO j = 1, nvex
	319	a(ix) = a(ix+inc)
	320	a(ix+inc) = 0.0_wp
	321	ix = ix + jump
	322	ENDDO
	323	IF ( MOD(n,2) /= 1 ) THEN
	324	iz = istart + (n+1) * inc
	325	DO j = 1, nvex
	326	a(iz) = 0.0_wp
	327	iz = iz + jump
	328	ENDDO
	329	ENDIF
	330	istart = istart + nvex * jump
	331	nvex = nfft
[1]	332
[3725]	333	ENDDO
[1]	334
[3725]	335	ENDIF
[1]	336
[3725]	337	END SUBROUTINE fft991cy
[1]	338
[1682]	339	!------------------------------------------------------------------------------!
	340	! Description:
	341	! ------------
	342	!> Performs one pass through data as part of
	343	!> multiple real fft (fourier analysis) routine.
	344	!>
	345	!> Method:
	346	!>
	347	!> a is first real input vector
	348	!> equivalence b(1) with a(ifaclainc1+1)
	349	!> c is first real output vector
	350	!> equivalence d(1) with c(la*inc2+1)
	351	!> trigs is a precalculated list of sines & cosines
	352	!> inc1 is the addressing increment for a
	353	!> inc2 is the addressing increment for c
	354	!> inc3 is the increment between input vectors a
	355	!> inc4 is the increment between output vectors c
	356	!> lot is the number of vectors
	357	!> n is the length of the vectors
	358	!> ifac is the current factor of n
	359	!> la = n/(product of factors used so far)
	360	!> ierr is an error indicator:
	361	!> 0 - pass completed without error
	362	!> 1 - lot greater than nfft
	363	!> 2 - ifac not catered for
	364	!> 3 - ifac only catered for if la=n/ifac
	365	!------------------------------------------------------------------------------!
[3725]	366	SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
[1]	367
[1320]	368	USE kinds
[1]	369
[1320]	370	IMPLICIT NONE
	371
[1682]	372	INTEGER(iwp) :: ierr !<
	373	INTEGER(iwp) :: ifac !<
	374	INTEGER(iwp) :: inc1 !<
	375	INTEGER(iwp) :: inc2 !<
	376	INTEGER(iwp) :: inc3 !<
	377	INTEGER(iwp) :: inc4 !<
	378	INTEGER(iwp) :: la !<
	379	INTEGER(iwp) :: lot !<
	380	INTEGER(iwp) :: n !<
[1]	381
[3725]	382	!
	383	!-- Arrays are dimensioned with n
[1682]	384	REAL(wp) :: a(*) !<
	385	REAL(wp) :: b(*) !<
	386	REAL(wp) :: c(*) !<
	387	REAL(wp) :: d(*) !<
	388	REAL(wp) :: trigs(*) !<
[1320]	389
[1682]	390	REAL(wp) :: a0 !<
	391	REAL(wp) :: a1 !<
	392	REAL(wp) :: a10 !<
	393	REAL(wp) :: a11 !<
	394	REAL(wp) :: a2 !<
	395	REAL(wp) :: a20 !<
	396	REAL(wp) :: a21 !<
	397	REAL(wp) :: a3 !<
	398	REAL(wp) :: a4 !<
	399	REAL(wp) :: a5 !<
	400	REAL(wp) :: a6 !<
	401	REAL(wp) :: b0 !<
	402	REAL(wp) :: b1 !<
	403	REAL(wp) :: b10 !<
	404	REAL(wp) :: b11 !<
	405	REAL(wp) :: b2 !<
	406	REAL(wp) :: b20 !<
	407	REAL(wp) :: b21 !<
	408	REAL(wp) :: b3 !<
	409	REAL(wp) :: b4 !<
	410	REAL(wp) :: b5 !<
	411	REAL(wp) :: b6 !<
	412	REAL(wp) :: c1 !<
	413	REAL(wp) :: c2 !<
	414	REAL(wp) :: c3 !<
	415	REAL(wp) :: c4 !<
	416	REAL(wp) :: c5 !<
	417	REAL(wp) :: qrt5 !<
	418	REAL(wp) :: s1 !<
	419	REAL(wp) :: s2 !<
	420	REAL(wp) :: s3 !<
	421	REAL(wp) :: s4 !<
	422	REAL(wp) :: s5 !<
	423	REAL(wp) :: sin36 !<
	424	REAL(wp) :: sin45 !<
	425	REAL(wp) :: sin60 !<
	426	REAL(wp) :: sin72 !<
	427	REAL(wp) :: z !<
	428	REAL(wp) :: zqrt5 !<
	429	REAL(wp) :: zsin36 !<
	430	REAL(wp) :: zsin45 !<
	431	REAL(wp) :: zsin60 !<
	432	REAL(wp) :: zsin72 !<
[1]	433
[1682]	434	INTEGER(iwp) :: i !<
	435	INTEGER(iwp) :: ia !<
	436	INTEGER(iwp) :: ib !<
	437	INTEGER(iwp) :: ibase !<
	438	INTEGER(iwp) :: ic !<
	439	INTEGER(iwp) :: id !<
	440	INTEGER(iwp) :: ie !<
	441	INTEGER(iwp) :: if !<
	442	INTEGER(iwp) :: ig !<
	443	INTEGER(iwp) :: igo !<
	444	INTEGER(iwp) :: ih !<
	445	INTEGER(iwp) :: iink !<
	446	INTEGER(iwp) :: ijk !<
	447	INTEGER(iwp) :: ijump !<
	448	INTEGER(iwp) :: j !<
	449	INTEGER(iwp) :: ja !<
	450	INTEGER(iwp) :: jb !<
	451	INTEGER(iwp) :: jbase !<
	452	INTEGER(iwp) :: jc !<
	453	INTEGER(iwp) :: jd !<
	454	INTEGER(iwp) :: je !<
	455	INTEGER(iwp) :: jf !<
	456	INTEGER(iwp) :: jink !<
	457	INTEGER(iwp) :: k !<
	458	INTEGER(iwp) :: kb !<
	459	INTEGER(iwp) :: kc !<
	460	INTEGER(iwp) :: kd !<
	461	INTEGER(iwp) :: ke !<
	462	INTEGER(iwp) :: kf !<
	463	INTEGER(iwp) :: kstop !<
	464	INTEGER(iwp) :: l !<
	465	INTEGER(iwp) :: m !<
[1320]	466
[1]	467
[3725]	468	DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
	469	qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
[1]	470
	471
[3725]	472	ierr = 0
[1]	473
[3725]	474	IF ( lot > nfft ) THEN
	475	ierr = 1
	476	RETURN
	477	ENDIF
[1]	478
[3725]	479	m = n / ifac
	480	iink = la * inc1
	481	jink = la * inc2
	482	ijump = (ifac-1) * iink
	483	kstop = ( n-ifac ) / ( 2*ifac )
	484
[1]	485	ibase = 0
	486	jbase = 0
	487	igo = ifac - 1
[3725]	488	IF ( igo == 7 ) igo = 6
[1]	489
[3725]	490	IF (igo < 1 .OR. igo > 6 ) THEN
	491	ierr = 2
	492	RETURN
	493	ENDIF
[1]	494
	495
[3725]	496	SELECT CASE ( igo )
[1]	497
[3725]	498	!
	499	!-- Coding for factor 2
	500	CASE ( 1 )
	501	ia = 1
	502	ib = ia + iink
	503	ja = 1
	504	jb = ja + (2m-la) inc2
[1]	505
[3725]	506	IF ( la /= m ) THEN
[1]	507
[3725]	508	DO l = 1, la
	509	i = ibase
	510	j = jbase
	511	!DIR$ IVDEP
	512	!CDIR NODEP
	513	!OCL NOVREC
	514	DO ijk = 1, lot
	515	c(ja+j) = a(ia+i) + a(ib+i)
	516	c(jb+j) = a(ia+i) - a(ib+i)
	517	i = i + inc3
	518	j = j + inc4
	519	ENDDO
	520	ibase = ibase + inc1
	521	jbase = jbase + inc2
	522	ENDDO
[1]	523
[3725]	524	ja = ja + jink
	525	jink = 2 * jink
	526	jb = jb - jink
	527	ibase = ibase + ijump
	528	ijump = 2 * ijump + iink
[1]	529
[3725]	530	IF ( ja /= jb ) THEN
[1]	531
[3725]	532	DO k = la, kstop, la
	533	kb = k + k
	534	c1 = trigs(kb+1)
	535	s1 = trigs(kb+2)
	536	jbase = 0
	537	DO l = 1, la
	538	i = ibase
	539	j = jbase
	540	!DIR$ IVDEP
	541	!CDIR NODEP
	542	!OCL NOVREC
	543	DO ijk = 1, lot
	544	c(ja+j) = a(ia+i) + (c1a(ib+i)+s1b(ib+i))
	545	c(jb+j) = a(ia+i) - (c1a(ib+i)+s1b(ib+i))
	546	d(ja+j) = (c1b(ib+i)-s1a(ib+i)) + b(ia+i)
	547	d(jb+j) = (c1b(ib+i)-s1a(ib+i)) - b(ia+i)
	548	i = i + inc3
	549	j = j + inc4
	550	ENDDO
	551	ibase = ibase + inc1
	552	jbase = jbase + inc2
	553	ENDDO
[1]	554
[3725]	555	ibase = ibase + ijump
	556	ja = ja + jink
	557	jb = jb - jink
	558	ENDDO
[1]	559
[3725]	560	IF ( ja > jb ) RETURN
[1]	561
[3725]	562	ENDIF
[1]	563
[3725]	564	jbase = 0
	565	DO l = 1, la
	566	i = ibase
	567	j = jbase
	568	!DIR$ IVDEP
	569	!CDIR NODEP
	570	!OCL NOVREC
	571	DO ijk = 1, lot
	572	c(ja+j) = a(ia+i)
	573	d(ja+j) = -a(ib+i)
	574	i = i + inc3
	575	j = j + inc4
[1]	576
[3725]	577	ENDDO
	578	ibase = ibase + inc1
	579	jbase = jbase + inc2
	580	ENDDO
[1]	581
[3725]	582	ELSE
[1]	583
[3725]	584	z = 1.0_wp/REAL(n,KIND=wp)
	585	DO l = 1, la
	586	i = ibase
	587	j = jbase
	588	!DIR$ IVDEP
	589	!CDIR NODEP
	590	!OCL NOVREC
	591	DO ijk = 1, lot
	592	c(ja+j) = z*(a(ia+i)+a(ib+i))
	593	c(jb+j) = z*(a(ia+i)-a(ib+i))
	594	i = i + inc3
	595	j = j + inc4
	596	ENDDO
	597	ibase = ibase + inc1
	598	jbase = jbase + inc2
	599	ENDDO
[1]	600
[3725]	601	ENDIF
[1]	602
[3725]	603	!
	604	!-- Coding for factor 3
	605	CASE ( 2 )
[1]	606
[3725]	607	ia = 1
	608	ib = ia + iink
	609	ic = ib + iink
	610	ja = 1
	611	jb = ja + (2m-la) inc2
	612	jc = jb
[1]	613
[3725]	614	IF ( la /= m ) THEN
[1]	615
[3725]	616	DO l = 1, la
	617	i = ibase
	618	j = jbase
	619	!DIR$ IVDEP
	620	!CDIR NODEP
	621	!OCL NOVREC
	622	DO ijk = 1, lot
	623	c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
	624	c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
	625	d(jb+j) = sin60*(a(ic+i)-a(ib+i))
	626	i = i + inc3
	627	j = j + inc4
	628	ENDDO
	629	ibase = ibase + inc1
	630	jbase = jbase + inc2
	631	ENDDO
[1]	632
[3725]	633	ja = ja + jink
	634	jink = 2 * jink
	635	jb = jb + jink
	636	jc = jc - jink
	637	ibase = ibase + ijump
	638	ijump = 2 * ijump + iink
[1]	639
[3725]	640	IF ( ja /= jc ) THEN
[1]	641
[3725]	642	DO k = la, kstop, la
	643	kb = k + k
	644	kc = kb + kb
	645	c1 = trigs(kb+1)
	646	s1 = trigs(kb+2)
	647	c2 = trigs(kc+1)
	648	s2 = trigs(kc+2)
[1]	649
[3725]	650	jbase = 0
	651	DO l = 1, la
	652	i = ibase
	653	j = jbase
	654	!DIR$ IVDEP
	655	!CDIR NODEP
	656	!OCL NOVREC
	657	DO ijk = 1, lot
	658	a1 = (c1a(ib+i)+s1b(ib+i)) + (c2a(ic+i)+s2b(ic+i))
	659	b1 = (c1b(ib+i)-s1a(ib+i)) + (c2b(ic+i)-s2a(ic+i))
	660	a2 = a(ia+i) - 0.5_wp*a1
	661	b2 = b(ia+i) - 0.5_wp*b1
	662	a3 = sin60((c1a(ib+i)+s1b(ib+i))-(c2a(ic+i)+s2*b(ic+i)))
	663	b3 = sin60((c1b(ib+i)-s1a(ib+i))-(c2b(ic+i)-s2*a(ic+i)))
	664	c(ja+j) = a(ia+i) + a1
	665	d(ja+j) = b(ia+i) + b1
	666	c(jb+j) = a2 + b3
	667	d(jb+j) = b2 - a3
	668	c(jc+j) = a2 - b3
	669	d(jc+j) = -(b2+a3)
	670	i = i + inc3
	671	j = j + inc4
	672	ENDDO
	673	ibase = ibase + inc1
	674	jbase = jbase + inc2
	675	ENDDO
	676	ibase = ibase + ijump
	677	ja = ja + jink
	678	jb = jb + jink
	679	jc = jc - jink
	680	ENDDO
[1]	681
[3725]	682	IF ( ja > jc ) RETURN
[1]	683
[3725]	684	ENDIF
[1]	685
[3725]	686	jbase = 0
	687	DO l = 1, la
	688	i = ibase
	689	j = jbase
	690	!DIR$ IVDEP
	691	!CDIR NODEP
	692	!OCL NOVREC
	693	DO ijk = 1, lot
	694	c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
	695	d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
	696	c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
	697	i = i + inc3
	698	j = j + inc4
	699	ENDDO
	700	ibase = ibase + inc1
	701	jbase = jbase + inc2
	702	ENDDO
[1]	703
[3725]	704	ELSE
[1]	705
[3725]	706	z = 1.0_wp / REAL( n, KIND=wp )
	707	zsin60 = z*sin60
	708	DO l = 1, la
	709	i = ibase
	710	j = jbase
	711	!DIR$ IVDEP
	712	!CDIR NODEP
	713	!OCL NOVREC
	714	DO ijk = 1, lot
	715	c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
	716	c(jb+j) = z(a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))
	717	d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
	718	i = i + inc3
	719	j = j + inc4
	720	ENDDO
	721	ibase = ibase + inc1
	722	jbase = jbase + inc2
	723	ENDDO
[1]	724
[3725]	725	ENDIF
	726
	727	!
	728	!-- Coding for factor 4
	729	CASE ( 3 )
	730
	731	ia = 1
	732	ib = ia + iink
	733	ic = ib + iink
	734	id = ic + iink
	735	ja = 1
	736	jb = ja + (2m-la) inc2
	737	jc = jb + 2minc2
	738	jd = jb
	739
	740	IF ( la /= m ) THEN
	741
	742	DO l = 1, la
	743
	744	i = ibase
	745	j = jbase
	746	!DIR$ IVDEP
	747	!CDIR NODEP
	748	!OCL NOVREC
	749	DO ijk = 1, lot
	750	c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
	751	c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
	752	c(jb+j) = a(ia+i) - a(ic+i)
	753	d(jb+j) = a(id+i) - a(ib+i)
	754	i = i + inc3
	755	j = j + inc4
	756	ENDDO
	757	ibase = ibase + inc1
	758	jbase = jbase + inc2
	759	ENDDO
	760
	761	ja = ja + jink
	762	jink = 2 * jink
	763	jb = jb + jink
	764	jc = jc - jink
	765	jd = jd - jink
	766	ibase = ibase + ijump
	767	ijump = 2 * ijump + iink
	768
	769	IF ( jb /= jc ) THEN
	770
	771	DO k = la, kstop, la
	772	kb = k + k
	773	kc = kb + kb
	774	kd = kc + kb
	775	c1 = trigs(kb+1)
	776	s1 = trigs(kb+2)
	777	c2 = trigs(kc+1)
	778	s2 = trigs(kc+2)
	779	c3 = trigs(kd+1)
	780	s3 = trigs(kd+2)
	781	jbase = 0
	782	DO l = 1, la
	783	i = ibase
	784	j = jbase
	785	!DIR$ IVDEP
	786	!CDIR NODEP
	787	!OCL NOVREC
	788	DO ijk = 1, lot
	789	a0 = a(ia+i) + (c2a(ic+i)+s2b(ic+i))
	790	a2 = a(ia+i) - (c2a(ic+i)+s2b(ic+i))
	791	a1 = (c1a(ib+i)+s1b(ib+i)) + (c3a(id+i)+s3b(id+i))
	792	a3 = (c1a(ib+i)+s1b(ib+i)) - (c3a(id+i)+s3b(id+i))
	793	b0 = b(ia+i) + (c2b(ic+i)-s2a(ic+i))
	794	b2 = b(ia+i) - (c2b(ic+i)-s2a(ic+i))
	795	b1 = (c1b(ib+i)-s1a(ib+i)) + (c3b(id+i)-s3a(id+i))
	796	b3 = (c1b(ib+i)-s1a(ib+i)) - (c3b(id+i)-s3a(id+i))
	797	c(ja+j) = a0 + a1
	798	c(jc+j) = a0 - a1
	799	d(ja+j) = b0 + b1
	800	d(jc+j) = b1 - b0
	801	c(jb+j) = a2 + b3
	802	c(jd+j) = a2 - b3
	803	d(jb+j) = b2 - a3
	804	d(jd+j) = -(b2+a3)
	805	i = i + inc3
	806	j = j + inc4
	807	ENDDO
	808	ibase = ibase + inc1
	809	jbase = jbase + inc2
	810	ENDDO
	811	ibase = ibase + ijump
	812	ja = ja + jink
	813	jb = jb + jink
	814	jc = jc - jink
	815	jd = jd - jink
	816	ENDDO
	817
	818	IF ( jb > jc ) RETURN
	819
	820	ENDIF
	821
	822	sin45 = SQRT( 0.5_wp )
	823	jbase = 0
	824	DO l = 1, la
	825	i = ibase
	826	j = jbase
	827	!DIR$ IVDEP
	828	!CDIR NODEP
	829	!OCL NOVREC
	830	DO ijk = 1, lot
	831	c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
	832	c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
	833	d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
	834	d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
	835	i = i + inc3
	836	j = j + inc4
	837	ENDDO
	838	ibase = ibase + inc1
	839	jbase = jbase + inc2
	840	ENDDO
	841
	842	ELSE
	843
	844	z = 1.0_wp / REAL( n, KIND=wp )
	845	DO l = 1, la
	846	i = ibase
	847	j = jbase
	848	!DIR$ IVDEP
	849	!CDIR NODEP
	850	!OCL NOVREC
	851	DO ijk = 1, lot
	852	c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
	853	c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
	854	c(jb+j) = z*(a(ia+i)-a(ic+i))
	855	d(jb+j) = z*(a(id+i)-a(ib+i))
	856	i = i + inc3
	857	j = j + inc4
	858	ENDDO
	859	ibase = ibase + inc1
	860	jbase = jbase + inc2
	861	ENDDO
	862
	863	ENDIF
	864
	865	!
	866	!-- Coding for factor 5
	867	CASE ( 4 )
	868
	869	ia = 1
	870	ib = ia + iink
	871	ic = ib + iink
	872	id = ic + iink
	873	ie = id + iink
	874	ja = 1
	875	jb = ja + (2m-la) inc2
	876	jc = jb + 2minc2
	877	jd = jc
	878	je = jb
	879
	880	IF ( la /= m ) THEN
	881
	882	DO l = 1, la
	883	i = ibase
	884	j = jbase
	885	!DIR$ IVDEP
	886	!CDIR NODEP
	887	!OCL NOVREC
	888	DO ijk = 1, lot
	889	a1 = a(ib+i) + a(ie+i)
	890	a3 = a(ib+i) - a(ie+i)
	891	a2 = a(ic+i) + a(id+i)
	892	a4 = a(ic+i) - a(id+i)
	893	a5 = a(ia+i) - 0.25_wp*(a1+a2)
	894	a6 = qrt5*(a1-a2)
	895	c(ja+j) = a(ia+i) + (a1+a2)
	896	c(jb+j) = a5 + a6
	897	c(jc+j) = a5 - a6
	898	d(jb+j) = -sin72a3 - sin36a4
	899	d(jc+j) = -sin36a3 + sin72a4
	900	i = i + inc3
	901	j = j + inc4
	902	ENDDO
	903	ibase = ibase + inc1
	904	jbase = jbase + inc2
	905	ENDDO
	906
	907	ja = ja + jink
	908	jink = 2 * jink
	909	jb = jb + jink
	910	jc = jc + jink
	911	jd = jd - jink
	912	je = je - jink
	913	ibase = ibase + ijump
	914	ijump = 2 * ijump + iink
	915
	916	IF ( jb /= jd ) THEN
	917
	918	DO k = la, kstop, la
	919	kb = k + k
	920	kc = kb + kb
	921	kd = kc + kb
	922	ke = kd + kb
	923	c1 = trigs(kb+1)
	924	s1 = trigs(kb+2)
	925	c2 = trigs(kc+1)
	926	s2 = trigs(kc+2)
	927	c3 = trigs(kd+1)
	928	s3 = trigs(kd+2)
	929	c4 = trigs(ke+1)
	930	s4 = trigs(ke+2)
	931	jbase = 0
	932	DO l = 1, la
	933	i = ibase
	934	j = jbase
	935	!DIR$ IVDEP
	936	!CDIR NODEP
	937	!OCL NOVREC
	938	DO ijk = 1, lot
	939	a1 = (c1a(ib+i)+s1b(ib+i)) + (c4a(ie+i)+s4b(ie+i))
	940	a3 = (c1a(ib+i)+s1b(ib+i)) - (c4a(ie+i)+s4b(ie+i))
	941	a2 = (c2a(ic+i)+s2b(ic+i)) + (c3a(id+i)+s3b(id+i))
	942	a4 = (c2a(ic+i)+s2b(ic+i)) - (c3a(id+i)+s3b(id+i))
	943	b1 = (c1b(ib+i)-s1a(ib+i)) + (c4b(ie+i)-s4a(ie+i))
	944	b3 = (c1b(ib+i)-s1a(ib+i)) - (c4b(ie+i)-s4a(ie+i))
	945	b2 = (c2b(ic+i)-s2a(ic+i)) + (c3b(id+i)-s3a(id+i))
	946	b4 = (c2b(ic+i)-s2a(ic+i)) - (c3b(id+i)-s3a(id+i))
	947	a5 = a(ia+i) - 0.25_wp*(a1+a2)
	948	a6 = qrt5*(a1-a2)
	949	b5 = b(ia+i) - 0.25_wp*(b1+b2)
	950	b6 = qrt5*(b1-b2)
	951	a10 = a5 + a6
	952	a20 = a5 - a6
	953	b10 = b5 + b6
	954	b20 = b5 - b6
	955	a11 = sin72b3 + sin36b4
	956	a21 = sin36b3 - sin72b4
	957	b11 = sin72a3 + sin36a4
	958	b21 = sin36a3 - sin72a4
	959	c(ja+j) = a(ia+i) + (a1+a2)
	960	c(jb+j) = a10 + a11
	961	c(je+j) = a10 - a11
	962	c(jc+j) = a20 + a21
	963	c(jd+j) = a20 - a21
	964	d(ja+j) = b(ia+i) + (b1+b2)
	965	d(jb+j) = b10 - b11
	966	d(je+j) = -(b10+b11)
	967	d(jc+j) = b20 - b21
	968	d(jd+j) = -(b20+b21)
	969	i = i + inc3
	970	j = j + inc4
	971	ENDDO
	972	ibase = ibase + inc1
	973	jbase = jbase + inc2
	974	ENDDO
	975	ibase = ibase + ijump
	976	ja = ja + jink
	977	jb = jb + jink
	978	jc = jc + jink
	979	jd = jd - jink
	980	je = je - jink
	981	ENDDO
	982
	983	IF ( jb > jd ) RETURN
	984
	985	ENDIF
	986
	987	jbase = 0
	988	DO l = 1, la
	989	i = ibase
	990	j = jbase
	991	!DIR$ IVDEP
	992	!CDIR NODEP
	993	!OCL NOVREC
	994	DO ijk = 1, lot
	995	a1 = a(ib+i) + a(ie+i)
	996	a3 = a(ib+i) - a(ie+i)
	997	a2 = a(ic+i) + a(id+i)
	998	a4 = a(ic+i) - a(id+i)
	999	a5 = a(ia+i) + 0.25_wp*(a3-a4)
	1000	a6 = qrt5*(a3+a4)
	1001	c(ja+j) = a5 + a6
	1002	c(jb+j) = a5 - a6
	1003	c(jc+j) = a(ia+i) - (a3-a4)
	1004	d(ja+j) = -sin36a1 - sin72a2
	1005	d(jb+j) = -sin72a1 + sin36a2
	1006	i = i + inc3
	1007	j = j + inc4
	1008	ENDDO
	1009	ibase = ibase + inc1
	1010	jbase = jbase + inc2
	1011	ENDDO
	1012
	1013	ELSE
	1014
	1015	z = 1.0_wp / REAL( n, KIND=wp )
	1016	zqrt5 = z * qrt5
	1017	zsin36 = z * sin36
	1018	zsin72 = z * sin72
	1019	DO l = 1, la
	1020	i = ibase
	1021	j = jbase
	1022	!DIR$ IVDEP
	1023	!CDIR NODEP
	1024	!OCL NOVREC
	1025	DO ijk = 1, lot
	1026	a1 = a(ib+i) + a(ie+i)
	1027	a3 = a(ib+i) - a(ie+i)
	1028	a2 = a(ic+i) + a(id+i)
	1029	a4 = a(ic+i) - a(id+i)
	1030	a5 = z(a(ia+i)-0.25_wp(a1+a2))
	1031	a6 = zqrt5*(a1-a2)
	1032	c(ja+j) = z*(a(ia+i)+(a1+a2))
	1033	c(jb+j) = a5 + a6
	1034	c(jc+j) = a5 - a6
	1035	d(jb+j) = -zsin72a3 - zsin36a4
	1036	d(jc+j) = -zsin36a3 + zsin72a4
	1037	i = i + inc3
	1038	j = j + inc4
	1039	ENDDO
	1040	ibase = ibase + inc1
	1041	jbase = jbase + inc2
	1042	ENDDO
	1043
	1044	ENDIF
	1045
	1046	!
	1047	!-- Coding for factor 6
	1048	CASE ( 5 )
	1049
	1050	ia = 1
	1051	ib = ia + iink
	1052	ic = ib + iink
	1053	id = ic + iink
	1054	ie = id + iink
	1055	if = ie + iink
	1056	ja = 1
	1057	jb = ja + (2m-la) inc2
	1058	jc = jb + 2minc2
	1059	jd = jc + 2minc2
	1060	je = jc
	1061	jf = jb
	1062
	1063	IF ( la /= m ) THEN
	1064
	1065	DO l = 1, la
	1066	i = ibase
	1067	j = jbase
	1068	!DIR$ IVDEP
	1069	!CDIR NODEP
	1070	!OCL NOVREC
	1071	DO ijk = 1, lot
	1072	a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
	1073	c(ja+j) = (a(ia+i)+a(id+i)) + a11
	1074	c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
	1075	d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
	1076	a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
	1077	c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
	1078	d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
	1079	c(jd+j) = (a(ia+i)-a(id+i)) + a11
	1080	i = i + inc3
	1081	j = j + inc4
	1082	ENDDO
	1083	ibase = ibase + inc1
	1084	jbase = jbase + inc2
	1085	ENDDO
	1086	ja = ja + jink
	1087	jink = 2 * jink
	1088	jb = jb + jink
	1089	jc = jc + jink
	1090	jd = jd - jink
	1091	je = je - jink
	1092	jf = jf - jink
	1093	ibase = ibase + ijump
	1094	ijump = 2 * ijump + iink
	1095
	1096	IF ( jc /= jd ) THEN
	1097
	1098	DO k = la, kstop, la
	1099	kb = k + k
	1100	kc = kb + kb
	1101	kd = kc + kb
	1102	ke = kd + kb
	1103	kf = ke + kb
	1104	c1 = trigs(kb+1)
	1105	s1 = trigs(kb+2)
	1106	c2 = trigs(kc+1)
	1107	s2 = trigs(kc+2)
	1108	c3 = trigs(kd+1)
	1109	s3 = trigs(kd+2)
	1110	c4 = trigs(ke+1)
	1111	s4 = trigs(ke+2)
	1112	c5 = trigs(kf+1)
	1113	s5 = trigs(kf+2)
	1114	jbase = 0
	1115	DO l = 1, la
	1116	i = ibase
	1117	j = jbase
	1118	!DIR$ IVDEP
	1119	!CDIR NODEP
	1120	!OCL NOVREC
	1121	DO ijk = 1, lot
	1122	a1 = c1a(ib+i) + s1b(ib+i)
	1123	b1 = c1b(ib+i) - s1a(ib+i)
	1124	a2 = c2a(ic+i) + s2b(ic+i)
	1125	b2 = c2b(ic+i) - s2a(ic+i)
	1126	a3 = c3a(id+i) + s3b(id+i)
	1127	b3 = c3b(id+i) - s3a(id+i)
	1128	a4 = c4a(ie+i) + s4b(ie+i)
	1129	b4 = c4b(ie+i) - s4a(ie+i)
	1130	a5 = c5a(if+i) + s5b(if+i)
	1131	b5 = c5b(if+i) - s5a(if+i)
	1132	a11 = (a2+a5) + (a1+a4)
	1133	a20 = (a(ia+i)+a3) - 0.5_wp*a11
	1134	a21 = sin60*((a2+a5)-(a1+a4))
	1135	b11 = (b2+b5) + (b1+b4)
	1136	b20 = (b(ia+i)+b3) - 0.5_wp*b11
	1137	b21 = sin60*((b2+b5)-(b1+b4))
	1138	c(ja+j) = (a(ia+i)+a3) + a11
	1139	d(ja+j) = (b(ia+i)+b3) + b11
	1140	c(jc+j) = a20 - b21
	1141	d(jc+j) = a21 + b20
	1142	c(je+j) = a20 + b21
	1143	d(je+j) = a21 - b20
	1144	a11 = (a2-a5) + (a4-a1)
	1145	a20 = (a(ia+i)-a3) - 0.5_wp*a11
	1146	a21 = sin60*((a4-a1)-(a2-a5))
	1147	b11 = (b5-b2) - (b4-b1)
	1148	b20 = (b3-b(ia+i)) - 0.5_wp*b11
	1149	b21 = sin60*((b5-b2)+(b4-b1))
	1150	c(jb+j) = a20 - b21
	1151	d(jb+j) = a21 - b20
	1152	c(jd+j) = a11 + (a(ia+i)-a3)
	1153	d(jd+j) = b11 + (b3-b(ia+i))
	1154	c(jf+j) = a20 + b21
	1155	d(jf+j) = a21 + b20
	1156	i = i + inc3
	1157	j = j + inc4
	1158	ENDDO
	1159	ibase = ibase + inc1
	1160	jbase = jbase + inc2
	1161	ENDDO
	1162	ibase = ibase + ijump
	1163	ja = ja + jink
	1164	jb = jb + jink
	1165	jc = jc + jink
	1166	jd = jd - jink
	1167	je = je - jink
	1168	jf = jf - jink
	1169	ENDDO
	1170
	1171	IF ( jc > jd ) RETURN
	1172
	1173	ENDIF
	1174
	1175	jbase = 0
	1176	DO l = 1, la
	1177	i = ibase
	1178	j = jbase
	1179	!DIR$ IVDEP
	1180	!CDIR NODEP
	1181	!OCL NOVREC
	1182	DO ijk = 1, lot
	1183	c(ja+j) = (a(ia+i)+0.5_wp(a(ic+i)-a(ie+i))) + sin60(a(ib+i)-a(if+i))
	1184	d(ja+j) = -(a(id+i)+0.5_wp(a(ib+i)+a(if+i))) - sin60(a(ic+i)+a(ie+i))
	1185	c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
	1186	d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
	1187	c(jc+j) = (a(ia+i)+0.5_wp(a(ic+i)-a(ie+i))) - sin60(a(ib+i)-a(if+i))
	1188	d(jc+j) = -(a(id+i)+0.5_wp(a(ib+i)+a(if+i))) + sin60(a(ic+i)+a(ie+i))
	1189	i = i + inc3
	1190	j = j + inc4
	1191	ENDDO
	1192	ibase = ibase + inc1
	1193	jbase = jbase + inc2
	1194	ENDDO
	1195
	1196	ELSE
	1197
	1198	z = 1.0_wp/REAL(n,KIND=wp)
	1199	zsin60 = z*sin60
	1200	DO l = 1, la
	1201	i = ibase
	1202	j = jbase
	1203	!DIR$ IVDEP
	1204	!CDIR NODEP
	1205	!OCL NOVREC
	1206	DO ijk = 1, lot
	1207	a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
	1208	c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
	1209	c(jc+j) = z((a(ia+i)+a(id+i))-0.5_wpa11)
	1210	d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
	1211	a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
	1212	c(jb+j) = z((a(ia+i)-a(id+i))-0.5_wpa11)
	1213	d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
	1214	c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
	1215	i = i + inc3
	1216	j = j + inc4
	1217	ENDDO
	1218	ibase = ibase + inc1
	1219	jbase = jbase + inc2
	1220	ENDDO
	1221
	1222	ENDIF
	1223
	1224	!
	1225	!-- Coding for factor 8
	1226	CASE ( 6 )
	1227
	1228	IF ( la /= m ) THEN
	1229	ierr = 3
	1230	RETURN
	1231	ENDIF
	1232
	1233	ia = 1
	1234	ib = ia + iink
	1235	ic = ib + iink
	1236	id = ic + iink
	1237	ie = id + iink
	1238	if = ie + iink
	1239	ig = if + iink
	1240	ih = ig + iink
	1241	ja = 1
	1242	jb = ja + la * inc2
	1243	jc = jb + 2minc2
	1244	jd = jc + 2minc2
	1245	je = jd + 2minc2
	1246	z = 1.0_wp / REAL( n, KIND=wp )
	1247	zsin45 = z * SQRT( 0.5_wp )
	1248
	1249	DO l = 1, la
	1250	i = ibase
	1251	j = jbase
	1252
	1253	!DIR$ IVDEP
	1254	!CDIR NODEP
	1255	!OCL NOVREC
	1256	DO ijk = 1, lot
	1257	c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i))))
	1258	c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i))))
	1259	c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
	1260	d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
	1261	c(jb+j) = z(a(ia+i)-a(ie+i)) + zsin45((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
	1262	c(jd+j) = z(a(ia+i)-a(ie+i)) - zsin45((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
	1263	d(jb+j) = zsin45((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z(a(ig+i)-a(ic+i))
	1264	d(jd+j) = zsin45((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z(a(ig+i)-a(ic+i))
	1265	i = i + inc3
	1266	j = j + inc4
	1267	ENDDO
	1268	ibase = ibase + inc1
	1269	jbase = jbase + inc2
	1270	ENDDO
	1271
	1272	END SELECT
	1273
	1274	END SUBROUTINE qpassm
	1275
[1682]	1276	!------------------------------------------------------------------------------!
	1277	! Description:
	1278	! ------------
[3725]	1279	!> Same as qpassm, but for backward fft
[1682]	1280	!------------------------------------------------------------------------------!
[3725]	1281	SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
[1]	1282
[1320]	1283	USE kinds
	1284
[1]	1285	IMPLICIT NONE
	1286
[1682]	1287	INTEGER(iwp) :: ierr !<
	1288	INTEGER(iwp) :: ifac !<
	1289	INTEGER(iwp) :: inc1 !<
	1290	INTEGER(iwp) :: inc2 !<
	1291	INTEGER(iwp) :: inc3 !<
	1292	INTEGER(iwp) :: inc4 !<
	1293	INTEGER(iwp) :: la !<
	1294	INTEGER(iwp) :: lot !<
	1295	INTEGER(iwp) :: n !<
[3725]	1296	!
	1297	!-- Arrays are dimensioned with n
[1682]	1298	REAL(wp) :: a(*) !<
	1299	REAL(wp) :: b(*) !<
	1300	REAL(wp) :: c(*) !<
	1301	REAL(wp) :: d(*) !<
	1302	REAL(wp) :: trigs(*) !<
[1]	1303
[1682]	1304	REAL(wp) :: c1 !<
	1305	REAL(wp) :: c2 !<
	1306	REAL(wp) :: c3 !<
	1307	REAL(wp) :: c4 !<
	1308	REAL(wp) :: c5 !<
	1309	REAL(wp) :: qqrt5 !<
	1310	REAL(wp) :: qrt5 !<
	1311	REAL(wp) :: s1 !<
	1312	REAL(wp) :: s2 !<
	1313	REAL(wp) :: s3 !<
	1314	REAL(wp) :: s4 !<
	1315	REAL(wp) :: s5 !<
	1316	REAL(wp) :: sin36 !<
	1317	REAL(wp) :: sin45 !<
	1318	REAL(wp) :: sin60 !<
	1319	REAL(wp) :: sin72 !<
	1320	REAL(wp) :: ssin36 !<
	1321	REAL(wp) :: ssin45 !<
	1322	REAL(wp) :: ssin60 !<
	1323	REAL(wp) :: ssin72 !<
[1]	1324
[1682]	1325	INTEGER(iwp) :: i !<
	1326	INTEGER(iwp) :: ia !<
	1327	INTEGER(iwp) :: ib !<
	1328	INTEGER(iwp) :: ibase !<
	1329	INTEGER(iwp) :: ic !<
	1330	INTEGER(iwp) :: id !<
	1331	INTEGER(iwp) :: ie !<
	1332	INTEGER(iwp) :: if !<
	1333	INTEGER(iwp) :: igo !<
	1334	INTEGER(iwp) :: iink !<
	1335	INTEGER(iwp) :: ijk !<
	1336	INTEGER(iwp) :: j !<
	1337	INTEGER(iwp) :: ja !<
	1338	INTEGER(iwp) :: jb !<
	1339	INTEGER(iwp) :: jbase !<
	1340	INTEGER(iwp) :: jc !<
	1341	INTEGER(iwp) :: jd !<
	1342	INTEGER(iwp) :: je !<
	1343	INTEGER(iwp) :: jf !<
	1344	INTEGER(iwp) :: jg !<
	1345	INTEGER(iwp) :: jh !<
	1346	INTEGER(iwp) :: jink !<
	1347	INTEGER(iwp) :: jump !<
	1348	INTEGER(iwp) :: k !<
	1349	INTEGER(iwp) :: kb !<
	1350	INTEGER(iwp) :: kc !<
	1351	INTEGER(iwp) :: kd !<
	1352	INTEGER(iwp) :: ke !<
	1353	INTEGER(iwp) :: kf !<
	1354	INTEGER(iwp) :: kstop !<
	1355	INTEGER(iwp) :: l !<
	1356	INTEGER(iwp) :: m !<
[1320]	1357
[1682]	1358	REAL(wp) :: a10(nfft) !<
	1359	REAL(wp) :: a11(nfft) !<
	1360	REAL(wp) :: a20(nfft) !<
	1361	REAL(wp) :: a21(nfft) !<
	1362	REAL(wp) :: b10(nfft) !<
	1363	REAL(wp) :: b11(nfft) !<
	1364	REAL(wp) :: b20(nfft) !<
	1365	REAL(wp) :: b21(nfft) !<
[1]	1366
	1367
[3725]	1368	DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
	1369	qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
[1]	1370
	1371
[3725]	1372	ierr = 0
[1]	1373
[3725]	1374	m = n / ifac
	1375	iink = la * inc1
	1376	jink = la * inc2
	1377	jump = (ifac-1) * jink
	1378	kstop = (n-ifac) / (2*ifac)
[1]	1379
[3725]	1380	IF ( lot > nfft ) THEN
	1381	ierr = 1
	1382	RETURN
	1383	ENDIF
[1]	1384	ibase = 0
	1385	jbase = 0
	1386	igo = ifac - 1
[3725]	1387	IF ( igo == 7 ) igo = 6
	1388	IF ( igo < 1 .OR. igo > 6 ) THEN
	1389	ierr = 2
	1390	RETURN
	1391	ENDIF
[1]	1392
	1393
[3725]	1394	SELECT CASE ( igo )
	1395	!
	1396	!-- Coding for factor 2
	1397	CASE ( 1 )
[1]	1398
[3725]	1399	ia = 1
	1400	ib = ia + (2m-la) inc1
	1401	ja = 1
	1402	jb = ja + jink
[1]	1403
[3725]	1404	IF ( la /= m ) THEN
[1]	1405
[3725]	1406	DO l = 1, la
	1407	i = ibase
	1408	j = jbase
	1409	!DIR$ IVDEP
	1410	!CDIR NODEP
	1411	!OCL NOVREC
	1412	DO ijk = 1, lot
	1413	c(ja+j) = a(ia+i) + a(ib+i)
	1414	c(jb+j) = a(ia+i) - a(ib+i)
	1415	i = i + inc3
	1416	j = j + inc4
	1417	ENDDO
	1418	ibase = ibase + inc1
	1419	jbase = jbase + inc2
	1420	ENDDO
[1]	1421
[3725]	1422	ia = ia + iink
	1423	iink = 2 * iink
	1424	ib = ib - iink
	1425	ibase = 0
	1426	jbase = jbase + jump
	1427	jump = 2 * jump + jink
[1]	1428
[3725]	1429	IF ( ia /= ib ) THEN
[1]	1430
[3725]	1431	DO k = la, kstop, la
	1432	kb = k + k
	1433	c1 = trigs(kb+1)
	1434	s1 = trigs(kb+2)
	1435	ibase = 0
	1436	DO l = 1, la
	1437	i = ibase
	1438	j = jbase
	1439	!DIR$ IVDEP
	1440	!CDIR NODEP
	1441	!OCL NOVREC
	1442	DO ijk = 1, lot
	1443	c(ja+j) = a(ia+i) + a(ib+i)
	1444	d(ja+j) = b(ia+i) - b(ib+i)
	1445	c(jb+j) = c1(a(ia+i)-a(ib+i)) - s1(b(ia+i)+b(ib+i))
	1446	d(jb+j) = s1(a(ia+i)-a(ib+i)) + c1(b(ia+i)+b(ib+i))
	1447	i = i + inc3
	1448	j = j + inc4
	1449	ENDDO
	1450	ibase = ibase + inc1
	1451	jbase = jbase + inc2
	1452	ENDDO
	1453	ia = ia + iink
	1454	ib = ib - iink
	1455	jbase = jbase + jump
	1456	ENDDO
[1]	1457
[3725]	1458	IF ( ia > ib ) RETURN
[1]	1459
[3725]	1460	ENDIF
[1]	1461
[3725]	1462	ibase = 0
	1463	DO l = 1, la
	1464	i = ibase
	1465	j = jbase
	1466	!DIR$ IVDEP
	1467	!CDIR NODEP
	1468	!OCL NOVREC
	1469	DO ijk = 1, lot
	1470	c(ja+j) = a(ia+i)
	1471	c(jb+j) = -b(ia+i)
	1472	i = i + inc3
	1473	j = j + inc4
	1474	ENDDO
	1475	ibase = ibase + inc1
	1476	jbase = jbase + inc2
	1477	ENDDO
[1]	1478
[3725]	1479	ELSE
[1]	1480
[3725]	1481	DO l = 1, la
	1482	i = ibase
	1483	j = jbase
	1484	!DIR$ IVDEP
	1485	!CDIR NODEP
	1486	!OCL NOVREC
	1487	DO ijk = 1, lot
	1488	c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
	1489	c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
	1490	i = i + inc3
	1491	j = j + inc4
	1492	ENDDO
	1493	ibase = ibase + inc1
	1494	jbase = jbase + inc2
	1495	ENDDO
[1]	1496
[3725]	1497	ENDIF
[1]	1498
[3725]	1499	!
	1500	!-- Coding for factor 3
	1501	CASE ( 2 )
[1]	1502
[3725]	1503	ia = 1
	1504	ib = ia + (2m-la) inc1
	1505	ic = ib
	1506	ja = 1
	1507	jb = ja + jink
	1508	jc = jb + jink
[1]	1509
[3725]	1510	IF ( la /= m ) THEN
[1]	1511
[3725]	1512	DO l = 1, la
	1513	i = ibase
	1514	j = jbase
	1515	!DIR$ IVDEP
	1516	!CDIR NODEP
	1517	!OCL NOVREC
	1518	DO ijk = 1, lot
	1519	c(ja+j) = a(ia+i) + a(ib+i)
	1520	c(jb+j) = (a(ia+i)-0.5_wpa(ib+i)) - (sin60(b(ib+i)))
	1521	c(jc+j) = (a(ia+i)-0.5_wpa(ib+i)) + (sin60(b(ib+i)))
	1522	i = i + inc3
	1523	j = j + inc4
	1524	ENDDO
	1525	ibase = ibase + inc1
	1526	jbase = jbase + inc2
	1527	ENDDO
	1528	ia = ia + iink
	1529	iink = 2 * iink
	1530	ib = ib + iink
	1531	ic = ic - iink
	1532	jbase = jbase + jump
	1533	jump = 2 * jump + jink
[1]	1534
[3725]	1535	IF ( ia /= ic ) THEN
[1]	1536
[3725]	1537	DO k = la, kstop, la
	1538	kb = k + k
	1539	kc = kb + kb
	1540	c1 = trigs(kb+1)
	1541	s1 = trigs(kb+2)
	1542	c2 = trigs(kc+1)
	1543	s2 = trigs(kc+2)
	1544	ibase = 0
	1545	DO l = 1, la
	1546	i = ibase
	1547	j = jbase
	1548	!DIR$ IVDEP
	1549	!CDIR NODEP
	1550	!OCL NOVREC
	1551	DO ijk = 1, lot
	1552	c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
	1553	d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
	1554	c(jb+j) = c1((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
	1555	- s1((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
	1556	d(jb+j) = s1((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
	1557	+ c1((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
	1558	c(jc+j) = c2((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
	1559	- s2((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
	1560	d(jc+j) = s2((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
	1561	+ c2((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
	1562	i = i + inc3
	1563	j = j + inc4
	1564	ENDDO
	1565	ibase = ibase + inc1
	1566	jbase = jbase + inc2
	1567	ENDDO
	1568	ia = ia + iink
	1569	ib = ib + iink
	1570	ic = ic - iink
	1571	jbase = jbase + jump
	1572	ENDDO
[1]	1573
[3725]	1574	IF ( ia > ic ) RETURN
[1]	1575
[3725]	1576	ENDIF
[1]	1577
[3725]	1578	ibase = 0
	1579	DO l = 1, la
	1580	i = ibase
	1581	j = jbase
	1582	!DIR$ IVDEP
	1583	!CDIR NODEP
	1584	!OCL NOVREC
	1585	DO ijk = 1, lot
	1586	c(ja+j) = a(ia+i) + a(ib+i)
	1587	c(jb+j) = (0.5_wpa(ia+i)-a(ib+i)) - (sin60b(ia+i))
	1588	c(jc+j) = -(0.5_wpa(ia+i)-a(ib+i)) - (sin60b(ia+i))
	1589	i = i + inc3
	1590	j = j + inc4
	1591	ENDDO
	1592	ibase = ibase + inc1
	1593	jbase = jbase + inc2
	1594	ENDDO
[1]	1595
[3725]	1596	ELSE
[1]	1597
[3725]	1598	ssin60 = 2.0_wp * sin60
	1599	DO l = 1, la
	1600	i = ibase
	1601	j = jbase
	1602	!DIR$ IVDEP
	1603	!CDIR NODEP
	1604	!OCL NOVREC
	1605	DO ijk = 1, lot
	1606	c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
	1607	c(jb+j) = (2.0_wpa(ia+i)-a(ib+i)) - (ssin60b(ib+i))
	1608	c(jc+j) = (2.0_wpa(ia+i)-a(ib+i)) + (ssin60b(ib+i))
	1609	i = i + inc3
	1610	j = j + inc4
	1611	ENDDO
	1612	ibase = ibase + inc1
	1613	jbase = jbase + inc2
	1614	ENDDO
[1]	1615
[3725]	1616	ENDIF
[1]	1617
[3725]	1618	!
	1619	!-- Coding for factor 4
	1620	CASE ( 3 )
[1]	1621
[3725]	1622	ia = 1
	1623	ib = ia + (2m-la) inc1
	1624	ic = ib + 2minc1
	1625	id = ib
	1626	ja = 1
	1627	jb = ja + jink
	1628	jc = jb + jink
	1629	jd = jc + jink
[1]	1630
[3725]	1631	IF ( la /= m ) THEN
[1]	1632
[3725]	1633	DO l = 1, la
	1634	i = ibase
	1635	j = jbase
	1636	!DIR$ IVDEP
	1637	!CDIR NODEP
	1638	!OCL NOVREC
	1639	DO ijk = 1, lot
	1640	c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
	1641	c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
	1642	c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
	1643	c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
	1644	i = i + inc3
	1645	j = j + inc4
	1646	ENDDO
	1647	ibase = ibase + inc1
	1648	jbase = jbase + inc2
	1649	ENDDO
	1650	ia = ia + iink
	1651	iink = 2 * iink
	1652	ib = ib + iink
	1653	ic = ic - iink
	1654	id = id - iink
	1655	jbase = jbase + jump
	1656	jump = 2 * jump + jink
[1]	1657
[3725]	1658	IF ( ib /= ic ) THEN
[1]	1659
[3725]	1660	DO k = la, kstop, la
	1661	kb = k + k
	1662	kc = kb + kb
	1663	kd = kc + kb
	1664	c1 = trigs(kb+1)
	1665	s1 = trigs(kb+2)
	1666	c2 = trigs(kc+1)
	1667	s2 = trigs(kc+2)
	1668	c3 = trigs(kd+1)
	1669	s3 = trigs(kd+2)
	1670	ibase = 0
	1671	DO l = 1, la
	1672	i = ibase
	1673	j = jbase
	1674	!DIR$ IVDEP
	1675	!CDIR NODEP
	1676	!OCL NOVREC
	1677	DO ijk = 1, lot
	1678	c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
	1679	d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
	1680	c(jc+j) = c2((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2((b(ia+i)-b(ic+i))&
	1681	-(b(ib+i)-b(id+i)))
	1682	d(jc+j) = s2((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2((b(ia+i)-b(ic+i))&
	1683	-(b(ib+i)-b(id+i)))
	1684	c(jb+j) = c1((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1((b(ia+i)+b(ic+i))&
	1685	+(a(ib+i)-a(id+i)))
	1686	d(jb+j) = s1((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1((b(ia+i)+b(ic+i))&
	1687	+(a(ib+i)-a(id+i)))
	1688	c(jd+j) = c3((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3((b(ia+i)+b(ic+i))&
	1689	-(a(ib+i)-a(id+i)))
	1690	d(jd+j) = s3((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3((b(ia+i)+b(ic+i))&
	1691	-(a(ib+i)-a(id+i)))
	1692	i = i + inc3
	1693	j = j + inc4
	1694	ENDDO
	1695	ibase = ibase + inc1
	1696	jbase = jbase + inc2
	1697	ENDDO
	1698	ia = ia + iink
	1699	ib = ib + iink
	1700	ic = ic - iink
	1701	id = id - iink
	1702	jbase = jbase + jump
	1703	ENDDO
[1]	1704
[3725]	1705	IF ( ib > ic ) RETURN
[1]	1706
[3725]	1707	ENDIF
[1]	1708
[3725]	1709	ibase = 0
	1710	sin45 = SQRT( 0.5_wp )
	1711	DO l = 1, la
	1712	i = ibase
	1713	j = jbase
	1714	!DIR$ IVDEP
	1715	!CDIR NODEP
	1716	!OCL NOVREC
	1717	DO ijk = 1, lot
	1718	c(ja+j) = a(ia+i) + a(ib+i)
	1719	c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
	1720	c(jc+j) = b(ib+i) - b(ia+i)
	1721	c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
	1722	i = i + inc3
	1723	j = j + inc4
	1724	ENDDO
	1725	ibase = ibase + inc1
	1726	jbase = jbase + inc2
	1727	ENDDO
[1]	1728
[3725]	1729	ELSE
[1]	1730
[3725]	1731	DO l = 1, la
	1732	i = ibase
	1733	j = jbase
	1734	!DIR$ IVDEP
	1735	!CDIR NODEP
	1736	!OCL NOVREC
	1737	DO ijk = 1, lot
	1738	c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
	1739	c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
	1740	c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
	1741	c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
	1742	i = i + inc3
	1743	j = j + inc4
	1744	ENDDO
	1745	ibase = ibase + inc1
	1746	jbase = jbase + inc2
	1747	ENDDO
[1]	1748
[3725]	1749	ENDIF
[1]	1750
[3725]	1751	!
	1752	!-- Coding for factor 5
	1753	CASE ( 4 )
	1754
	1755	ia = 1
	1756	ib = ia + (2m-la) inc1
	1757	ic = ib + 2minc1
	1758	id = ic
	1759	ie = ib
	1760	ja = 1
	1761	jb = ja + jink
	1762	jc = jb + jink
	1763	jd = jc + jink
	1764	je = jd + jink
	1765
	1766	IF ( la /= m ) THEN
	1767
	1768	DO l = 1, la
	1769	i = ibase
	1770	j = jbase
	1771	!DIR$ IVDEP
	1772	!CDIR NODEP
	1773	!OCL NOVREC
	1774	DO ijk = 1, lot
	1775	c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
	1776	c(jb+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))+qrt5(a(ib+i)-a(ic+i))) - &
	1777	(sin72b(ib+i)+sin36b(ic+i))
	1778	c(jc+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))-qrt5(a(ib+i)-a(ic+i))) - &
	1779	(sin36b(ib+i)-sin72b(ic+i))
	1780	c(jd+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))-qrt5(a(ib+i)-a(ic+i))) + &
	1781	(sin36b(ib+i)-sin72b(ic+i))
	1782	c(je+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))+qrt5(a(ib+i)-a(ic+i))) + &
	1783	(sin72b(ib+i)+sin36b(ic+i))
	1784	i = i + inc3
	1785	j = j + inc4
	1786	ENDDO
	1787	ibase = ibase + inc1
	1788	jbase = jbase + inc2
	1789	ENDDO
	1790	ia = ia + iink
	1791	iink = 2 * iink
	1792	ib = ib + iink
	1793	ic = ic + iink
	1794	id = id - iink
	1795	ie = ie - iink
	1796	jbase = jbase + jump
	1797	jump = 2 * jump + jink
	1798
	1799	IF ( ib /= id ) THEN
	1800
	1801	DO k = la, kstop, la
	1802	kb = k + k
	1803	kc = kb + kb
	1804	kd = kc + kb
	1805	ke = kd + kb
	1806	c1 = trigs(kb+1)
	1807	s1 = trigs(kb+2)
	1808	c2 = trigs(kc+1)
	1809	s2 = trigs(kc+2)
	1810	c3 = trigs(kd+1)
	1811	s3 = trigs(kd+2)
	1812	c4 = trigs(ke+1)
	1813	s4 = trigs(ke+2)
	1814	ibase = 0
	1815	DO l = 1, la
	1816	i = ibase
	1817	j = jbase
	1818	!DIR$ IVDEP
	1819	!CDIR NODEP
	1820	!OCL NOVREC
	1821	DO ijk = 1, lot
	1822	a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
	1823	qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
	1824	a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
	1825	qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
	1826	b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
	1827	qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
	1828	b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
	1829	qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
	1830	a11(ijk) = sin72(b(ib+i)+b(ie+i)) + sin36(b(ic+i)+b(id+i))
	1831	a21(ijk) = sin36(b(ib+i)+b(ie+i)) - sin72(b(ic+i)+b(id+i))
	1832	b11(ijk) = sin72(a(ib+i)-a(ie+i)) + sin36(a(ic+i)-a(id+i))
	1833	b21(ijk) = sin36(a(ib+i)-a(ie+i)) - sin72(a(ic+i)-a(id+i))
	1834	c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
	1835	d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
	1836	c(jb+j) = c1(a10(ijk)-a11(ijk)) - s1(b10(ijk)+b11(ijk))
	1837	d(jb+j) = s1(a10(ijk)-a11(ijk)) + c1(b10(ijk)+b11(ijk))
	1838	c(je+j) = c4(a10(ijk)+a11(ijk)) - s4(b10(ijk)-b11(ijk))
	1839	d(je+j) = s4(a10(ijk)+a11(ijk)) + c4(b10(ijk)-b11(ijk))
	1840	c(jc+j) = c2(a20(ijk)-a21(ijk)) - s2(b20(ijk)+b21(ijk))
	1841	d(jc+j) = s2(a20(ijk)-a21(ijk)) + c2(b20(ijk)+b21(ijk))
	1842	c(jd+j) = c3(a20(ijk)+a21(ijk)) - s3(b20(ijk)-b21(ijk))
	1843	d(jd+j) = s3(a20(ijk)+a21(ijk)) + c3(b20(ijk)-b21(ijk))
	1844
	1845	i = i + inc3
	1846	j = j + inc4
	1847	ENDDO
	1848	ibase = ibase + inc1
	1849	jbase = jbase + inc2
	1850	ENDDO
	1851	ia = ia + iink
	1852	ib = ib + iink
	1853	ic = ic + iink
	1854	id = id - iink
	1855	ie = ie - iink
	1856	jbase = jbase + jump
	1857	ENDDO
	1858
	1859	IF ( ib > id ) RETURN
	1860
	1861	ENDIF
	1862
	1863	ibase = 0
	1864	DO l = 1, la
	1865	i = ibase
	1866	j = jbase
	1867	!DIR$ IVDEP
	1868	!CDIR NODEP
	1869	!OCL NOVREC
	1870	DO ijk = 1, lot
	1871	c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
	1872	c(jb+j) = (qrt5(a(ia+i)-a(ib+i))+(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
	1873	(sin36b(ia+i)+sin72b(ib+i))
	1874	c(je+j) = -(qrt5(a(ia+i)-a(ib+i))+(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
	1875	(sin36b(ia+i)+sin72b(ib+i))
	1876	c(jc+j) = (qrt5(a(ia+i)-a(ib+i))-(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
	1877	(sin72b(ia+i)-sin36b(ib+i))
	1878	c(jd+j) = -(qrt5(a(ia+i)-a(ib+i))-(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
	1879	(sin72b(ia+i)-sin36b(ib+i))
	1880	i = i + inc3
	1881	j = j + inc4
	1882	ENDDO
	1883	ibase = ibase + inc1
	1884	jbase = jbase + inc2
	1885	ENDDO
	1886
	1887	ELSE
	1888
	1889	qqrt5 = 2.0_wp * qrt5
	1890	ssin36 = 2.0_wp * sin36
	1891	ssin72 = 2.0_wp * sin72
	1892	DO l = 1, la
	1893	i = ibase
	1894	j = jbase
	1895	!DIR$ IVDEP
	1896	!CDIR NODEP
	1897	!OCL NOVREC
	1898	DO ijk = 1, lot
	1899	c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
	1900	c(jb+j) = (2.0_wp(a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i))) &
	1901	- (ssin72b(ib+i)+ssin36b(ic+i))
	1902	c(jc+j) = (2.0_wp(a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i))) &
	1903	- (ssin36b(ib+i)-ssin72b(ic+i))
	1904	c(jd+j) = (2.0_wp(a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i))) &
	1905	+ (ssin36b(ib+i)-ssin72b(ic+i))
	1906	c(je+j) = (2.0_wp(a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i))) &
	1907	+ (ssin72b(ib+i)+ssin36b(ic+i))
	1908	i = i + inc3
	1909	j = j + inc4
	1910	ENDDO
	1911	ibase = ibase + inc1
	1912	jbase = jbase + inc2
	1913	ENDDO
	1914
	1915	ENDIF
	1916
	1917	!
	1918	!-- Coding for factor 6
	1919	CASE ( 5 )
	1920
	1921	ia = 1
	1922	ib = ia + (2m-la) inc1
	1923	ic = ib + 2minc1
	1924	id = ic + 2minc1
	1925	ie = ic
	1926	if = ib
	1927	ja = 1
	1928	jb = ja + jink
	1929	jc = jb + jink
	1930	jd = jc + jink
	1931	je = jd + jink
	1932	jf = je + jink
	1933
	1934	IF ( la /= m ) THEN
	1935
	1936	DO l = 1, la
	1937	i = ibase
	1938	j = jbase
	1939	!DIR$ IVDEP
	1940	!CDIR NODEP
	1941	!OCL NOVREC
	1942	DO ijk = 1, lot
	1943	c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
	1944	c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
	1945	c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp(a(ib+i)-a(ic+i))) - (sin60(b(ib+i)+b(ic+i)))
	1946	c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp(a(ib+i)-a(ic+i))) + (sin60(b(ib+i)+b(ic+i)))
	1947	c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp(a(ib+i)+a(ic+i))) - (sin60(b(ib+i)-b(ic+i)))
	1948	c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp(a(ib+i)+a(ic+i))) + (sin60(b(ib+i)-b(ic+i)))
	1949	i = i + inc3
	1950	j = j + inc4
	1951	ENDDO
	1952	ibase = ibase + inc1
	1953	jbase = jbase + inc2
	1954	ENDDO
	1955	ia = ia + iink
	1956	iink = 2 * iink
	1957	ib = ib + iink
	1958	ic = ic + iink
	1959	id = id - iink
	1960	ie = ie - iink
	1961	if = if - iink
	1962	jbase = jbase + jump
	1963	jump = 2 * jump + jink
	1964
	1965	IF ( ic /= id ) THEN
	1966
	1967	DO k = la, kstop, la
	1968	kb = k + k
	1969	kc = kb + kb
	1970	kd = kc + kb
	1971	ke = kd + kb
	1972	kf = ke + kb
	1973	c1 = trigs(kb+1)
	1974	s1 = trigs(kb+2)
	1975	c2 = trigs(kc+1)
	1976	s2 = trigs(kc+2)
	1977	c3 = trigs(kd+1)
	1978	s3 = trigs(kd+2)
	1979	c4 = trigs(ke+1)
	1980	s4 = trigs(ke+2)
	1981	c5 = trigs(kf+1)
	1982	s5 = trigs(kf+2)
	1983	ibase = 0
	1984	DO l = 1, la
	1985	i = ibase
	1986	j = jbase
	1987	!DIR$ IVDEP
	1988	!CDIR NODEP
	1989	!OCL NOVREC
	1990	DO ijk = 1, lot
	1991	a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
	1992	a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
	1993	a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
	1994	b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
	1995	b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
	1996	b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
	1997
	1998	c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
	1999	d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
	2000	c(jc+j) = c2(a20(ijk)-b21(ijk)) - s2(b20(ijk)+a21(ijk))
	2001	d(jc+j) = s2(a20(ijk)-b21(ijk)) + c2(b20(ijk)+a21(ijk))
	2002	c(je+j) = c4(a20(ijk)+b21(ijk)) - s4(b20(ijk)-a21(ijk))
	2003	d(je+j) = s4(a20(ijk)+b21(ijk)) + c4(b20(ijk)-a21(ijk))
	2004
	2005	a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
	2006	b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
	2007	a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
	2008	a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
	2009	b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
	2010	b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
	2011
	2012	c(jd+j) = c3((a(ia+i)-a(id+i))+a11(ijk)) - s3((b(ia+i)+b(id+i))-b11(ijk))
	2013	d(jd+j) = s3((a(ia+i)-a(id+i))+a11(ijk)) + c3((b(ia+i)+b(id+i))-b11(ijk))
	2014	c(jb+j) = c1(a20(ijk)-b21(ijk)) - s1(b20(ijk)-a21(ijk))
	2015	d(jb+j) = s1(a20(ijk)-b21(ijk)) + c1(b20(ijk)-a21(ijk))
	2016	c(jf+j) = c5(a20(ijk)+b21(ijk)) - s5(b20(ijk)+a21(ijk))
	2017	d(jf+j) = s5(a20(ijk)+b21(ijk)) + c5(b20(ijk)+a21(ijk))
	2018
	2019	i = i + inc3
	2020	j = j + inc4
	2021	ENDDO
	2022	ibase = ibase + inc1
	2023	jbase = jbase + inc2
	2024	ENDDO
	2025	ia = ia + iink
	2026	ib = ib + iink
	2027	ic = ic + iink
	2028	id = id - iink
	2029	ie = ie - iink
	2030	if = if - iink
	2031	jbase = jbase + jump
	2032	ENDDO
	2033
	2034	IF ( ic > id ) RETURN
	2035
	2036	ENDIF
	2037
	2038	ibase = 0
	2039	DO l = 1, la
	2040	i = ibase
	2041	j = jbase
	2042	!DIR$ IVDEP
	2043	!CDIR NODEP
	2044	!OCL NOVREC
	2045	DO ijk = 1, lot
	2046	c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
	2047	c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
	2048	c(jb+j) = (sin60(a(ia+i)-a(ic+i))) - (0.5_wp(b(ia+i)+b(ic+i))+b(ib+i))
	2049	c(jf+j) = -(sin60(a(ia+i)-a(ic+i))) - (0.5_wp(b(ia+i)+b(ic+i))+b(ib+i))
	2050	c(jc+j) = sin60(b(ic+i)-b(ia+i)) + (0.5_wp(a(ia+i)+a(ic+i))-a(ib+i))
	2051	c(je+j) = sin60(b(ic+i)-b(ia+i)) - (0.5_wp(a(ia+i)+a(ic+i))-a(ib+i))
	2052	i = i + inc3
	2053	j = j + inc4
	2054	ENDDO
	2055	ibase = ibase + inc1
	2056	jbase = jbase + inc2
	2057	ENDDO
	2058
	2059	ELSE
	2060
	2061	ssin60 = 2.0_wp * sin60
	2062	DO l = 1, la
	2063	i = ibase
	2064	j = jbase
	2065	!DIR$ IVDEP
	2066	!CDIR NODEP
	2067	!OCL NOVREC
	2068	DO ijk = 1, lot
	2069	c(ja+j) = (2.0_wp(a(ia+i)+a(id+i))) + (2.0_wp(a(ib+i)+a(ic+i)))
	2070	c(jd+j) = (2.0_wp(a(ia+i)-a(id+i))) - (2.0_wp(a(ib+i)-a(ic+i)))
	2071	c(jb+j) = (2.0_wp(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60(b(ib+i)+b(ic+i)))
	2072	c(jf+j) = (2.0_wp(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60(b(ib+i)+b(ic+i)))
	2073	c(jc+j) = (2.0_wp(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60(b(ib+i)-b(ic+i)))
	2074	c(je+j) = (2.0_wp(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60(b(ib+i)-b(ic+i)))
	2075	i = i + inc3
	2076	j = j + inc4
	2077	ENDDO
	2078	ibase = ibase + inc1
	2079	jbase = jbase + inc2
	2080	ENDDO
	2081
	2082	ENDIF
	2083
	2084	!
	2085	!-- Coding for factor 8
	2086	CASE ( 6 )
	2087
	2088	IF ( la /= m ) THEN
	2089	ierr = 3
	2090	RETURN
	2091	ENDIF
	2092
	2093	ia = 1
	2094	ib = ia + la*inc1
	2095	ic = ib + 2lainc1
	2096	id = ic + 2lainc1
	2097	ie = id + 2lainc1
	2098	ja = 1
	2099	jb = ja + jink
	2100	jc = jb + jink
	2101	jd = jc + jink
	2102	je = jd + jink
	2103	jf = je + jink
	2104	jg = jf + jink
	2105	jh = jg + jink
	2106	ssin45 = SQRT( 2.0_wp )
	2107
	2108	DO l = 1, la
	2109	i = ibase
	2110	j = jbase
	2111	!DIR$ IVDEP
	2112	!CDIR NODEP
	2113	!OCL NOVREC
	2114	DO ijk = 1, lot
	2115	c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
	2116	c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
	2117	c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
	2118	c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
	2119	c(jb+j) = 2.0_wp((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i)))
	2120	c(jf+j) = 2.0_wp((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i)))
	2121	c(jd+j) = 2.0_wp((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i)))
	2122	c(jh+j) = 2.0_wp((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i)))
	2123	i = i + inc3
	2124	j = j + inc4
	2125	ENDDO
	2126	ibase = ibase + inc1
	2127	jbase = jbase + inc2
	2128	ENDDO
	2129
	2130	END SELECT
	2131
	2132	END SUBROUTINE rpassm
	2133
[1682]	2134	!------------------------------------------------------------------------------!
	2135	! Description:
	2136	! ------------
	2137	!> Computes factors of n & trigonometric functins required by fft99 & fft991cy
	2138	!> Method: Dimension trigs(n),ifax(1),jfax(10),lfax(7)
	2139	!> subroutine 'set99' - computes factors of n & trigonometric
	2140	!> functins required by fft99 & fft991cy
	2141	!------------------------------------------------------------------------------!
[3725]	2142	SUBROUTINE set99( trigs, ifax, n )
[1]	2143
	2144
[3725]	2145	USE control_parameters, &
[1320]	2146	ONLY: message_string
[1]	2147
[1320]	2148	USE kinds
	2149
[1]	2150	IMPLICIT NONE
	2151
[1682]	2152	INTEGER(iwp) :: n !<
[1]	2153
[1682]	2154	INTEGER(iwp) :: ifax(*) !<
	2155	REAL(wp) :: trigs(*) !<
[1]	2156
[1320]	2157
[1682]	2158	REAL(wp) :: angle !<
	2159	REAL(wp) :: del !<
	2160	INTEGER(iwp) :: i !<
	2161	INTEGER(iwp) :: ifac !<
	2162	INTEGER(iwp) :: ixxx !<
	2163	INTEGER(iwp) :: k !<
	2164	INTEGER(iwp) :: l !<
	2165	INTEGER(iwp) :: nfax !<
	2166	INTEGER(iwp) :: nhl !<
	2167	INTEGER(iwp) :: nil !<
	2168	INTEGER(iwp) :: nu !<
[1]	2169
[1682]	2170	INTEGER(iwp) :: jfax(10) !<
	2171	INTEGER(iwp) :: lfax(7) !<
[1]	2172
[3725]	2173	DATA lfax / 6, 8, 5, 4, 3, 2, 1 /
[1]	2174
	2175
	2176	ixxx = 1
	2177
[3725]	2178	del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND=wp )
[1]	2179	nil = 0
	2180	nhl = (n/2) - 1
[3725]	2181	DO k = nil, nhl
	2182	angle = REAL( k, KIND=wp ) * del
	2183	trigs(2*k+1) = COS( angle )
	2184	trigs(2*k+2) = SIN( angle )
	2185	ENDDO
[1]	2186
[3725]	2187	!
	2188	!-- Find factors of n (8,6,5,4,3,2; only one 8 allowed).
	2189	!-- Look for sixes first, store factors in descending order.
	2190	nu = n
[1]	2191	ifac = 6
[3725]	2192	k = 0
	2193	l = 1
[1]	2194
[3725]	2195	DO
[1]	2196
[3725]	2197	IF ( MOD( nu, ifac ) == 0 ) THEN
	2198	k = k + 1
	2199	jfax(k) = ifac
	2200	IF ( ifac == 8 .AND. k /= 1 ) THEN
	2201	jfax(1) = 8
	2202	jfax(k) = 6
	2203	ENDIF
	2204	nu = nu/ifac
	2205	IF ( nu == 1 ) EXIT
	2206	IF ( ifac /= 8 ) CYCLE
	2207	ENDIF
[1]	2208
[3725]	2209	l = l + 1
	2210	ifac = lfax(l)
	2211	IF (ifac < 2 ) THEN
	2212	message_string = 'number of gridpoints along x or/and y ' // &
	2213	'contain illegal factors' // &
	2214	'&only factors 2, 3, 5 are allowed'
	2215	CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
	2216	RETURN
	2217	ENDIF
	2218
	2219	ENDDO
	2220	!
	2221	!-- Now reverse order of factors
[1]	2222	nfax = k
	2223	ifax(1) = nfax
[3725]	2224	DO i = 1, nfax
[1]	2225	ifax(nfax+2-i) = jfax(i)
[3725]	2226	ENDDO
[1]	2227	ifax(10) = n
	2228
[3725]	2229	END SUBROUTINE set99
	2230
[1]	2231	END MODULE temperton_fft

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |