Home

Context Navigation

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

Last change on this file since 4338 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 70.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |