1*     fortran version of *dgpfa5* -
2*     radix-5 section of self-sorting, in-place,
3*        generalized pfa
4*
5*-------------------------------------------------------------------
6*
7      subroutine dgpfa5f(a,b,trigs,inc,jump,n,mm,lot,isign)
8      double precision a(*), b(*), trigs(*)
9      integer inc, jump, n, mm, lot, isign
10      double precision s, ax, bx, c1, c2, c3
11      double precision t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11
12      double precision u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11
13      double precision co1, co2, co3, co4, si1, si2, si3, si4
14      double precision aja, ajb, ajc, ajd, aje, bjb, bje, bjc
15      double precision bjd, bja, ajf, ajk, bjf, bjk, ajg, ajj
16      double precision ajh, aji, ajl, ajq, bjg, bjj, bjh, bji
17      double precision bjl, bjq, ajo, ajm, ajn, ajr, ajw, bjo
18      double precision bjm, bjn, bjr, bjw, ajt, ajs, ajx, ajp
19      double precision bjt, bjs, bjx, bjp, ajv, ajy, aju, bjv
20      double precision bjy, bju
21      data sin36/0.587785252292473d+0/,
22     *     sin72/0.951056516295154d+0/,
23     *      qrt5/0.559016994374947d+0/
24      data lvr/128/
25*
26*     ***************************************************************
27*     *                                                             *
28*     *  N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. *
29*     *  RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE *
30*     *  (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER.      *
31*     *                                                             *
32*     ***************************************************************
33*
34      n5 = 5 ** mm
35      inq = n / n5
36      jstepx = (n5-n) * inc
37      ninc = n * inc
38      ink = inc * inq
39      mu = mod(inq,5)
40      if (isign.eq.-1) mu = 5 - mu
41*
42      m = mm
43      mh = (m+1)/2
44      s = dfloat(isign)
45      c1 = qrt5
46      c2 = sin72
47      c3 = sin36
48      if (mu.eq.2.or.mu.eq.3) then
49         c1 = -c1
50         c2 = sin36
51         c3 = sin72
52      endif
53      if (mu.eq.3.or.mu.eq.4) c2 = -c2
54      if (mu.eq.2.or.mu.eq.4) c3 = -c3
55*
56      nblox = 1 + (lot-1)/lvr
57      left = lot
58      s = dfloat(isign)
59      istart = 1
60*
61*  loop on blocks of lvr transforms
62*  --------------------------------
63      do 500 nb = 1 , nblox
64*
65      if (left.le.lvr) then
66         nvex = left
67      else if (left.lt.(2*lvr)) then
68         nvex = left/2
69         nvex = nvex + mod(nvex,2)
70      else
71         nvex = lvr
72      endif
73      left = left - nvex
74*
75      la = 1
76*
77*  loop on type I radix-5 passes
78*  -----------------------------
79      do 160 ipass = 1 , mh
80      jstep = (n*inc) / (5*la)
81      jstepl = jstep - ninc
82      kk = 0
83*
84*  loop on k
85*  ---------
86      do 150 k = 0 , jstep-ink , ink
87*
88      if (k.gt.0) then
89      co1 = trigs(kk+1)
90      si1 = s*trigs(kk+2)
91      co2 = trigs(2*kk+1)
92      si2 = s*trigs(2*kk+2)
93      co3 = trigs(3*kk+1)
94      si3 = s*trigs(3*kk+2)
95      co4 = trigs(4*kk+1)
96      si4 = s*trigs(4*kk+2)
97      endif
98*
99*  loop along transform
100*  --------------------
101      do 140 jjj = k , (n-1)*inc , 5*jstep
102      ja = istart + jjj
103*
104*     "transverse" loop
105*     -----------------
106      do 135 nu = 1 , inq
107      jb = ja + jstepl
108      if (jb.lt.istart) jb = jb + ninc
109      jc = jb + jstepl
110      if (jc.lt.istart) jc = jc + ninc
111      jd = jc + jstepl
112      if (jd.lt.istart) jd = jd + ninc
113      je = jd + jstepl
114      if (je.lt.istart) je = je + ninc
115      j = 0
116*
117*  loop across transforms
118*  ----------------------
119      if (k.eq.0) then
120*
121cdir$ ivdep, shortloop
122      do 110 l = 1 , nvex
123      ajb = a(jb+j)
124      aje = a(je+j)
125      t1 = ajb + aje
126      ajc = a(jc+j)
127      ajd = a(jd+j)
128      t2 = ajc + ajd
129      t3 = ajb - aje
130      t4 = ajc - ajd
131      t5 = t1 + t2
132      t6 = c1 * ( t1 - t2 )
133      aja = a(ja+j)
134      t7 = aja - 0.25 * t5
135      a(ja+j) = aja + t5
136      t8 = t7 + t6
137      t9 = t7 - t6
138      t10 = c3 * t3 - c2 * t4
139      t11 = c2 * t3 + c3 * t4
140      bjb = b(jb+j)
141      bje = b(je+j)
142      u1 = bjb + bje
143      bjc = b(jc+j)
144      bjd = b(jd+j)
145      u2 = bjc + bjd
146      u3 = bjb - bje
147      u4 = bjc - bjd
148      u5 = u1 + u2
149      u6 = c1 * ( u1 - u2 )
150      bja = b(ja+j)
151      u7 = bja - 0.25 * u5
152      b(ja+j) = bja + u5
153      u8 = u7 + u6
154      u9 = u7 - u6
155      u10 = c3 * u3 - c2 * u4
156      u11 = c2 * u3 + c3 * u4
157      a(jb+j) = t8 - u11
158      b(jb+j) = u8 + t11
159      a(je+j) = t8 + u11
160      b(je+j) = u8 - t11
161      a(jc+j) = t9 - u10
162      b(jc+j) = u9 + t10
163      a(jd+j) = t9 + u10
164      b(jd+j) = u9 - t10
165      j = j + jump
166  110 continue
167*
168      else
169*
170cdir$ ivdep,shortloop
171      do 130 l = 1 , nvex
172      ajb = a(jb+j)
173      aje = a(je+j)
174      t1 = ajb + aje
175      ajc = a(jc+j)
176      ajd = a(jd+j)
177      t2 = ajc + ajd
178      t3 = ajb - aje
179      t4 = ajc - ajd
180      t5 = t1 + t2
181      t6 = c1 * ( t1 - t2 )
182      aja = a(ja+j)
183      t7 = aja - 0.25 * t5
184      a(ja+j) = aja + t5
185      t8 = t7 + t6
186      t9 = t7 - t6
187      t10 = c3 * t3 - c2 * t4
188      t11 = c2 * t3 + c3 * t4
189      bjb = b(jb+j)
190      bje = b(je+j)
191      u1 = bjb + bje
192      bjc = b(jc+j)
193      bjd = b(jd+j)
194      u2 = bjc + bjd
195      u3 = bjb - bje
196      u4 = bjc - bjd
197      u5 = u1 + u2
198      u6 = c1 * ( u1 - u2 )
199      bja = b(ja+j)
200      u7 = bja - 0.25 * u5
201      b(ja+j) = bja + u5
202      u8 = u7 + u6
203      u9 = u7 - u6
204      u10 = c3 * u3 - c2 * u4
205      u11 = c2 * u3 + c3 * u4
206      a(jb+j) = co1*(t8-u11) - si1*(u8+t11)
207      b(jb+j) = si1*(t8-u11) + co1*(u8+t11)
208      a(je+j) = co4*(t8+u11) - si4*(u8-t11)
209      b(je+j) = si4*(t8+u11) + co4*(u8-t11)
210      a(jc+j) = co2*(t9-u10) - si2*(u9+t10)
211      b(jc+j) = si2*(t9-u10) + co2*(u9+t10)
212      a(jd+j) = co3*(t9+u10) - si3*(u9-t10)
213      b(jd+j) = si3*(t9+u10) + co3*(u9-t10)
214      j = j + jump
215  130 continue
216*
217      endif
218*
219*-----( end of loop across transforms )
220*
221      ja = ja + jstepx
222      if (ja.lt.istart) ja = ja + ninc
223  135 continue
224  140 continue
225*-----( end of loop along transforms )
226      kk = kk + 2*la
227  150 continue
228*-----( end of loop on nonzero k )
229      la = 5*la
230  160 continue
231*-----( end of loop on type I radix-5 passes)
232*
233      if (n.eq.5) go to 490
234*
235*  loop on type II radix-5 passes
236*  ------------------------------
237  400 continue
238*
239      do 480 ipass = mh+1 , m
240      jstep = (n*inc) / (5*la)
241      jstepl = jstep - ninc
242      laincl = la * ink - ninc
243      kk = 0
244*
245*     loop on k
246*     ---------
247      do 470 k = 0 , jstep-ink , ink
248*
249      if (k.gt.0) then
250      co1 = trigs(kk+1)
251      si1 = s*trigs(kk+2)
252      co2 = trigs(2*kk+1)
253      si2 = s*trigs(2*kk+2)
254      co3 = trigs(3*kk+1)
255      si3 = s*trigs(3*kk+2)
256      co4 = trigs(4*kk+1)
257      si4 = s*trigs(4*kk+2)
258      endif
259*
260*  double loop along first transform in block
261*  ------------------------------------------
262      do 460 ll = k , (la-1)*ink , 5*jstep
263*
264      do 450 jjj = ll , (n-1)*inc , 5*la*ink
265      ja = istart + jjj
266*
267*     "transverse" loop
268*     -----------------
269      do 445 nu = 1 , inq
270      jb = ja + jstepl
271      if (jb.lt.istart) jb = jb + ninc
272      jc = jb + jstepl
273      if (jc.lt.istart) jc = jc + ninc
274      jd = jc + jstepl
275      if (jd.lt.istart) jd = jd + ninc
276      je = jd + jstepl
277      if (je.lt.istart) je = je + ninc
278      jf = ja + laincl
279      if (jf.lt.istart) jf = jf + ninc
280      jg = jf + jstepl
281      if (jg.lt.istart) jg = jg + ninc
282      jh = jg + jstepl
283      if (jh.lt.istart) jh = jh + ninc
284      ji = jh + jstepl
285      if (ji.lt.istart) ji = ji + ninc
286      jj = ji + jstepl
287      if (jj.lt.istart) jj = jj + ninc
288      jk = jf + laincl
289      if (jk.lt.istart) jk = jk + ninc
290      jl = jk + jstepl
291      if (jl.lt.istart) jl = jl + ninc
292      jm = jl + jstepl
293      if (jm.lt.istart) jm = jm + ninc
294      jn = jm + jstepl
295      if (jn.lt.istart) jn = jn + ninc
296      jo = jn + jstepl
297      if (jo.lt.istart) jo = jo + ninc
298      jp = jk + laincl
299      if (jp.lt.istart) jp = jp + ninc
300      jq = jp + jstepl
301      if (jq.lt.istart) jq = jq + ninc
302      jr = jq + jstepl
303      if (jr.lt.istart) jr = jr + ninc
304      js = jr + jstepl
305      if (js.lt.istart) js = js + ninc
306      jt = js + jstepl
307      if (jt.lt.istart) jt = jt + ninc
308      ju = jp + laincl
309      if (ju.lt.istart) ju = ju + ninc
310      jv = ju + jstepl
311      if (jv.lt.istart) jv = jv + ninc
312      jw = jv + jstepl
313      if (jw.lt.istart) jw = jw + ninc
314      jx = jw + jstepl
315      if (jx.lt.istart) jx = jx + ninc
316      jy = jx + jstepl
317      if (jy.lt.istart) jy = jy + ninc
318      j = 0
319*
320*  loop across transforms
321*  ----------------------
322      if (k.eq.0) then
323*
324cdir$ ivdep, shortloop
325      do 410 l = 1 , nvex
326      ajb = a(jb+j)
327      aje = a(je+j)
328      t1 = ajb + aje
329      ajc = a(jc+j)
330      ajd = a(jd+j)
331      t2 = ajc + ajd
332      t3 = ajb - aje
333      t4 = ajc - ajd
334      ajf = a(jf+j)
335      ajb =  ajf
336      t5 = t1 + t2
337      t6 = c1 * ( t1 - t2 )
338      aja = a(ja+j)
339      t7 = aja - 0.25 * t5
340      a(ja+j) = aja + t5
341      t8 = t7 + t6
342      t9 = t7 - t6
343      ajk = a(jk+j)
344      ajc =  ajk
345      t10 = c3 * t3 - c2 * t4
346      t11 = c2 * t3 + c3 * t4
347      bjb = b(jb+j)
348      bje = b(je+j)
349      u1 = bjb + bje
350      bjc = b(jc+j)
351      bjd = b(jd+j)
352      u2 = bjc + bjd
353      u3 = bjb - bje
354      u4 = bjc - bjd
355      bjf = b(jf+j)
356      bjb =  bjf
357      u5 = u1 + u2
358      u6 = c1 * ( u1 - u2 )
359      bja = b(ja+j)
360      u7 = bja - 0.25 * u5
361      b(ja+j) = bja + u5
362      u8 = u7 + u6
363      u9 = u7 - u6
364      bjk = b(jk+j)
365      bjc =  bjk
366      u10 = c3 * u3 - c2 * u4
367      u11 = c2 * u3 + c3 * u4
368      a(jf+j) = t8 - u11
369      b(jf+j) = u8 + t11
370      aje =  t8 + u11
371      bje =  u8 - t11
372      a(jk+j) = t9 - u10
373      b(jk+j) = u9 + t10
374      ajd =  t9 + u10
375      bjd =  u9 - t10
376*----------------------
377      ajg = a(jg+j)
378      ajj = a(jj+j)
379      t1 = ajg + ajj
380      ajh = a(jh+j)
381      aji = a(ji+j)
382      t2 = ajh + aji
383      t3 = ajg - ajj
384      t4 = ajh - aji
385      ajl = a(jl+j)
386      ajh =  ajl
387      t5 = t1 + t2
388      t6 = c1 * ( t1 - t2 )
389      t7 = ajb - 0.25 * t5
390      a(jb+j) = ajb + t5
391      t8 = t7 + t6
392      t9 = t7 - t6
393      ajq = a(jq+j)
394      aji =  ajq
395      t10 = c3 * t3 - c2 * t4
396      t11 = c2 * t3 + c3 * t4
397      bjg = b(jg+j)
398      bjj = b(jj+j)
399      u1 = bjg + bjj
400      bjh = b(jh+j)
401      bji = b(ji+j)
402      u2 = bjh + bji
403      u3 = bjg - bjj
404      u4 = bjh - bji
405      bjl = b(jl+j)
406      bjh =  bjl
407      u5 = u1 + u2
408      u6 = c1 * ( u1 - u2 )
409      u7 = bjb - 0.25 * u5
410      b(jb+j) = bjb + u5
411      u8 = u7 + u6
412      u9 = u7 - u6
413      bjq = b(jq+j)
414      bji =  bjq
415      u10 = c3 * u3 - c2 * u4
416      u11 = c2 * u3 + c3 * u4
417      a(jg+j) = t8 - u11
418      b(jg+j) = u8 + t11
419      ajj =  t8 + u11
420      bjj =  u8 - t11
421      a(jl+j) = t9 - u10
422      b(jl+j) = u9 + t10
423      a(jq+j) = t9 + u10
424      b(jq+j) = u9 - t10
425*----------------------
426      ajo = a(jo+j)
427      t1 = ajh + ajo
428      ajm = a(jm+j)
429      ajn = a(jn+j)
430      t2 = ajm + ajn
431      t3 = ajh - ajo
432      t4 = ajm - ajn
433      ajr = a(jr+j)
434      ajn =  ajr
435      t5 = t1 + t2
436      t6 = c1 * ( t1 - t2 )
437      t7 = ajc - 0.25 * t5
438      a(jc+j) = ajc + t5
439      t8 = t7 + t6
440      t9 = t7 - t6
441      ajw = a(jw+j)
442      ajo =  ajw
443      t10 = c3 * t3 - c2 * t4
444      t11 = c2 * t3 + c3 * t4
445      bjo = b(jo+j)
446      u1 = bjh + bjo
447      bjm = b(jm+j)
448      bjn = b(jn+j)
449      u2 = bjm + bjn
450      u3 = bjh - bjo
451      u4 = bjm - bjn
452      bjr = b(jr+j)
453      bjn =  bjr
454      u5 = u1 + u2
455      u6 = c1 * ( u1 - u2 )
456      u7 = bjc - 0.25 * u5
457      b(jc+j) = bjc + u5
458      u8 = u7 + u6
459      u9 = u7 - u6
460      bjw = b(jw+j)
461      bjo =  bjw
462      u10 = c3 * u3 - c2 * u4
463      u11 = c2 * u3 + c3 * u4
464      a(jh+j) = t8 - u11
465      b(jh+j) = u8 + t11
466      a(jw+j) = t8 + u11
467      b(jw+j) = u8 - t11
468      a(jm+j) = t9 - u10
469      b(jm+j) = u9 + t10
470      a(jr+j) = t9 + u10
471      b(jr+j) = u9 - t10
472*----------------------
473      ajt = a(jt+j)
474      t1 = aji + ajt
475      ajs = a(js+j)
476      t2 = ajn + ajs
477      t3 = aji - ajt
478      t4 = ajn - ajs
479      ajx = a(jx+j)
480      ajt =  ajx
481      t5 = t1 + t2
482      t6 = c1 * ( t1 - t2 )
483      ajp = a(jp+j)
484      t7 = ajp - 0.25 * t5
485      ax = ajp + t5
486      t8 = t7 + t6
487      t9 = t7 - t6
488      a(jp+j) = ajd
489      t10 = c3 * t3 - c2 * t4
490      t11 = c2 * t3 + c3 * t4
491      a(jd+j) = ax
492      bjt = b(jt+j)
493      u1 = bji + bjt
494      bjs = b(js+j)
495      u2 = bjn + bjs
496      u3 = bji - bjt
497      u4 = bjn - bjs
498      bjx = b(jx+j)
499      bjt =  bjx
500      u5 = u1 + u2
501      u6 = c1 * ( u1 - u2 )
502      bjp = b(jp+j)
503      u7 = bjp - 0.25 * u5
504      bx = bjp + u5
505      u8 = u7 + u6
506      u9 = u7 - u6
507      b(jp+j) = bjd
508      u10 = c3 * u3 - c2 * u4
509      u11 = c2 * u3 + c3 * u4
510      b(jd+j) = bx
511      a(ji+j) = t8 - u11
512      b(ji+j) = u8 + t11
513      a(jx+j) = t8 + u11
514      b(jx+j) = u8 - t11
515      a(jn+j) = t9 - u10
516      b(jn+j) = u9 + t10
517      a(js+j) = t9 + u10
518      b(js+j) = u9 - t10
519*----------------------
520      ajv = a(jv+j)
521      ajy = a(jy+j)
522      t1 = ajv + ajy
523      t2 = ajo + ajt
524      t3 = ajv - ajy
525      t4 = ajo - ajt
526      a(jv+j) = ajj
527      t5 = t1 + t2
528      t6 = c1 * ( t1 - t2 )
529      aju = a(ju+j)
530      t7 = aju - 0.25 * t5
531      ax = aju + t5
532      t8 = t7 + t6
533      t9 = t7 - t6
534      a(ju+j) = aje
535      t10 = c3 * t3 - c2 * t4
536      t11 = c2 * t3 + c3 * t4
537      a(je+j) = ax
538      bjv = b(jv+j)
539      bjy = b(jy+j)
540      u1 = bjv + bjy
541      u2 = bjo + bjt
542      u3 = bjv - bjy
543      u4 = bjo - bjt
544      b(jv+j) = bjj
545      u5 = u1 + u2
546      u6 = c1 * ( u1 - u2 )
547      bju = b(ju+j)
548      u7 = bju - 0.25 * u5
549      bx = bju + u5
550      u8 = u7 + u6
551      u9 = u7 - u6
552      b(ju+j) = bje
553      u10 = c3 * u3 - c2 * u4
554      u11 = c2 * u3 + c3 * u4
555      b(je+j) = bx
556      a(jj+j) = t8 - u11
557      b(jj+j) = u8 + t11
558      a(jy+j) = t8 + u11
559      b(jy+j) = u8 - t11
560      a(jo+j) = t9 - u10
561      b(jo+j) = u9 + t10
562      a(jt+j) = t9 + u10
563      b(jt+j) = u9 - t10
564      j = j + jump
565  410 continue
566*
567      else
568*
569cdir$ ivdep, shortloop
570      do 440 l = 1 , nvex
571      ajb = a(jb+j)
572      aje = a(je+j)
573      t1 = ajb + aje
574      ajc = a(jc+j)
575      ajd = a(jd+j)
576      t2 = ajc + ajd
577      t3 = ajb - aje
578      t4 = ajc - ajd
579      ajf = a(jf+j)
580      ajb =  ajf
581      t5 = t1 + t2
582      t6 = c1 * ( t1 - t2 )
583      aja = a(ja+j)
584      t7 = aja - 0.25 * t5
585      a(ja+j) = aja + t5
586      t8 = t7 + t6
587      t9 = t7 - t6
588      ajk = a(jk+j)
589      ajc =  ajk
590      t10 = c3 * t3 - c2 * t4
591      t11 = c2 * t3 + c3 * t4
592      bjb = b(jb+j)
593      bje = b(je+j)
594      u1 = bjb + bje
595      bjc = b(jc+j)
596      bjd = b(jd+j)
597      u2 = bjc + bjd
598      u3 = bjb - bje
599      u4 = bjc - bjd
600      bjf = b(jf+j)
601      bjb =  bjf
602      u5 = u1 + u2
603      u6 = c1 * ( u1 - u2 )
604      bja = b(ja+j)
605      u7 = bja - 0.25 * u5
606      b(ja+j) = bja + u5
607      u8 = u7 + u6
608      u9 = u7 - u6
609      bjk = b(jk+j)
610      bjc =  bjk
611      u10 = c3 * u3 - c2 * u4
612      u11 = c2 * u3 + c3 * u4
613      a(jf+j) = co1*(t8-u11) - si1*(u8+t11)
614      b(jf+j) = si1*(t8-u11) + co1*(u8+t11)
615      aje =  co4*(t8+u11) - si4*(u8-t11)
616      bje =  si4*(t8+u11) + co4*(u8-t11)
617      a(jk+j) = co2*(t9-u10) - si2*(u9+t10)
618      b(jk+j) = si2*(t9-u10) + co2*(u9+t10)
619      ajd =  co3*(t9+u10) - si3*(u9-t10)
620      bjd =  si3*(t9+u10) + co3*(u9-t10)
621*----------------------
622      ajg = a(jg+j)
623      ajj = a(jj+j)
624      t1 = ajg + ajj
625      ajh = a(jh+j)
626      aji = a(ji+j)
627      t2 = ajh + aji
628      t3 = ajg - ajj
629      t4 = ajh - aji
630      ajl = a(jl+j)
631      ajh =  ajl
632      t5 = t1 + t2
633      t6 = c1 * ( t1 - t2 )
634      t7 = ajb - 0.25 * t5
635      a(jb+j) = ajb + t5
636      t8 = t7 + t6
637      t9 = t7 - t6
638      ajq = a(jq+j)
639      aji =  ajq
640      t10 = c3 * t3 - c2 * t4
641      t11 = c2 * t3 + c3 * t4
642      bjg = b(jg+j)
643      bjj = b(jj+j)
644      u1 = bjg + bjj
645      bjh = b(jh+j)
646      bji = b(ji+j)
647      u2 = bjh + bji
648      u3 = bjg - bjj
649      u4 = bjh - bji
650      bjl = b(jl+j)
651      bjh =  bjl
652      u5 = u1 + u2
653      u6 = c1 * ( u1 - u2 )
654      u7 = bjb - 0.25 * u5
655      b(jb+j) = bjb + u5
656      u8 = u7 + u6
657      u9 = u7 - u6
658      bjq = b(jq+j)
659      bji =  bjq
660      u10 = c3 * u3 - c2 * u4
661      u11 = c2 * u3 + c3 * u4
662      a(jg+j) = co1*(t8-u11) - si1*(u8+t11)
663      b(jg+j) = si1*(t8-u11) + co1*(u8+t11)
664      ajj =  co4*(t8+u11) - si4*(u8-t11)
665      bjj =  si4*(t8+u11) + co4*(u8-t11)
666      a(jl+j) = co2*(t9-u10) - si2*(u9+t10)
667      b(jl+j) = si2*(t9-u10) + co2*(u9+t10)
668      a(jq+j) = co3*(t9+u10) - si3*(u9-t10)
669      b(jq+j) = si3*(t9+u10) + co3*(u9-t10)
670*----------------------
671      ajo = a(jo+j)
672      t1 = ajh + ajo
673      ajm = a(jm+j)
674      ajn = a(jn+j)
675      t2 = ajm + ajn
676      t3 = ajh - ajo
677      t4 = ajm - ajn
678      ajr = a(jr+j)
679      ajn =  ajr
680      t5 = t1 + t2
681      t6 = c1 * ( t1 - t2 )
682      t7 = ajc - 0.25 * t5
683      a(jc+j) = ajc + t5
684      t8 = t7 + t6
685      t9 = t7 - t6
686      ajw = a(jw+j)
687      ajo =  ajw
688      t10 = c3 * t3 - c2 * t4
689      t11 = c2 * t3 + c3 * t4
690      bjo = b(jo+j)
691      u1 = bjh + bjo
692      bjm = b(jm+j)
693      bjn = b(jn+j)
694      u2 = bjm + bjn
695      u3 = bjh - bjo
696      u4 = bjm - bjn
697      bjr = b(jr+j)
698      bjn =  bjr
699      u5 = u1 + u2
700      u6 = c1 * ( u1 - u2 )
701      u7 = bjc - 0.25 * u5
702      b(jc+j) = bjc + u5
703      u8 = u7 + u6
704      u9 = u7 - u6
705      bjw = b(jw+j)
706      bjo =  bjw
707      u10 = c3 * u3 - c2 * u4
708      u11 = c2 * u3 + c3 * u4
709      a(jh+j) = co1*(t8-u11) - si1*(u8+t11)
710      b(jh+j) = si1*(t8-u11) + co1*(u8+t11)
711      a(jw+j) = co4*(t8+u11) - si4*(u8-t11)
712      b(jw+j) = si4*(t8+u11) + co4*(u8-t11)
713      a(jm+j) = co2*(t9-u10) - si2*(u9+t10)
714      b(jm+j) = si2*(t9-u10) + co2*(u9+t10)
715      a(jr+j) = co3*(t9+u10) - si3*(u9-t10)
716      b(jr+j) = si3*(t9+u10) + co3*(u9-t10)
717*----------------------
718      ajt = a(jt+j)
719      t1 = aji + ajt
720      ajs = a(js+j)
721      t2 = ajn + ajs
722      t3 = aji - ajt
723      t4 = ajn - ajs
724      ajx = a(jx+j)
725      ajt =  ajx
726      t5 = t1 + t2
727      t6 = c1 * ( t1 - t2 )
728      ajp = a(jp+j)
729      t7 = ajp - 0.25 * t5
730      ax = ajp + t5
731      t8 = t7 + t6
732      t9 = t7 - t6
733      a(jp+j) = ajd
734      t10 = c3 * t3 - c2 * t4
735      t11 = c2 * t3 + c3 * t4
736      a(jd+j) = ax
737      bjt = b(jt+j)
738      u1 = bji + bjt
739      bjs = b(js+j)
740      u2 = bjn + bjs
741      u3 = bji - bjt
742      u4 = bjn - bjs
743      bjx = b(jx+j)
744      bjt =  bjx
745      u5 = u1 + u2
746      u6 = c1 * ( u1 - u2 )
747      bjp = b(jp+j)
748      u7 = bjp - 0.25 * u5
749      bx = bjp + u5
750      u8 = u7 + u6
751      u9 = u7 - u6
752      b(jp+j) = bjd
753      u10 = c3 * u3 - c2 * u4
754      u11 = c2 * u3 + c3 * u4
755      b(jd+j) = bx
756      a(ji+j) = co1*(t8-u11) - si1*(u8+t11)
757      b(ji+j) = si1*(t8-u11) + co1*(u8+t11)
758      a(jx+j) = co4*(t8+u11) - si4*(u8-t11)
759      b(jx+j) = si4*(t8+u11) + co4*(u8-t11)
760      a(jn+j) = co2*(t9-u10) - si2*(u9+t10)
761      b(jn+j) = si2*(t9-u10) + co2*(u9+t10)
762      a(js+j) = co3*(t9+u10) - si3*(u9-t10)
763      b(js+j) = si3*(t9+u10) + co3*(u9-t10)
764*----------------------
765      ajv = a(jv+j)
766      ajy = a(jy+j)
767      t1 = ajv + ajy
768      t2 = ajo + ajt
769      t3 = ajv - ajy
770      t4 = ajo - ajt
771      a(jv+j) = ajj
772      t5 = t1 + t2
773      t6 = c1 * ( t1 - t2 )
774      aju = a(ju+j)
775      t7 = aju - 0.25 * t5
776      ax = aju + t5
777      t8 = t7 + t6
778      t9 = t7 - t6
779      a(ju+j) = aje
780      t10 = c3 * t3 - c2 * t4
781      t11 = c2 * t3 + c3 * t4
782      a(je+j) = ax
783      bjv = b(jv+j)
784      bjy = b(jy+j)
785      u1 = bjv + bjy
786      u2 = bjo + bjt
787      u3 = bjv - bjy
788      u4 = bjo - bjt
789      b(jv+j) = bjj
790      u5 = u1 + u2
791      u6 = c1 * ( u1 - u2 )
792      bju = b(ju+j)
793      u7 = bju - 0.25 * u5
794      bx = bju + u5
795      u8 = u7 + u6
796      u9 = u7 - u6
797      b(ju+j) = bje
798      u10 = c3 * u3 - c2 * u4
799      u11 = c2 * u3 + c3 * u4
800      b(je+j) = bx
801      a(jj+j) = co1*(t8-u11) - si1*(u8+t11)
802      b(jj+j) = si1*(t8-u11) + co1*(u8+t11)
803      a(jy+j) = co4*(t8+u11) - si4*(u8-t11)
804      b(jy+j) = si4*(t8+u11) + co4*(u8-t11)
805      a(jo+j) = co2*(t9-u10) - si2*(u9+t10)
806      b(jo+j) = si2*(t9-u10) + co2*(u9+t10)
807      a(jt+j) = co3*(t9+u10) - si3*(u9-t10)
808      b(jt+j) = si3*(t9+u10) + co3*(u9-t10)
809      j = j + jump
810  440 continue
811*
812      endif
813*
814*-----(end of loop across transforms)
815*
816      ja = ja + jstepx
817      if (ja.lt.istart) ja = ja + ninc
818  445 continue
819  450 continue
820  460 continue
821*-----( end of double loop for this k )
822      kk = kk + 2*la
823  470 continue
824*-----( end of loop over values of k )
825      la = 5*la
826  480 continue
827*-----( end of loop on type II radix-5 passes )
828*-----( nvex transforms completed)
829  490 continue
830      istart = istart + nvex * jump
831  500 continue
832*-----( end of loop on blocks of transforms )
833*
834      return
835      end
836