MARLEY (Model of Argon Reaction Low Energy Yields) v2.0.0
A Monte Carlo event generator for tens-of-MeV neutrino interactions
Loading...
Searching...
No Matches
fftpack4.cc
1/*
2 * This file is based largely on the following software distribution:
3 *
4 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
5 *
6 * FFTPACK
7 *
8 * Reference
9 * P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
10 * (G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83.
11 *
12 * http://www.netlib.org/fftpack/
13 *
14 * Updated to single, double, and extended precision,
15 * and translated to ISO-Standard C/C++ (without aliasing)
16 * on 10 October 2005 by Andrew Fernandes <andrew_AT_fernandes.org>
17 *
18 * Version 4 April 1985
19 *
20 * A Package of Fortran Subprograms for the Fast Fourier
21 * Transform of Periodic and other Symmetric Sequences
22 *
23 * by
24 *
25 * Paul N Swarztrauber
26 *
27 * National Center for Atmospheric Research, Boulder, Colorado 80307,
28 *
29 * which is sponsored by the National Science Foundation
30 *
31 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
32 *
33 * There appears to be no explicit license for FFTPACK. However, the
34 * package has been incorporated verbatim into a large number of software
35 * systems over the years with numerous types of license without complaint
36 * from the original author; therefore it would appear
37 * that the code is effectively public domain. If you are in doubt,
38 * however, you will need to contact the author or the National Center
39 * for Atmospheric Research to be sure.
40 *
41 * All the changes from the original FFTPACK to the current file
42 * fall under the following BSD-style open-source license:
43 *
44 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
45 *
46 * Copyright (c) 2005, Andrew Fernandes ([email protected]);
47 * All rights reserved.
48 *
49 * Redistribution and use in source and binary forms, with or without
50 * modification, are permitted provided that the following conditions
51 * are met:
52 *
53 * - Redistributions of source code must retain the above copyright
54 * notice, this list of conditions and the following disclaimer.
55 *
56 * - Redistributions in binary form must reproduce the above copyright
57 * notice, this list of conditions and the following disclaimer in the
58 * documentation and/or other materials provided with the distribution.
59 *
60 * - Neither the name of the North Carolina State University nor the
61 * names of its contributors may be used to endorse or promote products
62 * derived from this software without specific prior written permission.
63 *
64 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
65 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
66 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
67 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
68 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
69 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
70 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
71 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
72 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
73 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
74 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
75 * POSSIBILITY OF SUCH DAMAGE.
76 *
77 */
78#include "builtin/fftpack4/fftpack4.h"
79
80#ifdef __cplusplus
81
82#include <cmath> /* the correct precision will be automatically selected */
83using std::cos;
84using std::sin;
85
86#else /* ! __cplusplus */
87
88#include <math.h> /* you must define/typedef the functions 'cos/cosf/cosl' and 'sin/sinf/sinl' as appropriate */
89/* real_t cos(real_t); */
90/* real_t sin(real_t); */
91
92#endif
93
94static void cfftb1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );
95static void cfftf1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );
96static void cffti1( integer_t *n, real_t *wa, integer_t *ifac );
97static void cosqb1( integer_t *n, real_t *x, real_t *w, real_t *xh, integer_t *ifac );
98static void cosqf1( integer_t *n, real_t *x, real_t *w, real_t *xh, integer_t *ifac );
99static void ezfft1( integer_t *n, real_t *wa, integer_t *ifac );
100static void passb( integer_t *nac, integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );
101static void passb2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );
102static void passb3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );
103static void passb4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );
104static void passb5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );
105static void passf( integer_t *nac, integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );
106static void passf2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );
107static void passf3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );
108static void passf4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );
109static void passf5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );
110static void radb2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );
111static void radb3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );
112static void radb4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );
113static void radb5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );
114static void radbg( integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );
115static void radf2( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1 );
116static void radf3( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2 );
117static void radf4( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3 );
118static void radf5( integer_t *ido, integer_t *l1, real_t *cc, real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3, real_t *wa4 );
119static void radfg( integer_t *ido, integer_t *ip, integer_t *l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch, real_t *ch2, real_t *wa );
120static void rfftb1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );
121static void rfftf1( integer_t *n, real_t *c__, real_t *ch, real_t *wa, integer_t *ifac );
122static void rffti1( integer_t *n, real_t *wa, integer_t *ifac );
123static void sint1( integer_t *n, real_t *war, real_t *was, real_t *xh, real_t *x, integer_t *ifac );
124
125/* Subroutine */ void cfftb(integer_t *n, real_t *c__, real_t *wsave,
126 integer_t *ifac)
127{
128 integer_t iw1;
129
130 /* Parameter adjustments */
131 --ifac;
132 --wsave;
133 --c__;
134
135 /* Function Body */
136 if (*n == 1) {
137 return;
138 }
139 iw1 = *n + *n + 1;
140 cfftb1(n, &c__[1], &wsave[1], &wsave[iw1], &ifac[1]);
141 return;
142} /* cfftb_ */
143
144/* Subroutine */ static void cfftb1(integer_t *n, real_t *c__, real_t *ch,
145 real_t *wa, integer_t *ifac)
146{
147 /* System generated locals */
148 integer_t i__1;
149
150 /* Local variables */
151 integer_t i__, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido,
152 idl1, idot;
153
154 /* Parameter adjustments */
155 --ifac;
156 --wa;
157 --ch;
158 --c__;
159
160 /* Function Body */
161 nf = ifac[2];
162 na = 0;
163 l1 = 1;
164 iw = 1;
165 i__1 = nf;
166 for (k1 = 1; k1 <= i__1; ++k1) {
167 ip = ifac[k1 + 2];
168 l2 = ip * l1;
169 ido = *n / l2;
170 idot = ido + ido;
171 idl1 = idot * l1;
172 if (ip != 4) {
173 goto L103;
174 }
175 ix2 = iw + idot;
176 ix3 = ix2 + idot;
177 if (na != 0) {
178 goto L101;
179 }
180 passb4(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
181 goto L102;
182L101:
183 passb4(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3]);
184L102:
185 na = 1 - na;
186 goto L115;
187L103:
188 if (ip != 2) {
189 goto L106;
190 }
191 if (na != 0) {
192 goto L104;
193 }
194 passb2(&idot, &l1, &c__[1], &ch[1], &wa[iw]);
195 goto L105;
196L104:
197 passb2(&idot, &l1, &ch[1], &c__[1], &wa[iw]);
198L105:
199 na = 1 - na;
200 goto L115;
201L106:
202 if (ip != 3) {
203 goto L109;
204 }
205 ix2 = iw + idot;
206 if (na != 0) {
207 goto L107;
208 }
209 passb3(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2]);
210 goto L108;
211L107:
212 passb3(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2]);
213L108:
214 na = 1 - na;
215 goto L115;
216L109:
217 if (ip != 5) {
218 goto L112;
219 }
220 ix2 = iw + idot;
221 ix3 = ix2 + idot;
222 ix4 = ix3 + idot;
223 if (na != 0) {
224 goto L110;
225 }
226 passb5(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
227 ix4]);
228 goto L111;
229L110:
230 passb5(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
231 ix4]);
232L111:
233 na = 1 - na;
234 goto L115;
235L112:
236 if (na != 0) {
237 goto L113;
238 }
239 passb(&nac, &idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1]
240 , &ch[1], &wa[iw]);
241 goto L114;
242L113:
243 passb(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c__[1],
244 &c__[1], &wa[iw]);
245L114:
246 if (nac != 0) {
247 na = 1 - na;
248 }
249L115:
250 l1 = l2;
251 iw += (ip - 1) * idot;
252/* L116: */
253 }
254 if (na == 0) {
255 return;
256 }
257 n2 = *n + *n;
258 i__1 = n2;
259 for (i__ = 1; i__ <= i__1; ++i__) {
260 c__[i__] = ch[i__];
261/* L117: */
262 }
263 return;
264} /* cfftb1_ */
265
266/* Subroutine */ void cfftf(integer_t *n, real_t *c__, real_t *wsave,
267 integer_t *ifac)
268{
269 integer_t iw1;
270
271 /* Parameter adjustments */
272 --ifac;
273 --wsave;
274 --c__;
275
276 /* Function Body */
277 if (*n == 1) {
278 return;
279 }
280 iw1 = *n + *n + 1;
281 cfftf1(n, &c__[1], &wsave[1], &wsave[iw1], &ifac[1]);
282 return;
283} /* cfftf_ */
284
285/* Subroutine */ static void cfftf1(integer_t *n, real_t *c__, real_t *ch,
286 real_t *wa, integer_t *ifac)
287{
288 /* System generated locals */
289 integer_t i__1;
290
291 /* Local variables */
292 integer_t i__, k1, l1, l2, n2, na, nf, ip, iw, ix2, ix3, ix4, nac, ido,
293 idl1, idot;
294
295 /* Parameter adjustments */
296 --ifac;
297 --wa;
298 --ch;
299 --c__;
300
301 /* Function Body */
302 nf = ifac[2];
303 na = 0;
304 l1 = 1;
305 iw = 1;
306 i__1 = nf;
307 for (k1 = 1; k1 <= i__1; ++k1) {
308 ip = ifac[k1 + 2];
309 l2 = ip * l1;
310 ido = *n / l2;
311 idot = ido + ido;
312 idl1 = idot * l1;
313 if (ip != 4) {
314 goto L103;
315 }
316 ix2 = iw + idot;
317 ix3 = ix2 + idot;
318 if (na != 0) {
319 goto L101;
320 }
321 passf4(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
322 goto L102;
323L101:
324 passf4(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3]);
325L102:
326 na = 1 - na;
327 goto L115;
328L103:
329 if (ip != 2) {
330 goto L106;
331 }
332 if (na != 0) {
333 goto L104;
334 }
335 passf2(&idot, &l1, &c__[1], &ch[1], &wa[iw]);
336 goto L105;
337L104:
338 passf2(&idot, &l1, &ch[1], &c__[1], &wa[iw]);
339L105:
340 na = 1 - na;
341 goto L115;
342L106:
343 if (ip != 3) {
344 goto L109;
345 }
346 ix2 = iw + idot;
347 if (na != 0) {
348 goto L107;
349 }
350 passf3(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2]);
351 goto L108;
352L107:
353 passf3(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2]);
354L108:
355 na = 1 - na;
356 goto L115;
357L109:
358 if (ip != 5) {
359 goto L112;
360 }
361 ix2 = iw + idot;
362 ix3 = ix2 + idot;
363 ix4 = ix3 + idot;
364 if (na != 0) {
365 goto L110;
366 }
367 passf5(&idot, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
368 ix4]);
369 goto L111;
370L110:
371 passf5(&idot, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
372 ix4]);
373L111:
374 na = 1 - na;
375 goto L115;
376L112:
377 if (na != 0) {
378 goto L113;
379 }
380 passf(&nac, &idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1]
381 , &ch[1], &wa[iw]);
382 goto L114;
383L113:
384 passf(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c__[1],
385 &c__[1], &wa[iw]);
386L114:
387 if (nac != 0) {
388 na = 1 - na;
389 }
390L115:
391 l1 = l2;
392 iw += (ip - 1) * idot;
393/* L116: */
394 }
395 if (na == 0) {
396 return;
397 }
398 n2 = *n + *n;
399 i__1 = n2;
400 for (i__ = 1; i__ <= i__1; ++i__) {
401 c__[i__] = ch[i__];
402/* L117: */
403 }
404 return;
405} /* cfftf1_ */
406
407/* Subroutine */ void cffti(integer_t *n, real_t *wsave, integer_t *ifac)
408{
409 integer_t iw1;
410
411 /* Parameter adjustments */
412 --ifac;
413 --wsave;
414
415 /* Function Body */
416 if (*n == 1) {
417 return;
418 }
419 iw1 = *n + *n + 1;
420 cffti1(n, &wsave[iw1], &ifac[1]);
421 return;
422} /* cffti_ */
423
424/* Subroutine */ static void cffti1(integer_t *n, real_t *wa, integer_t *ifac)
425{
426 /* Initialized data */
427
428 static integer_t ntryh[4] = { 3,4,2,5 };
429
430 /* System generated locals */
431 integer_t i__1, i__2, i__3;
432
433 /* Local variables */
434 integer_t i__, j, i1, k1, l1, l2, ib;
435 real_t fi;
436 integer_t ld, ii, nf, ip, nl, nq, nr;
437 real_t arg;
438 integer_t ido, ipm;
439 real_t tpi, argh;
440 integer_t idot, ntry=0;
441 real_t argld;
442
443 /* Parameter adjustments */
444 --ifac;
445 --wa;
446
447 /* Function Body */
448 nl = *n;
449 nf = 0;
450 j = 0;
451L101:
452 ++j;
453 if (j - 4 <= 0) {
454 goto L102;
455 } else {
456 goto L103;
457 }
458L102:
459 ntry = ntryh[j - 1];
460 goto L104;
461L103:
462 ntry += 2;
463L104:
464 nq = nl / ntry;
465 nr = nl - ntry * nq;
466 if (nr != 0) {
467 goto L101;
468 } else {
469 goto L105;
470 }
471L105:
472 ++nf;
473 ifac[nf + 2] = ntry;
474 nl = nq;
475 if (ntry != 2) {
476 goto L107;
477 }
478 if (nf == 1) {
479 goto L107;
480 }
481 i__1 = nf;
482 for (i__ = 2; i__ <= i__1; ++i__) {
483 ib = nf - i__ + 2;
484 ifac[ib + 2] = ifac[ib + 1];
485/* L106: */
486 }
487 ifac[3] = 2;
488L107:
489 if (nl != 1) {
490 goto L104;
491 }
492 ifac[1] = *n;
493 ifac[2] = nf;
494 tpi = REAL_CONSTANT(6.283185307179586476925286766559005768394338798750211619498891846);
495 argh = tpi / (real_t) (*n);
496 i__ = 2;
497 l1 = 1;
498 i__1 = nf;
499 for (k1 = 1; k1 <= i__1; ++k1) {
500 ip = ifac[k1 + 2];
501 ld = 0;
502 l2 = l1 * ip;
503 ido = *n / l2;
504 idot = ido + ido + 2;
505 ipm = ip - 1;
506 i__2 = ipm;
507 for (j = 1; j <= i__2; ++j) {
508 i1 = i__;
509 wa[i__ - 1] = REAL_CONSTANT(1.0);
510 wa[i__] = REAL_CONSTANT(0.0);
511 ld += l1;
512 fi = REAL_CONSTANT(0.0);
513 argld = (real_t) ld * argh;
514 i__3 = idot;
515 for (ii = 4; ii <= i__3; ii += 2) {
516 i__ += 2;
517 fi += REAL_CONSTANT(1.0);
518 arg = fi * argld;
519 wa[i__ - 1] = cos(arg);
520 wa[i__] = sin(arg);
521/* L108: */
522 }
523 if (ip <= 5) {
524 goto L109;
525 }
526 wa[i1 - 1] = wa[i__ - 1];
527 wa[i1] = wa[i__];
528L109:
529 ;
530 }
531 l1 = l2;
532/* L110: */
533 }
534 return;
535} /* cffti1_ */
536
537/* Subroutine */ void cosqb(integer_t *n, real_t *x, real_t *wsave,
538 integer_t *ifac)
539{
540 /* Initialized data */
541
542 static real_t tsqrt2 =
543 REAL_CONSTANT(2.82842712474619009760337744841939615713934375053896146353359476);
544
545 /* System generated locals */
546 integer_t i__1;
547
548 /* Local variables */
549 real_t x1;
550
551 /* Parameter adjustments */
552 --ifac;
553 --wsave;
554 --x;
555
556 /* Function Body */
557 if ((i__1 = *n - 2) < 0) {
558 goto L101;
559 } else if (i__1 == 0) {
560 goto L102;
561 } else {
562 goto L103;
563 }
564L101:
565 x[1] *= REAL_CONSTANT(4.0);
566 return;
567L102:
568 x1 = (x[1] + x[2]) * REAL_CONSTANT(4.0);
569 x[2] = tsqrt2 * (x[1] - x[2]);
570 x[1] = x1;
571 return;
572L103:
573 cosqb1(n, &x[1], &wsave[1], &wsave[*n + 1], &ifac[1]);
574 return;
575} /* cosqb_ */
576
577/* Subroutine */ static void cosqb1(integer_t *n, real_t *x, real_t *w,
578 real_t *xh, integer_t *ifac)
579{
580 /* System generated locals */
581 integer_t i__1;
582
583 /* Local variables */
584 integer_t i__, k, kc, np2, ns2;
585 real_t xim1;
586 integer_t modn;
587
588 /* Parameter adjustments */
589 --ifac;
590 --xh;
591 --w;
592 --x;
593
594 /* Function Body */
595 ns2 = (*n + 1) / 2;
596 np2 = *n + 2;
597 i__1 = *n;
598 for (i__ = 3; i__ <= i__1; i__ += 2) {
599 xim1 = x[i__ - 1] + x[i__];
600 x[i__] -= x[i__ - 1];
601 x[i__ - 1] = xim1;
602/* L101: */
603 }
604 x[1] += x[1];
605 modn = *n % 2;
606 if (modn == 0) {
607 x[*n] += x[*n];
608 }
609 rfftb(n, &x[1], &xh[1], &ifac[1]);
610 i__1 = ns2;
611 for (k = 2; k <= i__1; ++k) {
612 kc = np2 - k;
613 xh[k] = w[k - 1] * x[kc] + w[kc - 1] * x[k];
614 xh[kc] = w[k - 1] * x[k] - w[kc - 1] * x[kc];
615/* L102: */
616 }
617 if (modn == 0) {
618 x[ns2 + 1] = w[ns2] * (x[ns2 + 1] + x[ns2 + 1]);
619 }
620 i__1 = ns2;
621 for (k = 2; k <= i__1; ++k) {
622 kc = np2 - k;
623 x[k] = xh[k] + xh[kc];
624 x[kc] = xh[k] - xh[kc];
625/* L103: */
626 }
627 x[1] += x[1];
628 return;
629} /* cosqb1_ */
630
631/* Subroutine */ void cosqf(integer_t *n, real_t *x, real_t *wsave,
632 integer_t *ifac)
633{
634 /* Initialized data */
635
636 static real_t sqrt2 =
637 REAL_CONSTANT(1.41421356237309504880168872420969807856967187536948073176679738);
638
639 /* System generated locals */
640 integer_t i__1;
641
642 /* Local variables */
643 real_t tsqx;
644
645 /* Parameter adjustments */
646 --ifac;
647 --wsave;
648 --x;
649
650 /* Function Body */
651 if ((i__1 = *n - 2) < 0) {
652 goto L102;
653 } else if (i__1 == 0) {
654 goto L101;
655 } else {
656 goto L103;
657 }
658L101:
659 tsqx = sqrt2 * x[2];
660 x[2] = x[1] - tsqx;
661 x[1] += tsqx;
662L102:
663 return;
664L103:
665 cosqf1(n, &x[1], &wsave[1], &wsave[*n + 1], &ifac[1]);
666 return;
667} /* cosqf_ */
668
669/* Subroutine */ static void cosqf1(integer_t *n, real_t *x, real_t *w,
670 real_t *xh, integer_t *ifac)
671{
672 /* System generated locals */
673 integer_t i__1;
674
675 /* Local variables */
676 integer_t i__, k, kc, np2, ns2;
677 real_t xim1;
678 integer_t modn;
679
680 /* Parameter adjustments */
681 --ifac;
682 --xh;
683 --w;
684 --x;
685
686 /* Function Body */
687 ns2 = (*n + 1) / 2;
688 np2 = *n + 2;
689 i__1 = ns2;
690 for (k = 2; k <= i__1; ++k) {
691 kc = np2 - k;
692 xh[k] = x[k] + x[kc];
693 xh[kc] = x[k] - x[kc];
694/* L101: */
695 }
696 modn = *n % 2;
697 if (modn == 0) {
698 xh[ns2 + 1] = x[ns2 + 1] + x[ns2 + 1];
699 }
700 i__1 = ns2;
701 for (k = 2; k <= i__1; ++k) {
702 kc = np2 - k;
703 x[k] = w[k - 1] * xh[kc] + w[kc - 1] * xh[k];
704 x[kc] = w[k - 1] * xh[k] - w[kc - 1] * xh[kc];
705/* L102: */
706 }
707 if (modn == 0) {
708 x[ns2 + 1] = w[ns2] * xh[ns2 + 1];
709 }
710 rfftf(n, &x[1], &xh[1], &ifac[1]);
711 i__1 = *n;
712 for (i__ = 3; i__ <= i__1; i__ += 2) {
713 xim1 = x[i__ - 1] - x[i__];
714 x[i__] = x[i__ - 1] + x[i__];
715 x[i__ - 1] = xim1;
716/* L103: */
717 }
718 return;
719} /* cosqf1_ */
720
721/* Subroutine */ void cosqi(integer_t *n, real_t *wsave, integer_t *ifac)
722{
723 /* Initialized data */
724
725 static real_t pih =
726 REAL_CONSTANT(1.570796326794896619231321691639751442098584699687529104874722962);
727
728 /* System generated locals */
729 integer_t i__1;
730
731 /* Local variables */
732 integer_t k;
733 real_t fk, dt;
734
735 /* Parameter adjustments */
736 --ifac;
737 --wsave;
738
739 /* Function Body */
740 dt = pih / (real_t) (*n);
741 fk = REAL_CONSTANT(0.0);
742 i__1 = *n;
743 for (k = 1; k <= i__1; ++k) {
744 fk += REAL_CONSTANT(1.0);
745 wsave[k] = cos(fk * dt);
746/* L101: */
747 }
748 rffti(n, &wsave[*n + 1], &ifac[1]);
749 return;
750} /* cosqi_ */
751
752/* Subroutine */ void cost(integer_t *n, real_t *x, real_t *wsave,
753 integer_t *ifac)
754{
755 /* System generated locals */
756 integer_t i__1;
757
758 /* Local variables */
759 integer_t i__, k;
760 real_t c1, t1, t2;
761 integer_t kc;
762 real_t xi;
763 integer_t nm1, np1;
764 real_t x1h;
765 integer_t ns2;
766 real_t tx2, x1p3, xim2;
767 integer_t modn;
768
769 /* Parameter adjustments */
770 --ifac;
771 --wsave;
772 --x;
773
774 /* Function Body */
775 nm1 = *n - 1;
776 np1 = *n + 1;
777 ns2 = *n / 2;
778 if ((i__1 = *n - 2) < 0) {
779 goto L106;
780 } else if (i__1 == 0) {
781 goto L101;
782 } else {
783 goto L102;
784 }
785L101:
786 x1h = x[1] + x[2];
787 x[2] = x[1] - x[2];
788 x[1] = x1h;
789 return;
790L102:
791 if (*n > 3) {
792 goto L103;
793 }
794 x1p3 = x[1] + x[3];
795 tx2 = x[2] + x[2];
796 x[2] = x[1] - x[3];
797 x[1] = x1p3 + tx2;
798 x[3] = x1p3 - tx2;
799 return;
800L103:
801 c1 = x[1] - x[*n];
802 x[1] += x[*n];
803 i__1 = ns2;
804 for (k = 2; k <= i__1; ++k) {
805 kc = np1 - k;
806 t1 = x[k] + x[kc];
807 t2 = x[k] - x[kc];
808 c1 += wsave[kc] * t2;
809 t2 = wsave[k] * t2;
810 x[k] = t1 - t2;
811 x[kc] = t1 + t2;
812/* L104: */
813 }
814 modn = *n % 2;
815 if (modn != 0) {
816 x[ns2 + 1] += x[ns2 + 1];
817 }
818 rfftf(&nm1, &x[1], &wsave[*n + 1], &ifac[1]);
819 xim2 = x[2];
820 x[2] = c1;
821 i__1 = *n;
822 for (i__ = 4; i__ <= i__1; i__ += 2) {
823 xi = x[i__];
824 x[i__] = x[i__ - 2] - x[i__ - 1];
825 x[i__ - 1] = xim2;
826 xim2 = xi;
827/* L105: */
828 }
829 if (modn != 0) {
830 x[*n] = xim2;
831 }
832L106:
833 return;
834} /* cost_ */
835
836/* Subroutine */ void costi(integer_t *n, real_t *wsave, integer_t *ifac)
837{
838 /* Initialized data */
839
840 static real_t pi =
841 REAL_CONSTANT(3.141592653589793238462643383279502884197169399375158209749445923);
842
843 /* System generated locals */
844 integer_t i__1;
845
846 /* Local variables */
847 integer_t k, kc;
848 real_t fk, dt;
849 integer_t nm1, np1, ns2;
850
851 /* Parameter adjustments */
852 --ifac;
853 --wsave;
854
855 /* Function Body */
856 if (*n <= 3) {
857 return;
858 }
859 nm1 = *n - 1;
860 np1 = *n + 1;
861 ns2 = *n / 2;
862 dt = pi / (real_t) nm1;
863 fk = REAL_CONSTANT(0.0);
864 i__1 = ns2;
865 for (k = 2; k <= i__1; ++k) {
866 kc = np1 - k;
867 fk += REAL_CONSTANT(1.0);
868 wsave[k] = sin(fk * dt) * REAL_CONSTANT(2.0);
869 wsave[kc] = cos(fk * dt) * REAL_CONSTANT(2.0);
870/* L101: */
871 }
872 rffti(&nm1, &wsave[*n + 1], &ifac[1]);
873 return;
874} /* costi_ */
875
876/* Subroutine */ static void ezfft1(integer_t *n, real_t *wa, integer_t *ifac)
877{
878 /* Initialized data */
879
880 static integer_t ntryh[4] = { 4,2,3,5 };
881 static real_t tpi =
882 REAL_CONSTANT(6.283185307179586476925286766559005768394338798750211419498891846);
883
884 /* System generated locals */
885 integer_t i__1, i__2, i__3;
886
887 /* Local variables */
888 integer_t i__, j, k1, l1, l2, ib, ii, nf, ip, nl, is, nq, nr;
889 real_t ch1, sh1;
890 integer_t ido, ipm;
891 real_t dch1, ch1h, arg1, dsh1;
892 integer_t nfm1;
893 real_t argh;
894 integer_t ntry=0;
895
896 /* Parameter adjustments */
897 --ifac;
898 --wa;
899
900 /* Function Body */
901 nl = *n;
902 nf = 0;
903 j = 0;
904L101:
905 ++j;
906 if (j - 4 <= 0) {
907 goto L102;
908 } else {
909 goto L103;
910 }
911L102:
912 ntry = ntryh[j - 1];
913 goto L104;
914L103:
915 ntry += 2;
916L104:
917 nq = nl / ntry;
918 nr = nl - ntry * nq;
919 if (nr != 0) {
920 goto L101;
921 } else {
922 goto L105;
923 }
924L105:
925 ++nf;
926 ifac[nf + 2] = ntry;
927 nl = nq;
928 if (ntry != 2) {
929 goto L107;
930 }
931 if (nf == 1) {
932 goto L107;
933 }
934 i__1 = nf;
935 for (i__ = 2; i__ <= i__1; ++i__) {
936 ib = nf - i__ + 2;
937 ifac[ib + 2] = ifac[ib + 1];
938/* L106: */
939 }
940 ifac[3] = 2;
941L107:
942 if (nl != 1) {
943 goto L104;
944 }
945 ifac[1] = *n;
946 ifac[2] = nf;
947 argh = tpi / (real_t) (*n);
948 is = 0;
949 nfm1 = nf - 1;
950 l1 = 1;
951 if (nfm1 == 0) {
952 return;
953 }
954 i__1 = nfm1;
955 for (k1 = 1; k1 <= i__1; ++k1) {
956 ip = ifac[k1 + 2];
957 l2 = l1 * ip;
958 ido = *n / l2;
959 ipm = ip - 1;
960 arg1 = (real_t) l1 * argh;
961 ch1 = REAL_CONSTANT(1.0);
962 sh1 = REAL_CONSTANT(0.0);
963 dch1 = cos(arg1);
964 dsh1 = sin(arg1);
965 i__2 = ipm;
966 for (j = 1; j <= i__2; ++j) {
967 ch1h = dch1 * ch1 - dsh1 * sh1;
968 sh1 = dch1 * sh1 + dsh1 * ch1;
969 ch1 = ch1h;
970 i__ = is + 2;
971 wa[i__ - 1] = ch1;
972 wa[i__] = sh1;
973 if (ido < 5) {
974 goto L109;
975 }
976 i__3 = ido;
977 for (ii = 5; ii <= i__3; ii += 2) {
978 i__ += 2;
979 wa[i__ - 1] = ch1 * wa[i__ - 3] - sh1 * wa[i__ - 2];
980 wa[i__] = ch1 * wa[i__ - 2] + sh1 * wa[i__ - 3];
981/* L108: */
982 }
983L109:
984 is += ido;
985/* L110: */
986 }
987 l1 = l2;
988/* L111: */
989 }
990 return;
991} /* ezfft1_ */
992
993/* Subroutine */ void ezfftb(integer_t *n, real_t *r__, real_t *azero,
994 real_t *a, real_t *b, real_t *wsave, integer_t *ifac)
995{
996 /* System generated locals */
997 integer_t i__1;
998
999 /* Local variables */
1000 integer_t i__, ns2;
1001
1002 /* Parameter adjustments */
1003 --ifac;
1004 --wsave;
1005 --b;
1006 --a;
1007 --r__;
1008
1009 /* Function Body */
1010 if ((i__1 = *n - 2) < 0) {
1011 goto L101;
1012 } else if (i__1 == 0) {
1013 goto L102;
1014 } else {
1015 goto L103;
1016 }
1017L101:
1018 r__[1] = *azero;
1019 return;
1020L102:
1021 r__[1] = *azero + a[1];
1022 r__[2] = *azero - a[1];
1023 return;
1024L103:
1025 ns2 = (*n - 1) / 2;
1026 i__1 = ns2;
1027 for (i__ = 1; i__ <= i__1; ++i__) {
1028 r__[i__ * 2] = a[i__] * REAL_CONSTANT(0.5);
1029 r__[(i__ << 1) + 1] = b[i__] * REAL_CONSTANT(-0.5);
1030/* L104: */
1031 }
1032 r__[1] = *azero;
1033 if (*n % 2 == 0) {
1034 r__[*n] = a[ns2 + 1];
1035 }
1036 rfftb(n, &r__[1], &wsave[*n + 1], &ifac[1]);
1037 return;
1038} /* ezfftb_ */
1039
1040/* Subroutine */ void ezfftf(integer_t *n, real_t *r__, real_t *azero,
1041 real_t *a, real_t *b, real_t *wsave, integer_t *ifac)
1042{
1043 /* System generated locals */
1044 integer_t i__1;
1045
1046 /* Local variables */
1047 integer_t i__;
1048 real_t cf;
1049 integer_t ns2;
1050 real_t cfm;
1051 integer_t ns2m;
1052
1053/* VERSION 3 JUNE 1979 */
1054
1055 /* Parameter adjustments */
1056 --ifac;
1057 --wsave;
1058 --b;
1059 --a;
1060 --r__;
1061
1062 /* Function Body */
1063 if ((i__1 = *n - 2) < 0) {
1064 goto L101;
1065 } else if (i__1 == 0) {
1066 goto L102;
1067 } else {
1068 goto L103;
1069 }
1070L101:
1071 *azero = r__[1];
1072 return;
1073L102:
1074 *azero = (r__[1] + r__[2]) * REAL_CONSTANT(0.5);
1075 a[1] = (r__[1] - r__[2]) * REAL_CONSTANT(0.5);
1076 return;
1077L103:
1078 i__1 = *n;
1079 for (i__ = 1; i__ <= i__1; ++i__) {
1080 wsave[i__] = r__[i__];
1081/* L104: */
1082 }
1083 rfftf(n, &wsave[1], &wsave[*n + 1], &ifac[1]);
1084 cf = 2.0 / (real_t) (*n);
1085 cfm = -cf;
1086 *azero = cf * 0.5 * wsave[1];
1087 ns2 = (*n + 1) / 2;
1088 ns2m = ns2 - 1;
1089 i__1 = ns2m;
1090 for (i__ = 1; i__ <= i__1; ++i__) {
1091 a[i__] = cf * wsave[i__ * 2];
1092 b[i__] = cfm * wsave[(i__ << 1) + 1];
1093/* L105: */
1094 }
1095 if (*n % 2 == 1) {
1096 return;
1097 }
1098 a[ns2] = cf * 0.5 * wsave[*n];
1099 b[ns2] = REAL_CONSTANT(0.0);
1100 return;
1101} /* ezfftf_ */
1102
1103/* Subroutine */ void ezffti(integer_t *n, real_t *wsave, integer_t *ifac)
1104{
1105 /* Parameter adjustments */
1106 --ifac;
1107 --wsave;
1108
1109 /* Function Body */
1110 if (*n == 1) {
1111 return;
1112 }
1113 ezfft1(n, &wsave[(*n << 1) + 1], &ifac[1]);
1114 return;
1115} /* ezffti_ */
1116
1117/* Subroutine */ static void passb(integer_t *nac, integer_t *ido, integer_t *ip, integer_t *
1118 l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2,
1119 real_t *ch, real_t *ch2, real_t *wa)
1120{
1121 /* System generated locals */
1122 integer_t ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
1123 c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset,
1124 i__1, i__2, i__3;
1125
1126 /* Local variables */
1127 integer_t i__, j, k, l, jc, lc, ik;
1128/*
1129 integer_t nt;
1130*/
1131 integer_t idj, idl, inc, idp;
1132 real_t wai, war;
1133 integer_t ipp2, idij, idlj, idot, ipph;
1134
1135 /* Parameter adjustments */
1136 ch_dim1 = *ido;
1137 ch_dim2 = *l1;
1138 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
1139 ch -= ch_offset;
1140 c1_dim1 = *ido;
1141 c1_dim2 = *l1;
1142 c1_offset = 1 + c1_dim1 * (1 + c1_dim2);
1143 c1 -= c1_offset;
1144 cc_dim1 = *ido;
1145 cc_dim2 = *ip;
1146 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
1147 cc -= cc_offset;
1148 ch2_dim1 = *idl1;
1149 ch2_offset = 1 + ch2_dim1;
1150 ch2 -= ch2_offset;
1151 c2_dim1 = *idl1;
1152 c2_offset = 1 + c2_dim1;
1153 c2 -= c2_offset;
1154 --wa;
1155
1156 /* Function Body */
1157 idot = *ido / 2;
1158/*
1159 nt = *ip * *idl1;
1160*/
1161 ipp2 = *ip + 2;
1162 ipph = (*ip + 1) / 2;
1163 idp = *ip * *ido;
1164
1165 if (*ido < *l1) {
1166 goto L106;
1167 }
1168 i__1 = ipph;
1169 for (j = 2; j <= i__1; ++j) {
1170 jc = ipp2 - j;
1171 i__2 = *l1;
1172 for (k = 1; k <= i__2; ++k) {
1173 i__3 = *ido;
1174 for (i__ = 1; i__ <= i__3; ++i__) {
1175 ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1176 cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) *
1177 cc_dim1];
1178 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1179 cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) *
1180 cc_dim1];
1181/* L101: */
1182 }
1183/* L102: */
1184 }
1185/* L103: */
1186 }
1187 i__1 = *l1;
1188 for (k = 1; k <= i__1; ++k) {
1189 i__2 = *ido;
1190 for (i__ = 1; i__ <= i__2; ++i__) {
1191 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) *
1192 cc_dim1];
1193/* L104: */
1194 }
1195/* L105: */
1196 }
1197 goto L112;
1198L106:
1199 i__1 = ipph;
1200 for (j = 2; j <= i__1; ++j) {
1201 jc = ipp2 - j;
1202 i__2 = *ido;
1203 for (i__ = 1; i__ <= i__2; ++i__) {
1204 i__3 = *l1;
1205 for (k = 1; k <= i__3; ++k) {
1206 ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1207 cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) *
1208 cc_dim1];
1209 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1210 cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) *
1211 cc_dim1];
1212/* L107: */
1213 }
1214/* L108: */
1215 }
1216/* L109: */
1217 }
1218 i__1 = *ido;
1219 for (i__ = 1; i__ <= i__1; ++i__) {
1220 i__2 = *l1;
1221 for (k = 1; k <= i__2; ++k) {
1222 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) *
1223 cc_dim1];
1224/* L110: */
1225 }
1226/* L111: */
1227 }
1228L112:
1229 idl = 2 - *ido;
1230 inc = 0;
1231 i__1 = ipph;
1232 for (l = 2; l <= i__1; ++l) {
1233 lc = ipp2 - l;
1234 idl += *ido;
1235 i__2 = *idl1;
1236 for (ik = 1; ik <= i__2; ++ik) {
1237 c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] * ch2[ik
1238 + (ch2_dim1 << 1)];
1239 c2[ik + lc * c2_dim1] = wa[idl] * ch2[ik + *ip * ch2_dim1];
1240/* L113: */
1241 }
1242 idlj = idl;
1243 inc += *ido;
1244 i__2 = ipph;
1245 for (j = 3; j <= i__2; ++j) {
1246 jc = ipp2 - j;
1247 idlj += inc;
1248 if (idlj > idp) {
1249 idlj -= idp;
1250 }
1251 war = wa[idlj - 1];
1252 wai = wa[idlj];
1253 i__3 = *idl1;
1254 for (ik = 1; ik <= i__3; ++ik) {
1255 c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
1256 c2[ik + lc * c2_dim1] += wai * ch2[ik + jc * ch2_dim1];
1257/* L114: */
1258 }
1259/* L115: */
1260 }
1261/* L116: */
1262 }
1263 i__1 = ipph;
1264 for (j = 2; j <= i__1; ++j) {
1265 i__2 = *idl1;
1266 for (ik = 1; ik <= i__2; ++ik) {
1267 ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
1268/* L117: */
1269 }
1270/* L118: */
1271 }
1272 i__1 = ipph;
1273 for (j = 2; j <= i__1; ++j) {
1274 jc = ipp2 - j;
1275 i__2 = *idl1;
1276 for (ik = 2; ik <= i__2; ik += 2) {
1277 ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik +
1278 jc * c2_dim1];
1279 ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik +
1280 jc * c2_dim1];
1281 ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc *
1282 c2_dim1];
1283 ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc *
1284 c2_dim1];
1285/* L119: */
1286 }
1287/* L120: */
1288 }
1289 *nac = 1;
1290 if (*ido == 2) {
1291 return;
1292 }
1293 *nac = 0;
1294 i__1 = *idl1;
1295 for (ik = 1; ik <= i__1; ++ik) {
1296 c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
1297/* L121: */
1298 }
1299 i__1 = *ip;
1300 for (j = 2; j <= i__1; ++j) {
1301 i__2 = *l1;
1302 for (k = 1; k <= i__2; ++k) {
1303 c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) *
1304 ch_dim1 + 1];
1305 c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) *
1306 ch_dim1 + 2];
1307/* L122: */
1308 }
1309/* L123: */
1310 }
1311 if (idot > *l1) {
1312 goto L127;
1313 }
1314 idij = 0;
1315 i__1 = *ip;
1316 for (j = 2; j <= i__1; ++j) {
1317 idij += 2;
1318 i__2 = *ido;
1319 for (i__ = 4; i__ <= i__2; i__ += 2) {
1320 idij += 2;
1321 i__3 = *l1;
1322 for (k = 1; k <= i__3; ++k) {
1323 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[
1324 i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] *
1325 ch[i__ + (k + j * ch_dim2) * ch_dim1];
1326 c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__
1327 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ -
1328 1 + (k + j * ch_dim2) * ch_dim1];
1329/* L124: */
1330 }
1331/* L125: */
1332 }
1333/* L126: */
1334 }
1335 return;
1336L127:
1337 idj = 2 - *ido;
1338 i__1 = *ip;
1339 for (j = 2; j <= i__1; ++j) {
1340 idj += *ido;
1341 i__2 = *l1;
1342 for (k = 1; k <= i__2; ++k) {
1343 idij = idj;
1344 i__3 = *ido;
1345 for (i__ = 4; i__ <= i__3; i__ += 2) {
1346 idij += 2;
1347 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[
1348 i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] *
1349 ch[i__ + (k + j * ch_dim2) * ch_dim1];
1350 c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__
1351 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ -
1352 1 + (k + j * ch_dim2) * ch_dim1];
1353/* L128: */
1354 }
1355/* L129: */
1356 }
1357/* L130: */
1358 }
1359 return;
1360} /* passb_ */
1361
1362/* Subroutine */ static void passb2(integer_t *ido, integer_t *l1, real_t *cc,
1363 real_t *ch, real_t *wa1)
1364{
1365 /* System generated locals */
1366 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
1367
1368 /* Local variables */
1369 integer_t i__, k;
1370 real_t ti2, tr2;
1371
1372 /* Parameter adjustments */
1373 ch_dim1 = *ido;
1374 ch_dim2 = *l1;
1375 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
1376 ch -= ch_offset;
1377 cc_dim1 = *ido;
1378 cc_offset = 1 + cc_dim1 * 3;
1379 cc -= cc_offset;
1380 --wa1;
1381
1382 /* Function Body */
1383 if (*ido > 2) {
1384 goto L102;
1385 }
1386 i__1 = *l1;
1387 for (k = 1; k <= i__1; ++k) {
1388 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] +
1389 cc[((k << 1) + 2) * cc_dim1 + 1];
1390 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1
1391 + 1] - cc[((k << 1) + 2) * cc_dim1 + 1];
1392 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] +
1393 cc[((k << 1) + 2) * cc_dim1 + 2];
1394 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1
1395 + 2] - cc[((k << 1) + 2) * cc_dim1 + 2];
1396/* L101: */
1397 }
1398 return;
1399L102:
1400 i__1 = *l1;
1401 for (k = 1; k <= i__1; ++k) {
1402 i__2 = *ido;
1403 for (i__ = 2; i__ <= i__2; i__ += 2) {
1404 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + ((k << 1) +
1405 1) * cc_dim1] + cc[i__ - 1 + ((k << 1) + 2) * cc_dim1];
1406 tr2 = cc[i__ - 1 + ((k << 1) + 1) * cc_dim1] - cc[i__ - 1 + ((k <<
1407 1) + 2) * cc_dim1];
1408 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) *
1409 cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1];
1410 ti2 = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2)
1411 * cc_dim1];
1412 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * ti2 +
1413 wa1[i__] * tr2;
1414 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * tr2
1415 - wa1[i__] * ti2;
1416/* L103: */
1417 }
1418/* L104: */
1419 }
1420 return;
1421} /* passb2_ */
1422
1423/* Subroutine */ static void passb3(integer_t *ido, integer_t *l1, real_t *cc,
1424 real_t *ch, real_t *wa1, real_t *wa2)
1425{
1426 /* Initialized data */
1427
1428 static real_t taur = REAL_CONSTANT(-0.5);
1429 static real_t taui =
1430 REAL_CONSTANT(0.8660254037844386467637231707529361834710262690519031402790348975);
1431
1432 /* System generated locals */
1433 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
1434
1435 /* Local variables */
1436 integer_t i__, k;
1437 real_t ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
1438
1439 /* Parameter adjustments */
1440 ch_dim1 = *ido;
1441 ch_dim2 = *l1;
1442 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
1443 ch -= ch_offset;
1444 cc_dim1 = *ido;
1445 cc_offset = 1 + (cc_dim1 << 2);
1446 cc -= cc_offset;
1447 --wa1;
1448 --wa2;
1449
1450 /* Function Body */
1451 if (*ido != 2) {
1452 goto L102;
1453 }
1454 i__1 = *l1;
1455 for (k = 1; k <= i__1; ++k) {
1456 tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
1457 cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
1458 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
1459 ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
1460 ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
1461 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
1462 cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] - cc[(k * 3 + 3) *
1463 cc_dim1 + 1]);
1464 ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] - cc[(k * 3 + 3) *
1465 cc_dim1 + 2]);
1466 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
1467 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
1468 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
1469 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
1470/* L101: */
1471 }
1472 return;
1473L102:
1474 i__1 = *l1;
1475 for (k = 1; k <= i__1; ++k) {
1476 i__2 = *ido;
1477 for (i__ = 2; i__ <= i__2; i__ += 2) {
1478 tr2 = cc[i__ - 1 + (k * 3 + 2) * cc_dim1] + cc[i__ - 1 + (k * 3 +
1479 3) * cc_dim1];
1480 cr2 = cc[i__ - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
1481 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + (k * 3 + 1) *
1482 cc_dim1] + tr2;
1483 ti2 = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) *
1484 cc_dim1];
1485 ci2 = cc[i__ + (k * 3 + 1) * cc_dim1] + taur * ti2;
1486 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) *
1487 cc_dim1] + ti2;
1488 cr3 = taui * (cc[i__ - 1 + (k * 3 + 2) * cc_dim1] - cc[i__ - 1 + (
1489 k * 3 + 3) * cc_dim1]);
1490 ci3 = taui * (cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 +
1491 3) * cc_dim1]);
1492 dr2 = cr2 - ci3;
1493 dr3 = cr2 + ci3;
1494 di2 = ci2 + cr3;
1495 di3 = ci2 - cr3;
1496 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * di2 +
1497 wa1[i__] * dr2;
1498 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * dr2
1499 - wa1[i__] * di2;
1500 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * di3 + wa2[
1501 i__] * dr3;
1502 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * dr3 -
1503 wa2[i__] * di3;
1504/* L103: */
1505 }
1506/* L104: */
1507 }
1508 return;
1509} /* passb3_ */
1510
1511/* Subroutine */ static void passb4(integer_t *ido, integer_t *l1, real_t *cc,
1512 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3)
1513{
1514 /* System generated locals */
1515 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
1516
1517 /* Local variables */
1518 integer_t i__, k;
1519 real_t ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
1520 tr3, tr4;
1521
1522 /* Parameter adjustments */
1523 ch_dim1 = *ido;
1524 ch_dim2 = *l1;
1525 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
1526 ch -= ch_offset;
1527 cc_dim1 = *ido;
1528 cc_offset = 1 + cc_dim1 * 5;
1529 cc -= cc_offset;
1530 --wa1;
1531 --wa2;
1532 --wa3;
1533
1534 /* Function Body */
1535 if (*ido != 2) {
1536 goto L102;
1537 }
1538 i__1 = *l1;
1539 for (k = 1; k <= i__1; ++k) {
1540 ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] - cc[((k << 2) + 3) * cc_dim1
1541 + 2];
1542 ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] + cc[((k << 2) + 3) * cc_dim1
1543 + 2];
1544 tr4 = cc[((k << 2) + 4) * cc_dim1 + 2] - cc[((k << 2) + 2) * cc_dim1
1545 + 2];
1546 ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] + cc[((k << 2) + 4) * cc_dim1
1547 + 2];
1548 tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - cc[((k << 2) + 3) * cc_dim1
1549 + 1];
1550 tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + cc[((k << 2) + 3) * cc_dim1
1551 + 1];
1552 ti4 = cc[((k << 2) + 2) * cc_dim1 + 1] - cc[((k << 2) + 4) * cc_dim1
1553 + 1];
1554 tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] + cc[((k << 2) + 4) * cc_dim1
1555 + 1];
1556 ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
1557 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
1558 ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
1559 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
1560 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
1561 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
1562 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
1563 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
1564/* L101: */
1565 }
1566 return;
1567L102:
1568 i__1 = *l1;
1569 for (k = 1; k <= i__1; ++k) {
1570 i__2 = *ido;
1571 for (i__ = 2; i__ <= i__2; i__ += 2) {
1572 ti1 = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3)
1573 * cc_dim1];
1574 ti2 = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3)
1575 * cc_dim1];
1576 ti3 = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4)
1577 * cc_dim1];
1578 tr4 = cc[i__ + ((k << 2) + 4) * cc_dim1] - cc[i__ + ((k << 2) + 2)
1579 * cc_dim1];
1580 tr1 = cc[i__ - 1 + ((k << 2) + 1) * cc_dim1] - cc[i__ - 1 + ((k <<
1581 2) + 3) * cc_dim1];
1582 tr2 = cc[i__ - 1 + ((k << 2) + 1) * cc_dim1] + cc[i__ - 1 + ((k <<
1583 2) + 3) * cc_dim1];
1584 ti4 = cc[i__ - 1 + ((k << 2) + 2) * cc_dim1] - cc[i__ - 1 + ((k <<
1585 2) + 4) * cc_dim1];
1586 tr3 = cc[i__ - 1 + ((k << 2) + 2) * cc_dim1] + cc[i__ - 1 + ((k <<
1587 2) + 4) * cc_dim1];
1588 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
1589 cr3 = tr2 - tr3;
1590 ch[i__ + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
1591 ci3 = ti2 - ti3;
1592 cr2 = tr1 + tr4;
1593 cr4 = tr1 - tr4;
1594 ci2 = ti1 + ti4;
1595 ci4 = ti1 - ti4;
1596 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * cr2
1597 - wa1[i__] * ci2;
1598 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * ci2 +
1599 wa1[i__] * cr2;
1600 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * cr3 -
1601 wa2[i__] * ci3;
1602 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * ci3 + wa2[
1603 i__] * cr3;
1604 ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * cr4
1605 - wa3[i__] * ci4;
1606 ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * ci4 +
1607 wa3[i__] * cr4;
1608/* L103: */
1609 }
1610/* L104: */
1611 }
1612 return;
1613} /* passb4_ */
1614
1615/* Subroutine */ static void passb5(integer_t *ido, integer_t *l1, real_t *cc,
1616 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3,
1617 real_t *wa4)
1618{
1619 /* Initialized data */
1620
1621 static real_t tr11 =
1622 REAL_CONSTANT(0.3090169943749474241022934171828195886015458990288143106772431137);
1623 static real_t ti11 =
1624 REAL_CONSTANT(0.9510565162951535721164393337938214340569863412575022244730564442);
1625 static real_t tr12 =
1626 REAL_CONSTANT(-0.8090169943749474241022934171828190588601545899028814310677431135);
1627 static real_t ti12 =
1628 REAL_CONSTANT(0.5877852522924731291687059546390727685976524376431459107227248076);
1629
1630 /* System generated locals */
1631 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
1632
1633 /* Local variables */
1634 integer_t i__, k;
1635 real_t ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
1636 ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
1637
1638 /* Parameter adjustments */
1639 ch_dim1 = *ido;
1640 ch_dim2 = *l1;
1641 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
1642 ch -= ch_offset;
1643 cc_dim1 = *ido;
1644 cc_offset = 1 + cc_dim1 * 6;
1645 cc -= cc_offset;
1646 --wa1;
1647 --wa2;
1648 --wa3;
1649 --wa4;
1650
1651 /* Function Body */
1652 if (*ido != 2) {
1653 goto L102;
1654 }
1655 i__1 = *l1;
1656 for (k = 1; k <= i__1; ++k) {
1657 ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
1658 ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
1659 ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
1660 ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
1661 tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
1662 tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
1663 tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
1664 tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
1665 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + tr2
1666 + tr3;
1667 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] + ti2
1668 + ti3;
1669 cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
1670 ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
1671 cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
1672 ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
1673 cr5 = ti11 * tr5 + ti12 * tr4;
1674 ci5 = ti11 * ti5 + ti12 * ti4;
1675 cr4 = ti12 * tr5 - ti11 * tr4;
1676 ci4 = ti12 * ti5 - ti11 * ti4;
1677 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
1678 ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
1679 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
1680 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
1681 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
1682 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
1683 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
1684 ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
1685/* L101: */
1686 }
1687 return;
1688L102:
1689 i__1 = *l1;
1690 for (k = 1; k <= i__1; ++k) {
1691 i__2 = *ido;
1692 for (i__ = 2; i__ <= i__2; i__ += 2) {
1693 ti5 = cc[i__ + (k * 5 + 2) * cc_dim1] - cc[i__ + (k * 5 + 5) *
1694 cc_dim1];
1695 ti2 = cc[i__ + (k * 5 + 2) * cc_dim1] + cc[i__ + (k * 5 + 5) *
1696 cc_dim1];
1697 ti4 = cc[i__ + (k * 5 + 3) * cc_dim1] - cc[i__ + (k * 5 + 4) *
1698 cc_dim1];
1699 ti3 = cc[i__ + (k * 5 + 3) * cc_dim1] + cc[i__ + (k * 5 + 4) *
1700 cc_dim1];
1701 tr5 = cc[i__ - 1 + (k * 5 + 2) * cc_dim1] - cc[i__ - 1 + (k * 5 +
1702 5) * cc_dim1];
1703 tr2 = cc[i__ - 1 + (k * 5 + 2) * cc_dim1] + cc[i__ - 1 + (k * 5 +
1704 5) * cc_dim1];
1705 tr4 = cc[i__ - 1 + (k * 5 + 3) * cc_dim1] - cc[i__ - 1 + (k * 5 +
1706 4) * cc_dim1];
1707 tr3 = cc[i__ - 1 + (k * 5 + 3) * cc_dim1] + cc[i__ - 1 + (k * 5 +
1708 4) * cc_dim1];
1709 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + (k * 5 + 1) *
1710 cc_dim1] + tr2 + tr3;
1711 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 5 + 1) *
1712 cc_dim1] + ti2 + ti3;
1713 cr2 = cc[i__ - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 *
1714 tr3;
1715 ci2 = cc[i__ + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
1716 cr3 = cc[i__ - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 *
1717 tr3;
1718 ci3 = cc[i__ + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
1719 cr5 = ti11 * tr5 + ti12 * tr4;
1720 ci5 = ti11 * ti5 + ti12 * ti4;
1721 cr4 = ti12 * tr5 - ti11 * tr4;
1722 ci4 = ti12 * ti5 - ti11 * ti4;
1723 dr3 = cr3 - ci4;
1724 dr4 = cr3 + ci4;
1725 di3 = ci3 + cr4;
1726 di4 = ci3 - cr4;
1727 dr5 = cr2 + ci5;
1728 dr2 = cr2 - ci5;
1729 di5 = ci2 - cr5;
1730 di2 = ci2 + cr5;
1731 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * dr2
1732 - wa1[i__] * di2;
1733 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * di2 +
1734 wa1[i__] * dr2;
1735 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * dr3 -
1736 wa2[i__] * di3;
1737 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * di3 + wa2[
1738 i__] * dr3;
1739 ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * dr4
1740 - wa3[i__] * di4;
1741 ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * di4 +
1742 wa3[i__] * dr4;
1743 ch[i__ - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i__ - 1] * dr5 -
1744 wa4[i__] * di5;
1745 ch[i__ + (k + ch_dim2 * 5) * ch_dim1] = wa4[i__ - 1] * di5 + wa4[
1746 i__] * dr5;
1747/* L103: */
1748 }
1749/* L104: */
1750 }
1751 return;
1752} /* passb5_ */
1753
1754/* Subroutine */ static void passf(integer_t *nac, integer_t *ido, integer_t *ip, integer_t *
1755 l1, integer_t *idl1, real_t *cc, real_t *c1, real_t *c2,
1756 real_t *ch, real_t *ch2, real_t *wa)
1757{
1758 /* System generated locals */
1759 integer_t ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
1760 c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset,
1761 i__1, i__2, i__3;
1762
1763 /* Local variables */
1764 integer_t i__, j, k, l, jc, lc, ik;
1765/*
1766 integer_t nt;
1767*/
1768 integer_t idj, idl, inc, idp;
1769 real_t wai, war;
1770 integer_t ipp2, idij, idlj, idot, ipph;
1771
1772 /* Parameter adjustments */
1773 ch_dim1 = *ido;
1774 ch_dim2 = *l1;
1775 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
1776 ch -= ch_offset;
1777 c1_dim1 = *ido;
1778 c1_dim2 = *l1;
1779 c1_offset = 1 + c1_dim1 * (1 + c1_dim2);
1780 c1 -= c1_offset;
1781 cc_dim1 = *ido;
1782 cc_dim2 = *ip;
1783 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
1784 cc -= cc_offset;
1785 ch2_dim1 = *idl1;
1786 ch2_offset = 1 + ch2_dim1;
1787 ch2 -= ch2_offset;
1788 c2_dim1 = *idl1;
1789 c2_offset = 1 + c2_dim1;
1790 c2 -= c2_offset;
1791 --wa;
1792
1793 /* Function Body */
1794 idot = *ido / 2;
1795/*
1796 nt = *ip * *idl1;
1797*/
1798 ipp2 = *ip + 2;
1799 ipph = (*ip + 1) / 2;
1800 idp = *ip * *ido;
1801
1802 if (*ido < *l1) {
1803 goto L106;
1804 }
1805 i__1 = ipph;
1806 for (j = 2; j <= i__1; ++j) {
1807 jc = ipp2 - j;
1808 i__2 = *l1;
1809 for (k = 1; k <= i__2; ++k) {
1810 i__3 = *ido;
1811 for (i__ = 1; i__ <= i__3; ++i__) {
1812 ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1813 cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) *
1814 cc_dim1];
1815 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1816 cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) *
1817 cc_dim1];
1818/* L101: */
1819 }
1820/* L102: */
1821 }
1822/* L103: */
1823 }
1824 i__1 = *l1;
1825 for (k = 1; k <= i__1; ++k) {
1826 i__2 = *ido;
1827 for (i__ = 1; i__ <= i__2; ++i__) {
1828 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) *
1829 cc_dim1];
1830/* L104: */
1831 }
1832/* L105: */
1833 }
1834 goto L112;
1835L106:
1836 i__1 = ipph;
1837 for (j = 2; j <= i__1; ++j) {
1838 jc = ipp2 - j;
1839 i__2 = *ido;
1840 for (i__ = 1; i__ <= i__2; ++i__) {
1841 i__3 = *l1;
1842 for (k = 1; k <= i__3; ++k) {
1843 ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1844 cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) *
1845 cc_dim1];
1846 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k *
1847 cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) *
1848 cc_dim1];
1849/* L107: */
1850 }
1851/* L108: */
1852 }
1853/* L109: */
1854 }
1855 i__1 = *ido;
1856 for (i__ = 1; i__ <= i__1; ++i__) {
1857 i__2 = *l1;
1858 for (k = 1; k <= i__2; ++k) {
1859 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) *
1860 cc_dim1];
1861/* L110: */
1862 }
1863/* L111: */
1864 }
1865L112:
1866 idl = 2 - *ido;
1867 inc = 0;
1868 i__1 = ipph;
1869 for (l = 2; l <= i__1; ++l) {
1870 lc = ipp2 - l;
1871 idl += *ido;
1872 i__2 = *idl1;
1873 for (ik = 1; ik <= i__2; ++ik) {
1874 c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] * ch2[ik
1875 + (ch2_dim1 << 1)];
1876 c2[ik + lc * c2_dim1] = -wa[idl] * ch2[ik + *ip * ch2_dim1];
1877/* L113: */
1878 }
1879 idlj = idl;
1880 inc += *ido;
1881 i__2 = ipph;
1882 for (j = 3; j <= i__2; ++j) {
1883 jc = ipp2 - j;
1884 idlj += inc;
1885 if (idlj > idp) {
1886 idlj -= idp;
1887 }
1888 war = wa[idlj - 1];
1889 wai = wa[idlj];
1890 i__3 = *idl1;
1891 for (ik = 1; ik <= i__3; ++ik) {
1892 c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
1893 c2[ik + lc * c2_dim1] -= wai * ch2[ik + jc * ch2_dim1];
1894/* L114: */
1895 }
1896/* L115: */
1897 }
1898/* L116: */
1899 }
1900 i__1 = ipph;
1901 for (j = 2; j <= i__1; ++j) {
1902 i__2 = *idl1;
1903 for (ik = 1; ik <= i__2; ++ik) {
1904 ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
1905/* L117: */
1906 }
1907/* L118: */
1908 }
1909 i__1 = ipph;
1910 for (j = 2; j <= i__1; ++j) {
1911 jc = ipp2 - j;
1912 i__2 = *idl1;
1913 for (ik = 2; ik <= i__2; ik += 2) {
1914 ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik +
1915 jc * c2_dim1];
1916 ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik +
1917 jc * c2_dim1];
1918 ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc *
1919 c2_dim1];
1920 ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc *
1921 c2_dim1];
1922/* L119: */
1923 }
1924/* L120: */
1925 }
1926 *nac = 1;
1927 if (*ido == 2) {
1928 return;
1929 }
1930 *nac = 0;
1931 i__1 = *idl1;
1932 for (ik = 1; ik <= i__1; ++ik) {
1933 c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
1934/* L121: */
1935 }
1936 i__1 = *ip;
1937 for (j = 2; j <= i__1; ++j) {
1938 i__2 = *l1;
1939 for (k = 1; k <= i__2; ++k) {
1940 c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) *
1941 ch_dim1 + 1];
1942 c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) *
1943 ch_dim1 + 2];
1944/* L122: */
1945 }
1946/* L123: */
1947 }
1948 if (idot > *l1) {
1949 goto L127;
1950 }
1951 idij = 0;
1952 i__1 = *ip;
1953 for (j = 2; j <= i__1; ++j) {
1954 idij += 2;
1955 i__2 = *ido;
1956 for (i__ = 4; i__ <= i__2; i__ += 2) {
1957 idij += 2;
1958 i__3 = *l1;
1959 for (k = 1; k <= i__3; ++k) {
1960 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[
1961 i__ - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] *
1962 ch[i__ + (k + j * ch_dim2) * ch_dim1];
1963 c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__
1964 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ -
1965 1 + (k + j * ch_dim2) * ch_dim1];
1966/* L124: */
1967 }
1968/* L125: */
1969 }
1970/* L126: */
1971 }
1972 return;
1973L127:
1974 idj = 2 - *ido;
1975 i__1 = *ip;
1976 for (j = 2; j <= i__1; ++j) {
1977 idj += *ido;
1978 i__2 = *l1;
1979 for (k = 1; k <= i__2; ++k) {
1980 idij = idj;
1981 i__3 = *ido;
1982 for (i__ = 4; i__ <= i__3; i__ += 2) {
1983 idij += 2;
1984 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[
1985 i__ - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] *
1986 ch[i__ + (k + j * ch_dim2) * ch_dim1];
1987 c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__
1988 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ -
1989 1 + (k + j * ch_dim2) * ch_dim1];
1990/* L128: */
1991 }
1992/* L129: */
1993 }
1994/* L130: */
1995 }
1996 return;
1997} /* passf_ */
1998
1999/* Subroutine */ static void passf2(integer_t *ido, integer_t *l1, real_t *cc,
2000 real_t *ch, real_t *wa1)
2001{
2002 /* System generated locals */
2003 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2004
2005 /* Local variables */
2006 integer_t i__, k;
2007 real_t ti2, tr2;
2008
2009 /* Parameter adjustments */
2010 ch_dim1 = *ido;
2011 ch_dim2 = *l1;
2012 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2013 ch -= ch_offset;
2014 cc_dim1 = *ido;
2015 cc_offset = 1 + cc_dim1 * 3;
2016 cc -= cc_offset;
2017 --wa1;
2018
2019 /* Function Body */
2020 if (*ido > 2) {
2021 goto L102;
2022 }
2023 i__1 = *l1;
2024 for (k = 1; k <= i__1; ++k) {
2025 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] +
2026 cc[((k << 1) + 2) * cc_dim1 + 1];
2027 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1
2028 + 1] - cc[((k << 1) + 2) * cc_dim1 + 1];
2029 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] +
2030 cc[((k << 1) + 2) * cc_dim1 + 2];
2031 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1
2032 + 2] - cc[((k << 1) + 2) * cc_dim1 + 2];
2033/* L101: */
2034 }
2035 return;
2036L102:
2037 i__1 = *l1;
2038 for (k = 1; k <= i__1; ++k) {
2039 i__2 = *ido;
2040 for (i__ = 2; i__ <= i__2; i__ += 2) {
2041 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + ((k << 1) +
2042 1) * cc_dim1] + cc[i__ - 1 + ((k << 1) + 2) * cc_dim1];
2043 tr2 = cc[i__ - 1 + ((k << 1) + 1) * cc_dim1] - cc[i__ - 1 + ((k <<
2044 1) + 2) * cc_dim1];
2045 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) *
2046 cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1];
2047 ti2 = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2)
2048 * cc_dim1];
2049 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * ti2 -
2050 wa1[i__] * tr2;
2051 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * tr2
2052 + wa1[i__] * ti2;
2053/* L103: */
2054 }
2055/* L104: */
2056 }
2057 return;
2058} /* passf2_ */
2059
2060/* Subroutine */ static void passf3(integer_t *ido, integer_t *l1, real_t *cc,
2061 real_t *ch, real_t *wa1, real_t *wa2)
2062{
2063 /* Initialized data */
2064
2065 static real_t taur = REAL_CONSTANT(-0.5);
2066 static real_t taui =
2067 REAL_CONSTANT(-0.8660254037844386467637231707529361834740262690519031402790348975);
2068
2069 /* System generated locals */
2070 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2071
2072 /* Local variables */
2073 integer_t i__, k;
2074 real_t ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
2075
2076 /* Parameter adjustments */
2077 ch_dim1 = *ido;
2078 ch_dim2 = *l1;
2079 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2080 ch -= ch_offset;
2081 cc_dim1 = *ido;
2082 cc_offset = 1 + (cc_dim1 << 2);
2083 cc -= cc_offset;
2084 --wa1;
2085 --wa2;
2086
2087 /* Function Body */
2088 if (*ido != 2) {
2089 goto L102;
2090 }
2091 i__1 = *l1;
2092 for (k = 1; k <= i__1; ++k) {
2093 tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
2094 cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
2095 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
2096 ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
2097 ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
2098 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
2099 cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] - cc[(k * 3 + 3) *
2100 cc_dim1 + 1]);
2101 ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] - cc[(k * 3 + 3) *
2102 cc_dim1 + 2]);
2103 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
2104 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
2105 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
2106 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
2107/* L101: */
2108 }
2109 return;
2110L102:
2111 i__1 = *l1;
2112 for (k = 1; k <= i__1; ++k) {
2113 i__2 = *ido;
2114 for (i__ = 2; i__ <= i__2; i__ += 2) {
2115 tr2 = cc[i__ - 1 + (k * 3 + 2) * cc_dim1] + cc[i__ - 1 + (k * 3 +
2116 3) * cc_dim1];
2117 cr2 = cc[i__ - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
2118 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + (k * 3 + 1) *
2119 cc_dim1] + tr2;
2120 ti2 = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) *
2121 cc_dim1];
2122 ci2 = cc[i__ + (k * 3 + 1) * cc_dim1] + taur * ti2;
2123 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) *
2124 cc_dim1] + ti2;
2125 cr3 = taui * (cc[i__ - 1 + (k * 3 + 2) * cc_dim1] - cc[i__ - 1 + (
2126 k * 3 + 3) * cc_dim1]);
2127 ci3 = taui * (cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 +
2128 3) * cc_dim1]);
2129 dr2 = cr2 - ci3;
2130 dr3 = cr2 + ci3;
2131 di2 = ci2 + cr3;
2132 di3 = ci2 - cr3;
2133 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * di2 -
2134 wa1[i__] * dr2;
2135 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * dr2
2136 + wa1[i__] * di2;
2137 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * di3 - wa2[
2138 i__] * dr3;
2139 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * dr3 +
2140 wa2[i__] * di3;
2141/* L103: */
2142 }
2143/* L104: */
2144 }
2145 return;
2146} /* passf3_ */
2147
2148/* Subroutine */ static void passf4(integer_t *ido, integer_t *l1, real_t *cc,
2149 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3)
2150{
2151 /* System generated locals */
2152 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2153
2154 /* Local variables */
2155 integer_t i__, k;
2156 real_t ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
2157 tr3, tr4;
2158
2159 /* Parameter adjustments */
2160 ch_dim1 = *ido;
2161 ch_dim2 = *l1;
2162 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2163 ch -= ch_offset;
2164 cc_dim1 = *ido;
2165 cc_offset = 1 + cc_dim1 * 5;
2166 cc -= cc_offset;
2167 --wa1;
2168 --wa2;
2169 --wa3;
2170
2171 /* Function Body */
2172 if (*ido != 2) {
2173 goto L102;
2174 }
2175 i__1 = *l1;
2176 for (k = 1; k <= i__1; ++k) {
2177 ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] - cc[((k << 2) + 3) * cc_dim1
2178 + 2];
2179 ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] + cc[((k << 2) + 3) * cc_dim1
2180 + 2];
2181 tr4 = cc[((k << 2) + 2) * cc_dim1 + 2] - cc[((k << 2) + 4) * cc_dim1
2182 + 2];
2183 ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] + cc[((k << 2) + 4) * cc_dim1
2184 + 2];
2185 tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - cc[((k << 2) + 3) * cc_dim1
2186 + 1];
2187 tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + cc[((k << 2) + 3) * cc_dim1
2188 + 1];
2189 ti4 = cc[((k << 2) + 4) * cc_dim1 + 1] - cc[((k << 2) + 2) * cc_dim1
2190 + 1];
2191 tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] + cc[((k << 2) + 4) * cc_dim1
2192 + 1];
2193 ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
2194 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
2195 ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
2196 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
2197 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
2198 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
2199 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
2200 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
2201/* L101: */
2202 }
2203 return;
2204L102:
2205 i__1 = *l1;
2206 for (k = 1; k <= i__1; ++k) {
2207 i__2 = *ido;
2208 for (i__ = 2; i__ <= i__2; i__ += 2) {
2209 ti1 = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3)
2210 * cc_dim1];
2211 ti2 = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3)
2212 * cc_dim1];
2213 ti3 = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4)
2214 * cc_dim1];
2215 tr4 = cc[i__ + ((k << 2) + 2) * cc_dim1] - cc[i__ + ((k << 2) + 4)
2216 * cc_dim1];
2217 tr1 = cc[i__ - 1 + ((k << 2) + 1) * cc_dim1] - cc[i__ - 1 + ((k <<
2218 2) + 3) * cc_dim1];
2219 tr2 = cc[i__ - 1 + ((k << 2) + 1) * cc_dim1] + cc[i__ - 1 + ((k <<
2220 2) + 3) * cc_dim1];
2221 ti4 = cc[i__ - 1 + ((k << 2) + 4) * cc_dim1] - cc[i__ - 1 + ((k <<
2222 2) + 2) * cc_dim1];
2223 tr3 = cc[i__ - 1 + ((k << 2) + 2) * cc_dim1] + cc[i__ - 1 + ((k <<
2224 2) + 4) * cc_dim1];
2225 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
2226 cr3 = tr2 - tr3;
2227 ch[i__ + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
2228 ci3 = ti2 - ti3;
2229 cr2 = tr1 + tr4;
2230 cr4 = tr1 - tr4;
2231 ci2 = ti1 + ti4;
2232 ci4 = ti1 - ti4;
2233 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * cr2
2234 + wa1[i__] * ci2;
2235 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * ci2 -
2236 wa1[i__] * cr2;
2237 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * cr3 +
2238 wa2[i__] * ci3;
2239 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * ci3 - wa2[
2240 i__] * cr3;
2241 ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * cr4
2242 + wa3[i__] * ci4;
2243 ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * ci4 -
2244 wa3[i__] * cr4;
2245/* L103: */
2246 }
2247/* L104: */
2248 }
2249 return;
2250} /* passf4_ */
2251
2252/* Subroutine */ static void passf5(integer_t *ido, integer_t *l1, real_t *cc,
2253 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3,
2254 real_t *wa4)
2255{
2256 /* Initialized data */
2257
2258 static real_t tr11 =
2259 REAL_CONSTANT(0.3090169943749474241022934171828195886015458990288143106772431137);
2260 static real_t ti11 =
2261 REAL_CONSTANT(-0.9510565162951535721164393337938214340569863412575022244730564442);
2262 static real_t tr12 =
2263 REAL_CONSTANT(-0.8090169943749474241022934171828190588601545899028814310677431135);
2264 static real_t ti12 =
2265 REAL_CONSTANT(-0.5877852522924731291687059546390727685976524376431459107227248076);
2266
2267 /* System generated locals */
2268 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2269
2270 /* Local variables */
2271 integer_t i__, k;
2272 real_t ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
2273 ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
2274
2275 /* Parameter adjustments */
2276 ch_dim1 = *ido;
2277 ch_dim2 = *l1;
2278 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2279 ch -= ch_offset;
2280 cc_dim1 = *ido;
2281 cc_offset = 1 + cc_dim1 * 6;
2282 cc -= cc_offset;
2283 --wa1;
2284 --wa2;
2285 --wa3;
2286 --wa4;
2287
2288 /* Function Body */
2289 if (*ido != 2) {
2290 goto L102;
2291 }
2292 i__1 = *l1;
2293 for (k = 1; k <= i__1; ++k) {
2294 ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
2295 ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
2296 ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
2297 ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
2298 tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
2299 tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
2300 tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
2301 tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
2302 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + tr2
2303 + tr3;
2304 ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] + ti2
2305 + ti3;
2306 cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
2307 ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
2308 cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
2309 ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
2310 cr5 = ti11 * tr5 + ti12 * tr4;
2311 ci5 = ti11 * ti5 + ti12 * ti4;
2312 cr4 = ti12 * tr5 - ti11 * tr4;
2313 ci4 = ti12 * ti5 - ti11 * ti4;
2314 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
2315 ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
2316 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
2317 ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
2318 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
2319 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
2320 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
2321 ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
2322/* L101: */
2323 }
2324 return;
2325L102:
2326 i__1 = *l1;
2327 for (k = 1; k <= i__1; ++k) {
2328 i__2 = *ido;
2329 for (i__ = 2; i__ <= i__2; i__ += 2) {
2330 ti5 = cc[i__ + (k * 5 + 2) * cc_dim1] - cc[i__ + (k * 5 + 5) *
2331 cc_dim1];
2332 ti2 = cc[i__ + (k * 5 + 2) * cc_dim1] + cc[i__ + (k * 5 + 5) *
2333 cc_dim1];
2334 ti4 = cc[i__ + (k * 5 + 3) * cc_dim1] - cc[i__ + (k * 5 + 4) *
2335 cc_dim1];
2336 ti3 = cc[i__ + (k * 5 + 3) * cc_dim1] + cc[i__ + (k * 5 + 4) *
2337 cc_dim1];
2338 tr5 = cc[i__ - 1 + (k * 5 + 2) * cc_dim1] - cc[i__ - 1 + (k * 5 +
2339 5) * cc_dim1];
2340 tr2 = cc[i__ - 1 + (k * 5 + 2) * cc_dim1] + cc[i__ - 1 + (k * 5 +
2341 5) * cc_dim1];
2342 tr4 = cc[i__ - 1 + (k * 5 + 3) * cc_dim1] - cc[i__ - 1 + (k * 5 +
2343 4) * cc_dim1];
2344 tr3 = cc[i__ - 1 + (k * 5 + 3) * cc_dim1] + cc[i__ - 1 + (k * 5 +
2345 4) * cc_dim1];
2346 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + (k * 5 + 1) *
2347 cc_dim1] + tr2 + tr3;
2348 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 5 + 1) *
2349 cc_dim1] + ti2 + ti3;
2350 cr2 = cc[i__ - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 *
2351 tr3;
2352 ci2 = cc[i__ + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
2353 cr3 = cc[i__ - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 *
2354 tr3;
2355 ci3 = cc[i__ + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
2356 cr5 = ti11 * tr5 + ti12 * tr4;
2357 ci5 = ti11 * ti5 + ti12 * ti4;
2358 cr4 = ti12 * tr5 - ti11 * tr4;
2359 ci4 = ti12 * ti5 - ti11 * ti4;
2360 dr3 = cr3 - ci4;
2361 dr4 = cr3 + ci4;
2362 di3 = ci3 + cr4;
2363 di4 = ci3 - cr4;
2364 dr5 = cr2 + ci5;
2365 dr2 = cr2 - ci5;
2366 di5 = ci2 - cr5;
2367 di2 = ci2 + cr5;
2368 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * dr2
2369 + wa1[i__] * di2;
2370 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 1] * di2 -
2371 wa1[i__] * dr2;
2372 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * dr3 +
2373 wa2[i__] * di3;
2374 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 1] * di3 - wa2[
2375 i__] * dr3;
2376 ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * dr4
2377 + wa3[i__] * di4;
2378 ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 1] * di4 -
2379 wa3[i__] * dr4;
2380 ch[i__ - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i__ - 1] * dr5 +
2381 wa4[i__] * di5;
2382 ch[i__ + (k + ch_dim2 * 5) * ch_dim1] = wa4[i__ - 1] * di5 - wa4[
2383 i__] * dr5;
2384/* L103: */
2385 }
2386/* L104: */
2387 }
2388 return;
2389} /* passf5_ */
2390
2391/* Subroutine */ static void radb2(integer_t *ido, integer_t *l1, real_t *cc,
2392 real_t *ch, real_t *wa1)
2393{
2394 /* System generated locals */
2395 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2396
2397 /* Local variables */
2398 integer_t i__, k, ic;
2399 real_t ti2, tr2;
2400 integer_t idp2;
2401
2402 /* Parameter adjustments */
2403 ch_dim1 = *ido;
2404 ch_dim2 = *l1;
2405 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2406 ch -= ch_offset;
2407 cc_dim1 = *ido;
2408 cc_offset = 1 + cc_dim1 * 3;
2409 cc -= cc_offset;
2410 --wa1;
2411
2412 /* Function Body */
2413 i__1 = *l1;
2414 for (k = 1; k <= i__1; ++k) {
2415 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] +
2416 cc[*ido + ((k << 1) + 2) * cc_dim1];
2417 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1
2418 + 1] - cc[*ido + ((k << 1) + 2) * cc_dim1];
2419/* L101: */
2420 }
2421 if ((i__1 = *ido - 2) < 0) {
2422 goto L107;
2423 } else if (i__1 == 0) {
2424 goto L105;
2425 } else {
2426 goto L102;
2427 }
2428L102:
2429 idp2 = *ido + 2;
2430 i__1 = *l1;
2431 for (k = 1; k <= i__1; ++k) {
2432 i__2 = *ido;
2433 for (i__ = 3; i__ <= i__2; i__ += 2) {
2434 ic = idp2 - i__;
2435 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + ((k << 1) +
2436 1) * cc_dim1] + cc[ic - 1 + ((k << 1) + 2) * cc_dim1];
2437 tr2 = cc[i__ - 1 + ((k << 1) + 1) * cc_dim1] - cc[ic - 1 + ((k <<
2438 1) + 2) * cc_dim1];
2439 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) *
2440 cc_dim1] - cc[ic + ((k << 1) + 2) * cc_dim1];
2441 ti2 = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[ic + ((k << 1) + 2)
2442 * cc_dim1];
2443 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * tr2
2444 - wa1[i__ - 1] * ti2;
2445 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * ti2 +
2446 wa1[i__ - 1] * tr2;
2447/* L103: */
2448 }
2449/* L104: */
2450 }
2451 if (*ido % 2 == 1) {
2452 return;
2453 }
2454L105:
2455 i__1 = *l1;
2456 for (k = 1; k <= i__1; ++k) {
2457 ch[*ido + (k + ch_dim2) * ch_dim1] = cc[*ido + ((k << 1) + 1) *
2458 cc_dim1] + cc[*ido + ((k << 1) + 1) * cc_dim1];
2459 ch[*ido + (k + (ch_dim2 << 1)) * ch_dim1] = -(cc[((k << 1) + 2) *
2460 cc_dim1 + 1] + cc[((k << 1) + 2) * cc_dim1 + 1]);
2461/* L106: */
2462 }
2463L107:
2464 return;
2465} /* radb2_ */
2466
2467/* Subroutine */ static void radb3(integer_t *ido, integer_t *l1, real_t *cc,
2468 real_t *ch, real_t *wa1, real_t *wa2)
2469{
2470 /* Initialized data */
2471
2472 static real_t taur = REAL_CONSTANT(-0.5);
2473 static real_t taui =
2474 REAL_CONSTANT(0.8660254037844386467637231707529361834710262690519031402790348975);
2475
2476 /* System generated locals */
2477 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2478
2479 /* Local variables */
2480 integer_t i__, k, ic;
2481 real_t ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
2482 integer_t idp2;
2483
2484 /* Parameter adjustments */
2485 ch_dim1 = *ido;
2486 ch_dim2 = *l1;
2487 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2488 ch -= ch_offset;
2489 cc_dim1 = *ido;
2490 cc_offset = 1 + (cc_dim1 << 2);
2491 cc -= cc_offset;
2492 --wa1;
2493 --wa2;
2494
2495 /* Function Body */
2496 i__1 = *l1;
2497 for (k = 1; k <= i__1; ++k) {
2498 tr2 = cc[*ido + (k * 3 + 2) * cc_dim1] + cc[*ido + (k * 3 + 2) *
2499 cc_dim1];
2500 cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
2501 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
2502 ci3 = taui * (cc[(k * 3 + 3) * cc_dim1 + 1] + cc[(k * 3 + 3) *
2503 cc_dim1 + 1]);
2504 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
2505 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
2506/* L101: */
2507 }
2508 if (*ido == 1) {
2509 return;
2510 }
2511 idp2 = *ido + 2;
2512 i__1 = *l1;
2513 for (k = 1; k <= i__1; ++k) {
2514 i__2 = *ido;
2515 for (i__ = 3; i__ <= i__2; i__ += 2) {
2516 ic = idp2 - i__;
2517 tr2 = cc[i__ - 1 + (k * 3 + 3) * cc_dim1] + cc[ic - 1 + (k * 3 +
2518 2) * cc_dim1];
2519 cr2 = cc[i__ - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
2520 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + (k * 3 + 1) *
2521 cc_dim1] + tr2;
2522 ti2 = cc[i__ + (k * 3 + 3) * cc_dim1] - cc[ic + (k * 3 + 2) *
2523 cc_dim1];
2524 ci2 = cc[i__ + (k * 3 + 1) * cc_dim1] + taur * ti2;
2525 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) *
2526 cc_dim1] + ti2;
2527 cr3 = taui * (cc[i__ - 1 + (k * 3 + 3) * cc_dim1] - cc[ic - 1 + (
2528 k * 3 + 2) * cc_dim1]);
2529 ci3 = taui * (cc[i__ + (k * 3 + 3) * cc_dim1] + cc[ic + (k * 3 +
2530 2) * cc_dim1]);
2531 dr2 = cr2 - ci3;
2532 dr3 = cr2 + ci3;
2533 di2 = ci2 + cr3;
2534 di3 = ci2 - cr3;
2535 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * dr2
2536 - wa1[i__ - 1] * di2;
2537 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * di2 +
2538 wa1[i__ - 1] * dr2;
2539 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 2] * dr3 -
2540 wa2[i__ - 1] * di3;
2541 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 2] * di3 + wa2[
2542 i__ - 1] * dr3;
2543/* L102: */
2544 }
2545/* L103: */
2546 }
2547 return;
2548} /* radb3_ */
2549
2550/* Subroutine */ static void radb4(integer_t *ido, integer_t *l1, real_t *cc,
2551 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3)
2552{
2553 /* Initialized data */
2554
2555 static real_t sqrt2 =
2556 REAL_CONSTANT(1.41421356237309504880168872420969807856967187536948073176679738);
2557
2558 /* System generated locals */
2559 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2560
2561 /* Local variables */
2562 integer_t i__, k, ic;
2563 real_t ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
2564 tr3, tr4;
2565 integer_t idp2;
2566
2567 /* Parameter adjustments */
2568 ch_dim1 = *ido;
2569 ch_dim2 = *l1;
2570 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2571 ch -= ch_offset;
2572 cc_dim1 = *ido;
2573 cc_offset = 1 + cc_dim1 * 5;
2574 cc -= cc_offset;
2575 --wa1;
2576 --wa2;
2577 --wa3;
2578
2579 /* Function Body */
2580 i__1 = *l1;
2581 for (k = 1; k <= i__1; ++k) {
2582 tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - cc[*ido + ((k << 2) + 4) *
2583 cc_dim1];
2584 tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + cc[*ido + ((k << 2) + 4) *
2585 cc_dim1];
2586 tr3 = cc[*ido + ((k << 2) + 2) * cc_dim1] + cc[*ido + ((k << 2) + 2) *
2587 cc_dim1];
2588 tr4 = cc[((k << 2) + 3) * cc_dim1 + 1] + cc[((k << 2) + 3) * cc_dim1
2589 + 1];
2590 ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
2591 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 - tr4;
2592 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
2593 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 + tr4;
2594/* L101: */
2595 }
2596 if ((i__1 = *ido - 2) < 0) {
2597 goto L107;
2598 } else if (i__1 == 0) {
2599 goto L105;
2600 } else {
2601 goto L102;
2602 }
2603L102:
2604 idp2 = *ido + 2;
2605 i__1 = *l1;
2606 for (k = 1; k <= i__1; ++k) {
2607 i__2 = *ido;
2608 for (i__ = 3; i__ <= i__2; i__ += 2) {
2609 ic = idp2 - i__;
2610 ti1 = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[ic + ((k << 2) + 4)
2611 * cc_dim1];
2612 ti2 = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[ic + ((k << 2) + 4)
2613 * cc_dim1];
2614 ti3 = cc[i__ + ((k << 2) + 3) * cc_dim1] - cc[ic + ((k << 2) + 2)
2615 * cc_dim1];
2616 tr4 = cc[i__ + ((k << 2) + 3) * cc_dim1] + cc[ic + ((k << 2) + 2)
2617 * cc_dim1];
2618 tr1 = cc[i__ - 1 + ((k << 2) + 1) * cc_dim1] - cc[ic - 1 + ((k <<
2619 2) + 4) * cc_dim1];
2620 tr2 = cc[i__ - 1 + ((k << 2) + 1) * cc_dim1] + cc[ic - 1 + ((k <<
2621 2) + 4) * cc_dim1];
2622 ti4 = cc[i__ - 1 + ((k << 2) + 3) * cc_dim1] - cc[ic - 1 + ((k <<
2623 2) + 2) * cc_dim1];
2624 tr3 = cc[i__ - 1 + ((k << 2) + 3) * cc_dim1] + cc[ic - 1 + ((k <<
2625 2) + 2) * cc_dim1];
2626 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
2627 cr3 = tr2 - tr3;
2628 ch[i__ + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
2629 ci3 = ti2 - ti3;
2630 cr2 = tr1 - tr4;
2631 cr4 = tr1 + tr4;
2632 ci2 = ti1 + ti4;
2633 ci4 = ti1 - ti4;
2634 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * cr2
2635 - wa1[i__ - 1] * ci2;
2636 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * ci2 +
2637 wa1[i__ - 1] * cr2;
2638 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 2] * cr3 -
2639 wa2[i__ - 1] * ci3;
2640 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 2] * ci3 + wa2[
2641 i__ - 1] * cr3;
2642 ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 2] * cr4
2643 - wa3[i__ - 1] * ci4;
2644 ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 2] * ci4 +
2645 wa3[i__ - 1] * cr4;
2646/* L103: */
2647 }
2648/* L104: */
2649 }
2650 if (*ido % 2 == 1) {
2651 return;
2652 }
2653L105:
2654 i__1 = *l1;
2655 for (k = 1; k <= i__1; ++k) {
2656 ti1 = cc[((k << 2) + 2) * cc_dim1 + 1] + cc[((k << 2) + 4) * cc_dim1
2657 + 1];
2658 ti2 = cc[((k << 2) + 4) * cc_dim1 + 1] - cc[((k << 2) + 2) * cc_dim1
2659 + 1];
2660 tr1 = cc[*ido + ((k << 2) + 1) * cc_dim1] - cc[*ido + ((k << 2) + 3) *
2661 cc_dim1];
2662 tr2 = cc[*ido + ((k << 2) + 1) * cc_dim1] + cc[*ido + ((k << 2) + 3) *
2663 cc_dim1];
2664 ch[*ido + (k + ch_dim2) * ch_dim1] = tr2 + tr2;
2665 ch[*ido + (k + (ch_dim2 << 1)) * ch_dim1] = sqrt2 * (tr1 - ti1);
2666 ch[*ido + (k + ch_dim2 * 3) * ch_dim1] = ti2 + ti2;
2667 ch[*ido + (k + (ch_dim2 << 2)) * ch_dim1] = -sqrt2 * (tr1 + ti1);
2668/* L106: */
2669 }
2670L107:
2671 return;
2672} /* radb4_ */
2673
2674/* Subroutine */ static void radb5(integer_t *ido, integer_t *l1, real_t *cc,
2675 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3,
2676 real_t *wa4)
2677{
2678 /* Initialized data */
2679
2680 static real_t tr11 =
2681 REAL_CONSTANT(0.3090169943749474241022934171828195886015458990288143106772431137);
2682 static real_t ti11 =
2683 REAL_CONSTANT(0.9510565162951535721164393337938214340569863412575022244730564442);
2684 static real_t tr12 =
2685 REAL_CONSTANT(-0.8090169943749474241022934171828190588601545899028814310677431135);
2686 static real_t ti12 =
2687 REAL_CONSTANT(0.5877852522924731291687059546390727685976524376431459107227248076);
2688
2689 /* System generated locals */
2690 integer_t cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;
2691
2692 /* Local variables */
2693 integer_t i__, k, ic;
2694 real_t ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
2695 ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
2696 integer_t idp2;
2697
2698 /* Parameter adjustments */
2699 ch_dim1 = *ido;
2700 ch_dim2 = *l1;
2701 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2702 ch -= ch_offset;
2703 cc_dim1 = *ido;
2704 cc_offset = 1 + cc_dim1 * 6;
2705 cc -= cc_offset;
2706 --wa1;
2707 --wa2;
2708 --wa3;
2709 --wa4;
2710
2711 /* Function Body */
2712 i__1 = *l1;
2713 for (k = 1; k <= i__1; ++k) {
2714 ti5 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 3) * cc_dim1 + 1];
2715 ti4 = cc[(k * 5 + 5) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
2716 tr2 = cc[*ido + (k * 5 + 2) * cc_dim1] + cc[*ido + (k * 5 + 2) *
2717 cc_dim1];
2718 tr3 = cc[*ido + (k * 5 + 4) * cc_dim1] + cc[*ido + (k * 5 + 4) *
2719 cc_dim1];
2720 ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + tr2
2721 + tr3;
2722 cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
2723 cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
2724 ci5 = ti11 * ti5 + ti12 * ti4;
2725 ci4 = ti12 * ti5 - ti11 * ti4;
2726 ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
2727 ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
2728 ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
2729 ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
2730/* L101: */
2731 }
2732 if (*ido == 1) {
2733 return;
2734 }
2735 idp2 = *ido + 2;
2736 i__1 = *l1;
2737 for (k = 1; k <= i__1; ++k) {
2738 i__2 = *ido;
2739 for (i__ = 3; i__ <= i__2; i__ += 2) {
2740 ic = idp2 - i__;
2741 ti5 = cc[i__ + (k * 5 + 3) * cc_dim1] + cc[ic + (k * 5 + 2) *
2742 cc_dim1];
2743 ti2 = cc[i__ + (k * 5 + 3) * cc_dim1] - cc[ic + (k * 5 + 2) *
2744 cc_dim1];
2745 ti4 = cc[i__ + (k * 5 + 5) * cc_dim1] + cc[ic + (k * 5 + 4) *
2746 cc_dim1];
2747 ti3 = cc[i__ + (k * 5 + 5) * cc_dim1] - cc[ic + (k * 5 + 4) *
2748 cc_dim1];
2749 tr5 = cc[i__ - 1 + (k * 5 + 3) * cc_dim1] - cc[ic - 1 + (k * 5 +
2750 2) * cc_dim1];
2751 tr2 = cc[i__ - 1 + (k * 5 + 3) * cc_dim1] + cc[ic - 1 + (k * 5 +
2752 2) * cc_dim1];
2753 tr4 = cc[i__ - 1 + (k * 5 + 5) * cc_dim1] - cc[ic - 1 + (k * 5 +
2754 4) * cc_dim1];
2755 tr3 = cc[i__ - 1 + (k * 5 + 5) * cc_dim1] + cc[ic - 1 + (k * 5 +
2756 4) * cc_dim1];
2757 ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + (k * 5 + 1) *
2758 cc_dim1] + tr2 + tr3;
2759 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 5 + 1) *
2760 cc_dim1] + ti2 + ti3;
2761 cr2 = cc[i__ - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 *
2762 tr3;
2763 ci2 = cc[i__ + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
2764 cr3 = cc[i__ - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 *
2765 tr3;
2766 ci3 = cc[i__ + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
2767 cr5 = ti11 * tr5 + ti12 * tr4;
2768 ci5 = ti11 * ti5 + ti12 * ti4;
2769 cr4 = ti12 * tr5 - ti11 * tr4;
2770 ci4 = ti12 * ti5 - ti11 * ti4;
2771 dr3 = cr3 - ci4;
2772 dr4 = cr3 + ci4;
2773 di3 = ci3 + cr4;
2774 di4 = ci3 - cr4;
2775 dr5 = cr2 + ci5;
2776 dr2 = cr2 - ci5;
2777 di5 = ci2 - cr5;
2778 di2 = ci2 + cr5;
2779 ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * dr2
2780 - wa1[i__ - 1] * di2;
2781 ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * di2 +
2782 wa1[i__ - 1] * dr2;
2783 ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 2] * dr3 -
2784 wa2[i__ - 1] * di3;
2785 ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = wa2[i__ - 2] * di3 + wa2[
2786 i__ - 1] * dr3;
2787 ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 2] * dr4
2788 - wa3[i__ - 1] * di4;
2789 ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i__ - 2] * di4 +
2790 wa3[i__ - 1] * dr4;
2791 ch[i__ - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i__ - 2] * dr5 -
2792 wa4[i__ - 1] * di5;
2793 ch[i__ + (k + ch_dim2 * 5) * ch_dim1] = wa4[i__ - 2] * di5 + wa4[
2794 i__ - 1] * dr5;
2795/* L102: */
2796 }
2797/* L103: */
2798 }
2799 return;
2800} /* radb5_ */
2801
2802/* Subroutine */ static void radbg(integer_t *ido, integer_t *ip, integer_t *l1, integer_t *
2803 idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch,
2804 real_t *ch2, real_t *wa)
2805{
2806 /* Initialized data */
2807
2808 static real_t tpi =
2809 REAL_CONSTANT(6.283185307179586476925286766559005768394338798750116419498891846);
2810
2811 /* System generated locals */
2812 integer_t ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
2813 c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset,
2814 i__1, i__2, i__3;
2815
2816 /* Local variables */
2817 integer_t i__, j, k, l, j2, ic, jc, lc, ik, is;
2818 real_t dc2, ai1, ai2, ar1, ar2, ds2;
2819 integer_t nbd;
2820 real_t dcp, arg, dsp, ar1h, ar2h;
2821 integer_t idp2, ipp2, idij, ipph;
2822
2823 /* Parameter adjustments */
2824 ch_dim1 = *ido;
2825 ch_dim2 = *l1;
2826 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
2827 ch -= ch_offset;
2828 c1_dim1 = *ido;
2829 c1_dim2 = *l1;
2830 c1_offset = 1 + c1_dim1 * (1 + c1_dim2);
2831 c1 -= c1_offset;
2832 cc_dim1 = *ido;
2833 cc_dim2 = *ip;
2834 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
2835 cc -= cc_offset;
2836 ch2_dim1 = *idl1;
2837 ch2_offset = 1 + ch2_dim1;
2838 ch2 -= ch2_offset;
2839 c2_dim1 = *idl1;
2840 c2_offset = 1 + c2_dim1;
2841 c2 -= c2_offset;
2842 --wa;
2843
2844 /* Function Body */
2845 arg = tpi / (real_t) (*ip);
2846 dcp = cos(arg);
2847 dsp = sin(arg);
2848 idp2 = *ido + 2;
2849 nbd = (*ido - 1) / 2;
2850 ipp2 = *ip + 2;
2851 ipph = (*ip + 1) / 2;
2852 if (*ido < *l1) {
2853 goto L103;
2854 }
2855 i__1 = *l1;
2856 for (k = 1; k <= i__1; ++k) {
2857 i__2 = *ido;
2858 for (i__ = 1; i__ <= i__2; ++i__) {
2859 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) *
2860 cc_dim1];
2861/* L101: */
2862 }
2863/* L102: */
2864 }
2865 goto L106;
2866L103:
2867 i__1 = *ido;
2868 for (i__ = 1; i__ <= i__1; ++i__) {
2869 i__2 = *l1;
2870 for (k = 1; k <= i__2; ++k) {
2871 ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) *
2872 cc_dim1];
2873/* L104: */
2874 }
2875/* L105: */
2876 }
2877L106:
2878 i__1 = ipph;
2879 for (j = 2; j <= i__1; ++j) {
2880 jc = ipp2 - j;
2881 j2 = j + j;
2882 i__2 = *l1;
2883 for (k = 1; k <= i__2; ++k) {
2884 ch[(k + j * ch_dim2) * ch_dim1 + 1] = cc[*ido + (j2 - 2 + k *
2885 cc_dim2) * cc_dim1] + cc[*ido + (j2 - 2 + k * cc_dim2) *
2886 cc_dim1];
2887 ch[(k + jc * ch_dim2) * ch_dim1 + 1] = cc[(j2 - 1 + k * cc_dim2) *
2888 cc_dim1 + 1] + cc[(j2 - 1 + k * cc_dim2) * cc_dim1 + 1];
2889/* L107: */
2890 }
2891/* L108: */
2892 }
2893 if (*ido == 1) {
2894 goto L116;
2895 }
2896 if (nbd < *l1) {
2897 goto L112;
2898 }
2899 i__1 = ipph;
2900 for (j = 2; j <= i__1; ++j) {
2901 jc = ipp2 - j;
2902 i__2 = *l1;
2903 for (k = 1; k <= i__2; ++k) {
2904 i__3 = *ido;
2905 for (i__ = 3; i__ <= i__3; i__ += 2) {
2906 ic = idp2 - i__;
2907 ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1] = cc[i__ - 1 + ((j
2908 << 1) - 1 + k * cc_dim2) * cc_dim1] + cc[ic - 1 + ((j
2909 << 1) - 2 + k * cc_dim2) * cc_dim1];
2910 ch[i__ - 1 + (k + jc * ch_dim2) * ch_dim1] = cc[i__ - 1 + ((j
2911 << 1) - 1 + k * cc_dim2) * cc_dim1] - cc[ic - 1 + ((j
2912 << 1) - 2 + k * cc_dim2) * cc_dim1];
2913 ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + ((j << 1) -
2914 1 + k * cc_dim2) * cc_dim1] - cc[ic + ((j << 1) - 2 +
2915 k * cc_dim2) * cc_dim1];
2916 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + ((j << 1) -
2917 1 + k * cc_dim2) * cc_dim1] + cc[ic + ((j << 1) - 2 +
2918 k * cc_dim2) * cc_dim1];
2919/* L109: */
2920 }
2921/* L110: */
2922 }
2923/* L111: */
2924 }
2925 goto L116;
2926L112:
2927 i__1 = ipph;
2928 for (j = 2; j <= i__1; ++j) {
2929 jc = ipp2 - j;
2930 i__2 = *ido;
2931 for (i__ = 3; i__ <= i__2; i__ += 2) {
2932 ic = idp2 - i__;
2933 i__3 = *l1;
2934 for (k = 1; k <= i__3; ++k) {
2935 ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1] = cc[i__ - 1 + ((j
2936 << 1) - 1 + k * cc_dim2) * cc_dim1] + cc[ic - 1 + ((j
2937 << 1) - 2 + k * cc_dim2) * cc_dim1];
2938 ch[i__ - 1 + (k + jc * ch_dim2) * ch_dim1] = cc[i__ - 1 + ((j
2939 << 1) - 1 + k * cc_dim2) * cc_dim1] - cc[ic - 1 + ((j
2940 << 1) - 2 + k * cc_dim2) * cc_dim1];
2941 ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + ((j << 1) -
2942 1 + k * cc_dim2) * cc_dim1] - cc[ic + ((j << 1) - 2 +
2943 k * cc_dim2) * cc_dim1];
2944 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + ((j << 1) -
2945 1 + k * cc_dim2) * cc_dim1] + cc[ic + ((j << 1) - 2 +
2946 k * cc_dim2) * cc_dim1];
2947/* L113: */
2948 }
2949/* L114: */
2950 }
2951/* L115: */
2952 }
2953L116:
2954 ar1 = REAL_CONSTANT(1.0);
2955 ai1 = REAL_CONSTANT(0.0);
2956 i__1 = ipph;
2957 for (l = 2; l <= i__1; ++l) {
2958 lc = ipp2 - l;
2959 ar1h = dcp * ar1 - dsp * ai1;
2960 ai1 = dcp * ai1 + dsp * ar1;
2961 ar1 = ar1h;
2962 i__2 = *idl1;
2963 for (ik = 1; ik <= i__2; ++ik) {
2964 c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + ar1 * ch2[ik + (
2965 ch2_dim1 << 1)];
2966 c2[ik + lc * c2_dim1] = ai1 * ch2[ik + *ip * ch2_dim1];
2967/* L117: */
2968 }
2969 dc2 = ar1;
2970 ds2 = ai1;
2971 ar2 = ar1;
2972 ai2 = ai1;
2973 i__2 = ipph;
2974 for (j = 3; j <= i__2; ++j) {
2975 jc = ipp2 - j;
2976 ar2h = dc2 * ar2 - ds2 * ai2;
2977 ai2 = dc2 * ai2 + ds2 * ar2;
2978 ar2 = ar2h;
2979 i__3 = *idl1;
2980 for (ik = 1; ik <= i__3; ++ik) {
2981 c2[ik + l * c2_dim1] += ar2 * ch2[ik + j * ch2_dim1];
2982 c2[ik + lc * c2_dim1] += ai2 * ch2[ik + jc * ch2_dim1];
2983/* L118: */
2984 }
2985/* L119: */
2986 }
2987/* L120: */
2988 }
2989 i__1 = ipph;
2990 for (j = 2; j <= i__1; ++j) {
2991 i__2 = *idl1;
2992 for (ik = 1; ik <= i__2; ++ik) {
2993 ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
2994/* L121: */
2995 }
2996/* L122: */
2997 }
2998 i__1 = ipph;
2999 for (j = 2; j <= i__1; ++j) {
3000 jc = ipp2 - j;
3001 i__2 = *l1;
3002 for (k = 1; k <= i__2; ++k) {
3003 ch[(k + j * ch_dim2) * ch_dim1 + 1] = c1[(k + j * c1_dim2) *
3004 c1_dim1 + 1] - c1[(k + jc * c1_dim2) * c1_dim1 + 1];
3005 ch[(k + jc * ch_dim2) * ch_dim1 + 1] = c1[(k + j * c1_dim2) *
3006 c1_dim1 + 1] + c1[(k + jc * c1_dim2) * c1_dim1 + 1];
3007/* L123: */
3008 }
3009/* L124: */
3010 }
3011 if (*ido == 1) {
3012 goto L132;
3013 }
3014 if (nbd < *l1) {
3015 goto L128;
3016 }
3017 i__1 = ipph;
3018 for (j = 2; j <= i__1; ++j) {
3019 jc = ipp2 - j;
3020 i__2 = *l1;
3021 for (k = 1; k <= i__2; ++k) {
3022 i__3 = *ido;
3023 for (i__ = 3; i__ <= i__3; i__ += 2) {
3024 ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1] = c1[i__ - 1 + (k +
3025 j * c1_dim2) * c1_dim1] - c1[i__ + (k + jc * c1_dim2)
3026 * c1_dim1];
3027 ch[i__ - 1 + (k + jc * ch_dim2) * ch_dim1] = c1[i__ - 1 + (k
3028 + j * c1_dim2) * c1_dim1] + c1[i__ + (k + jc *
3029 c1_dim2) * c1_dim1];
3030 ch[i__ + (k + j * ch_dim2) * ch_dim1] = c1[i__ + (k + j *
3031 c1_dim2) * c1_dim1] + c1[i__ - 1 + (k + jc * c1_dim2)
3032 * c1_dim1];
3033 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = c1[i__ + (k + j *
3034 c1_dim2) * c1_dim1] - c1[i__ - 1 + (k + jc * c1_dim2)
3035 * c1_dim1];
3036/* L125: */
3037 }
3038/* L126: */
3039 }
3040/* L127: */
3041 }
3042 goto L132;
3043L128:
3044 i__1 = ipph;
3045 for (j = 2; j <= i__1; ++j) {
3046 jc = ipp2 - j;
3047 i__2 = *ido;
3048 for (i__ = 3; i__ <= i__2; i__ += 2) {
3049 i__3 = *l1;
3050 for (k = 1; k <= i__3; ++k) {
3051 ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1] = c1[i__ - 1 + (k +
3052 j * c1_dim2) * c1_dim1] - c1[i__ + (k + jc * c1_dim2)
3053 * c1_dim1];
3054 ch[i__ - 1 + (k + jc * ch_dim2) * ch_dim1] = c1[i__ - 1 + (k
3055 + j * c1_dim2) * c1_dim1] + c1[i__ + (k + jc *
3056 c1_dim2) * c1_dim1];
3057 ch[i__ + (k + j * ch_dim2) * ch_dim1] = c1[i__ + (k + j *
3058 c1_dim2) * c1_dim1] + c1[i__ - 1 + (k + jc * c1_dim2)
3059 * c1_dim1];
3060 ch[i__ + (k + jc * ch_dim2) * ch_dim1] = c1[i__ + (k + j *
3061 c1_dim2) * c1_dim1] - c1[i__ - 1 + (k + jc * c1_dim2)
3062 * c1_dim1];
3063/* L129: */
3064 }
3065/* L130: */
3066 }
3067/* L131: */
3068 }
3069L132:
3070 if (*ido == 1) {
3071 return;
3072 }
3073 i__1 = *idl1;
3074 for (ik = 1; ik <= i__1; ++ik) {
3075 c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
3076/* L133: */
3077 }
3078 i__1 = *ip;
3079 for (j = 2; j <= i__1; ++j) {
3080 i__2 = *l1;
3081 for (k = 1; k <= i__2; ++k) {
3082 c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) *
3083 ch_dim1 + 1];
3084/* L134: */
3085 }
3086/* L135: */
3087 }
3088 if (nbd > *l1) {
3089 goto L139;
3090 }
3091 is = -(*ido);
3092 i__1 = *ip;
3093 for (j = 2; j <= i__1; ++j) {
3094 is += *ido;
3095 idij = is;
3096 i__2 = *ido;
3097 for (i__ = 3; i__ <= i__2; i__ += 2) {
3098 idij += 2;
3099 i__3 = *l1;
3100 for (k = 1; k <= i__3; ++k) {
3101 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[
3102 i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] *
3103 ch[i__ + (k + j * ch_dim2) * ch_dim1];
3104 c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__
3105 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ -
3106 1 + (k + j * ch_dim2) * ch_dim1];
3107/* L136: */
3108 }
3109/* L137: */
3110 }
3111/* L138: */
3112 }
3113 goto L143;
3114L139:
3115 is = -(*ido);
3116 i__1 = *ip;
3117 for (j = 2; j <= i__1; ++j) {
3118 is += *ido;
3119 i__2 = *l1;
3120 for (k = 1; k <= i__2; ++k) {
3121 idij = is;
3122 i__3 = *ido;
3123 for (i__ = 3; i__ <= i__3; i__ += 2) {
3124 idij += 2;
3125 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[
3126 i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] *
3127 ch[i__ + (k + j * ch_dim2) * ch_dim1];
3128 c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__
3129 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ -
3130 1 + (k + j * ch_dim2) * ch_dim1];
3131/* L140: */
3132 }
3133/* L141: */
3134 }
3135/* L142: */
3136 }
3137L143:
3138 return;
3139} /* radbg_ */
3140
3141/* Subroutine */ static void radf2(integer_t *ido, integer_t *l1, real_t *cc,
3142 real_t *ch, real_t *wa1)
3143{
3144 /* System generated locals */
3145 integer_t ch_dim1, ch_offset, cc_dim1, cc_dim2, cc_offset, i__1, i__2;
3146
3147 /* Local variables */
3148 integer_t i__, k, ic;
3149 real_t ti2, tr2;
3150 integer_t idp2;
3151
3152 /* Parameter adjustments */
3153 ch_dim1 = *ido;
3154 ch_offset = 1 + ch_dim1 * 3;
3155 ch -= ch_offset;
3156 cc_dim1 = *ido;
3157 cc_dim2 = *l1;
3158 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
3159 cc -= cc_offset;
3160 --wa1;
3161
3162 /* Function Body */
3163 i__1 = *l1;
3164 for (k = 1; k <= i__1; ++k) {
3165 ch[((k << 1) + 1) * ch_dim1 + 1] = cc[(k + cc_dim2) * cc_dim1 + 1] +
3166 cc[(k + (cc_dim2 << 1)) * cc_dim1 + 1];
3167 ch[*ido + ((k << 1) + 2) * ch_dim1] = cc[(k + cc_dim2) * cc_dim1 + 1]
3168 - cc[(k + (cc_dim2 << 1)) * cc_dim1 + 1];
3169/* L101: */
3170 }
3171 if ((i__1 = *ido - 2) < 0) {
3172 goto L107;
3173 } else if (i__1 == 0) {
3174 goto L105;
3175 } else {
3176 goto L102;
3177 }
3178L102:
3179 idp2 = *ido + 2;
3180 i__1 = *l1;
3181 for (k = 1; k <= i__1; ++k) {
3182 i__2 = *ido;
3183 for (i__ = 3; i__ <= i__2; i__ += 2) {
3184 ic = idp2 - i__;
3185 tr2 = wa1[i__ - 2] * cc[i__ - 1 + (k + (cc_dim2 << 1)) * cc_dim1]
3186 + wa1[i__ - 1] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1];
3187 ti2 = wa1[i__ - 2] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1] -
3188 wa1[i__ - 1] * cc[i__ - 1 + (k + (cc_dim2 << 1)) *
3189 cc_dim1];
3190 ch[i__ + ((k << 1) + 1) * ch_dim1] = cc[i__ + (k + cc_dim2) *
3191 cc_dim1] + ti2;
3192 ch[ic + ((k << 1) + 2) * ch_dim1] = ti2 - cc[i__ + (k + cc_dim2) *
3193 cc_dim1];
3194 ch[i__ - 1 + ((k << 1) + 1) * ch_dim1] = cc[i__ - 1 + (k +
3195 cc_dim2) * cc_dim1] + tr2;
3196 ch[ic - 1 + ((k << 1) + 2) * ch_dim1] = cc[i__ - 1 + (k + cc_dim2)
3197 * cc_dim1] - tr2;
3198/* L103: */
3199 }
3200/* L104: */
3201 }
3202 if (*ido % 2 == 1) {
3203 return;
3204 }
3205L105:
3206 i__1 = *l1;
3207 for (k = 1; k <= i__1; ++k) {
3208 ch[((k << 1) + 2) * ch_dim1 + 1] = -cc[*ido + (k + (cc_dim2 << 1)) *
3209 cc_dim1];
3210 ch[*ido + ((k << 1) + 1) * ch_dim1] = cc[*ido + (k + cc_dim2) *
3211 cc_dim1];
3212/* L106: */
3213 }
3214L107:
3215 return;
3216} /* radf2_ */
3217
3218/* Subroutine */ static void radf3(integer_t *ido, integer_t *l1, real_t *cc,
3219 real_t *ch, real_t *wa1, real_t *wa2)
3220{
3221 /* Initialized data */
3222
3223 static real_t taur = REAL_CONSTANT(-0.5);
3224 static real_t taui =
3225 REAL_CONSTANT(0.8660254037844386467637231707529361834710262690519031402790348975);
3226
3227 /* System generated locals */
3228 integer_t ch_dim1, ch_offset, cc_dim1, cc_dim2, cc_offset, i__1, i__2;
3229
3230 /* Local variables */
3231 integer_t i__, k, ic;
3232 real_t ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
3233 integer_t idp2;
3234
3235 /* Parameter adjustments */
3236 ch_dim1 = *ido;
3237 ch_offset = 1 + (ch_dim1 << 2);
3238 ch -= ch_offset;
3239 cc_dim1 = *ido;
3240 cc_dim2 = *l1;
3241 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
3242 cc -= cc_offset;
3243 --wa1;
3244 --wa2;
3245
3246 /* Function Body */
3247 i__1 = *l1;
3248 for (k = 1; k <= i__1; ++k) {
3249 cr2 = cc[(k + (cc_dim2 << 1)) * cc_dim1 + 1] + cc[(k + cc_dim2 * 3) *
3250 cc_dim1 + 1];
3251 ch[(k * 3 + 1) * ch_dim1 + 1] = cc[(k + cc_dim2) * cc_dim1 + 1] + cr2;
3252 ch[(k * 3 + 3) * ch_dim1 + 1] = taui * (cc[(k + cc_dim2 * 3) *
3253 cc_dim1 + 1] - cc[(k + (cc_dim2 << 1)) * cc_dim1 + 1]);
3254 ch[*ido + (k * 3 + 2) * ch_dim1] = cc[(k + cc_dim2) * cc_dim1 + 1] +
3255 taur * cr2;
3256/* L101: */
3257 }
3258 if (*ido == 1) {
3259 return;
3260 }
3261 idp2 = *ido + 2;
3262 i__1 = *l1;
3263 for (k = 1; k <= i__1; ++k) {
3264 i__2 = *ido;
3265 for (i__ = 3; i__ <= i__2; i__ += 2) {
3266 ic = idp2 - i__;
3267 dr2 = wa1[i__ - 2] * cc[i__ - 1 + (k + (cc_dim2 << 1)) * cc_dim1]
3268 + wa1[i__ - 1] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1];
3269 di2 = wa1[i__ - 2] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1] -
3270 wa1[i__ - 1] * cc[i__ - 1 + (k + (cc_dim2 << 1)) *
3271 cc_dim1];
3272 dr3 = wa2[i__ - 2] * cc[i__ - 1 + (k + cc_dim2 * 3) * cc_dim1] +
3273 wa2[i__ - 1] * cc[i__ + (k + cc_dim2 * 3) * cc_dim1];
3274 di3 = wa2[i__ - 2] * cc[i__ + (k + cc_dim2 * 3) * cc_dim1] - wa2[
3275 i__ - 1] * cc[i__ - 1 + (k + cc_dim2 * 3) * cc_dim1];
3276 cr2 = dr2 + dr3;
3277 ci2 = di2 + di3;
3278 ch[i__ - 1 + (k * 3 + 1) * ch_dim1] = cc[i__ - 1 + (k + cc_dim2) *
3279 cc_dim1] + cr2;
3280 ch[i__ + (k * 3 + 1) * ch_dim1] = cc[i__ + (k + cc_dim2) *
3281 cc_dim1] + ci2;
3282 tr2 = cc[i__ - 1 + (k + cc_dim2) * cc_dim1] + taur * cr2;
3283 ti2 = cc[i__ + (k + cc_dim2) * cc_dim1] + taur * ci2;
3284 tr3 = taui * (di2 - di3);
3285 ti3 = taui * (dr3 - dr2);
3286 ch[i__ - 1 + (k * 3 + 3) * ch_dim1] = tr2 + tr3;
3287 ch[ic - 1 + (k * 3 + 2) * ch_dim1] = tr2 - tr3;
3288 ch[i__ + (k * 3 + 3) * ch_dim1] = ti2 + ti3;
3289 ch[ic + (k * 3 + 2) * ch_dim1] = ti3 - ti2;
3290/* L102: */
3291 }
3292/* L103: */
3293 }
3294 return;
3295} /* radf3_ */
3296
3297/* Subroutine */ static void radf4(integer_t *ido, integer_t *l1, real_t *cc,
3298 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3)
3299{
3300 /* Initialized data */
3301
3302 static real_t hsqt2 =
3303 REAL_CONSTANT(0.70710678118654752440084436210484903928483593768474036588339869);
3304
3305 /* System generated locals */
3306 integer_t cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_offset, i__1, i__2;
3307
3308 /* Local variables */
3309 integer_t i__, k, ic;
3310 real_t ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
3311 tr3, tr4;
3312 integer_t idp2;
3313
3314 /* Parameter adjustments */
3315 ch_dim1 = *ido;
3316 ch_offset = 1 + ch_dim1 * 5;
3317 ch -= ch_offset;
3318 cc_dim1 = *ido;
3319 cc_dim2 = *l1;
3320 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
3321 cc -= cc_offset;
3322 --wa1;
3323 --wa2;
3324 --wa3;
3325
3326 /* Function Body */
3327 i__1 = *l1;
3328 for (k = 1; k <= i__1; ++k) {
3329 tr1 = cc[(k + (cc_dim2 << 1)) * cc_dim1 + 1] + cc[(k + (cc_dim2 << 2))
3330 * cc_dim1 + 1];
3331 tr2 = cc[(k + cc_dim2) * cc_dim1 + 1] + cc[(k + cc_dim2 * 3) *
3332 cc_dim1 + 1];
3333 ch[((k << 2) + 1) * ch_dim1 + 1] = tr1 + tr2;
3334 ch[*ido + ((k << 2) + 4) * ch_dim1] = tr2 - tr1;
3335 ch[*ido + ((k << 2) + 2) * ch_dim1] = cc[(k + cc_dim2) * cc_dim1 + 1]
3336 - cc[(k + cc_dim2 * 3) * cc_dim1 + 1];
3337 ch[((k << 2) + 3) * ch_dim1 + 1] = cc[(k + (cc_dim2 << 2)) * cc_dim1
3338 + 1] - cc[(k + (cc_dim2 << 1)) * cc_dim1 + 1];
3339/* L101: */
3340 }
3341 if ((i__1 = *ido - 2) < 0) {
3342 goto L107;
3343 } else if (i__1 == 0) {
3344 goto L105;
3345 } else {
3346 goto L102;
3347 }
3348L102:
3349 idp2 = *ido + 2;
3350 i__1 = *l1;
3351 for (k = 1; k <= i__1; ++k) {
3352 i__2 = *ido;
3353 for (i__ = 3; i__ <= i__2; i__ += 2) {
3354 ic = idp2 - i__;
3355 cr2 = wa1[i__ - 2] * cc[i__ - 1 + (k + (cc_dim2 << 1)) * cc_dim1]
3356 + wa1[i__ - 1] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1];
3357 ci2 = wa1[i__ - 2] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1] -
3358 wa1[i__ - 1] * cc[i__ - 1 + (k + (cc_dim2 << 1)) *
3359 cc_dim1];
3360 cr3 = wa2[i__ - 2] * cc[i__ - 1 + (k + cc_dim2 * 3) * cc_dim1] +
3361 wa2[i__ - 1] * cc[i__ + (k + cc_dim2 * 3) * cc_dim1];
3362 ci3 = wa2[i__ - 2] * cc[i__ + (k + cc_dim2 * 3) * cc_dim1] - wa2[
3363 i__ - 1] * cc[i__ - 1 + (k + cc_dim2 * 3) * cc_dim1];
3364 cr4 = wa3[i__ - 2] * cc[i__ - 1 + (k + (cc_dim2 << 2)) * cc_dim1]
3365 + wa3[i__ - 1] * cc[i__ + (k + (cc_dim2 << 2)) * cc_dim1];
3366 ci4 = wa3[i__ - 2] * cc[i__ + (k + (cc_dim2 << 2)) * cc_dim1] -
3367 wa3[i__ - 1] * cc[i__ - 1 + (k + (cc_dim2 << 2)) *
3368 cc_dim1];
3369 tr1 = cr2 + cr4;
3370 tr4 = cr4 - cr2;
3371 ti1 = ci2 + ci4;
3372 ti4 = ci2 - ci4;
3373 ti2 = cc[i__ + (k + cc_dim2) * cc_dim1] + ci3;
3374 ti3 = cc[i__ + (k + cc_dim2) * cc_dim1] - ci3;
3375 tr2 = cc[i__ - 1 + (k + cc_dim2) * cc_dim1] + cr3;
3376 tr3 = cc[i__ - 1 + (k + cc_dim2) * cc_dim1] - cr3;
3377 ch[i__ - 1 + ((k << 2) + 1) * ch_dim1] = tr1 + tr2;
3378 ch[ic - 1 + ((k << 2) + 4) * ch_dim1] = tr2 - tr1;
3379 ch[i__ + ((k << 2) + 1) * ch_dim1] = ti1 + ti2;
3380 ch[ic + ((k << 2) + 4) * ch_dim1] = ti1 - ti2;
3381 ch[i__ - 1 + ((k << 2) + 3) * ch_dim1] = ti4 + tr3;
3382 ch[ic - 1 + ((k << 2) + 2) * ch_dim1] = tr3 - ti4;
3383 ch[i__ + ((k << 2) + 3) * ch_dim1] = tr4 + ti3;
3384 ch[ic + ((k << 2) + 2) * ch_dim1] = tr4 - ti3;
3385/* L103: */
3386 }
3387/* L104: */
3388 }
3389 if (*ido % 2 == 1) {
3390 return;
3391 }
3392L105:
3393 i__1 = *l1;
3394 for (k = 1; k <= i__1; ++k) {
3395 ti1 = -hsqt2 * (cc[*ido + (k + (cc_dim2 << 1)) * cc_dim1] + cc[*ido +
3396 (k + (cc_dim2 << 2)) * cc_dim1]);
3397 tr1 = hsqt2 * (cc[*ido + (k + (cc_dim2 << 1)) * cc_dim1] - cc[*ido + (
3398 k + (cc_dim2 << 2)) * cc_dim1]);
3399 ch[*ido + ((k << 2) + 1) * ch_dim1] = tr1 + cc[*ido + (k + cc_dim2) *
3400 cc_dim1];
3401 ch[*ido + ((k << 2) + 3) * ch_dim1] = cc[*ido + (k + cc_dim2) *
3402 cc_dim1] - tr1;
3403 ch[((k << 2) + 2) * ch_dim1 + 1] = ti1 - cc[*ido + (k + cc_dim2 * 3) *
3404 cc_dim1];
3405 ch[((k << 2) + 4) * ch_dim1 + 1] = ti1 + cc[*ido + (k + cc_dim2 * 3) *
3406 cc_dim1];
3407/* L106: */
3408 }
3409L107:
3410 return;
3411} /* radf4_ */
3412
3413/* Subroutine */ static void radf5(integer_t *ido, integer_t *l1, real_t *cc,
3414 real_t *ch, real_t *wa1, real_t *wa2, real_t *wa3,
3415 real_t *wa4)
3416{
3417 /* Initialized data */
3418
3419 static real_t tr11 =
3420 REAL_CONSTANT(0.3090169943749474241022934171828195886015458990288143106772431137);
3421 static real_t ti11 =
3422 REAL_CONSTANT(0.9510565162951535721164393337938214340569863412575022244730564442);
3423 static real_t tr12 =
3424 REAL_CONSTANT(-0.8090169943749474241022934171828190588601545899028814310677431135);
3425 static real_t ti12 =
3426 REAL_CONSTANT(0.5877852522924731291687059546390727685976524376431459107227248076);
3427
3428 /* System generated locals */
3429 integer_t cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_offset, i__1, i__2;
3430
3431 /* Local variables */
3432 integer_t i__, k, ic;
3433 real_t ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3,
3434 dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
3435 integer_t idp2;
3436
3437 /* Parameter adjustments */
3438 ch_dim1 = *ido;
3439 ch_offset = 1 + ch_dim1 * 6;
3440 ch -= ch_offset;
3441 cc_dim1 = *ido;
3442 cc_dim2 = *l1;
3443 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
3444 cc -= cc_offset;
3445 --wa1;
3446 --wa2;
3447 --wa3;
3448 --wa4;
3449
3450 /* Function Body */
3451 i__1 = *l1;
3452 for (k = 1; k <= i__1; ++k) {
3453 cr2 = cc[(k + cc_dim2 * 5) * cc_dim1 + 1] + cc[(k + (cc_dim2 << 1)) *
3454 cc_dim1 + 1];
3455 ci5 = cc[(k + cc_dim2 * 5) * cc_dim1 + 1] - cc[(k + (cc_dim2 << 1)) *
3456 cc_dim1 + 1];
3457 cr3 = cc[(k + (cc_dim2 << 2)) * cc_dim1 + 1] + cc[(k + cc_dim2 * 3) *
3458 cc_dim1 + 1];
3459 ci4 = cc[(k + (cc_dim2 << 2)) * cc_dim1 + 1] - cc[(k + cc_dim2 * 3) *
3460 cc_dim1 + 1];
3461 ch[(k * 5 + 1) * ch_dim1 + 1] = cc[(k + cc_dim2) * cc_dim1 + 1] + cr2
3462 + cr3;
3463 ch[*ido + (k * 5 + 2) * ch_dim1] = cc[(k + cc_dim2) * cc_dim1 + 1] +
3464 tr11 * cr2 + tr12 * cr3;
3465 ch[(k * 5 + 3) * ch_dim1 + 1] = ti11 * ci5 + ti12 * ci4;
3466 ch[*ido + (k * 5 + 4) * ch_dim1] = cc[(k + cc_dim2) * cc_dim1 + 1] +
3467 tr12 * cr2 + tr11 * cr3;
3468 ch[(k * 5 + 5) * ch_dim1 + 1] = ti12 * ci5 - ti11 * ci4;
3469/* L101: */
3470 }
3471 if (*ido == 1) {
3472 return;
3473 }
3474 idp2 = *ido + 2;
3475 i__1 = *l1;
3476 for (k = 1; k <= i__1; ++k) {
3477 i__2 = *ido;
3478 for (i__ = 3; i__ <= i__2; i__ += 2) {
3479 ic = idp2 - i__;
3480 dr2 = wa1[i__ - 2] * cc[i__ - 1 + (k + (cc_dim2 << 1)) * cc_dim1]
3481 + wa1[i__ - 1] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1];
3482 di2 = wa1[i__ - 2] * cc[i__ + (k + (cc_dim2 << 1)) * cc_dim1] -
3483 wa1[i__ - 1] * cc[i__ - 1 + (k + (cc_dim2 << 1)) *
3484 cc_dim1];
3485 dr3 = wa2[i__ - 2] * cc[i__ - 1 + (k + cc_dim2 * 3) * cc_dim1] +
3486 wa2[i__ - 1] * cc[i__ + (k + cc_dim2 * 3) * cc_dim1];
3487 di3 = wa2[i__ - 2] * cc[i__ + (k + cc_dim2 * 3) * cc_dim1] - wa2[
3488 i__ - 1] * cc[i__ - 1 + (k + cc_dim2 * 3) * cc_dim1];
3489 dr4 = wa3[i__ - 2] * cc[i__ - 1 + (k + (cc_dim2 << 2)) * cc_dim1]
3490 + wa3[i__ - 1] * cc[i__ + (k + (cc_dim2 << 2)) * cc_dim1];
3491 di4 = wa3[i__ - 2] * cc[i__ + (k + (cc_dim2 << 2)) * cc_dim1] -
3492 wa3[i__ - 1] * cc[i__ - 1 + (k + (cc_dim2 << 2)) *
3493 cc_dim1];
3494 dr5 = wa4[i__ - 2] * cc[i__ - 1 + (k + cc_dim2 * 5) * cc_dim1] +
3495 wa4[i__ - 1] * cc[i__ + (k + cc_dim2 * 5) * cc_dim1];
3496 di5 = wa4[i__ - 2] * cc[i__ + (k + cc_dim2 * 5) * cc_dim1] - wa4[
3497 i__ - 1] * cc[i__ - 1 + (k + cc_dim2 * 5) * cc_dim1];
3498 cr2 = dr2 + dr5;
3499 ci5 = dr5 - dr2;
3500 cr5 = di2 - di5;
3501 ci2 = di2 + di5;
3502 cr3 = dr3 + dr4;
3503 ci4 = dr4 - dr3;
3504 cr4 = di3 - di4;
3505 ci3 = di3 + di4;
3506 ch[i__ - 1 + (k * 5 + 1) * ch_dim1] = cc[i__ - 1 + (k + cc_dim2) *
3507 cc_dim1] + cr2 + cr3;
3508 ch[i__ + (k * 5 + 1) * ch_dim1] = cc[i__ + (k + cc_dim2) *
3509 cc_dim1] + ci2 + ci3;
3510 tr2 = cc[i__ - 1 + (k + cc_dim2) * cc_dim1] + tr11 * cr2 + tr12 *
3511 cr3;
3512 ti2 = cc[i__ + (k + cc_dim2) * cc_dim1] + tr11 * ci2 + tr12 * ci3;
3513 tr3 = cc[i__ - 1 + (k + cc_dim2) * cc_dim1] + tr12 * cr2 + tr11 *
3514 cr3;
3515 ti3 = cc[i__ + (k + cc_dim2) * cc_dim1] + tr12 * ci2 + tr11 * ci3;
3516 tr5 = ti11 * cr5 + ti12 * cr4;
3517 ti5 = ti11 * ci5 + ti12 * ci4;
3518 tr4 = ti12 * cr5 - ti11 * cr4;
3519 ti4 = ti12 * ci5 - ti11 * ci4;
3520 ch[i__ - 1 + (k * 5 + 3) * ch_dim1] = tr2 + tr5;
3521 ch[ic - 1 + (k * 5 + 2) * ch_dim1] = tr2 - tr5;
3522 ch[i__ + (k * 5 + 3) * ch_dim1] = ti2 + ti5;
3523 ch[ic + (k * 5 + 2) * ch_dim1] = ti5 - ti2;
3524 ch[i__ - 1 + (k * 5 + 5) * ch_dim1] = tr3 + tr4;
3525 ch[ic - 1 + (k * 5 + 4) * ch_dim1] = tr3 - tr4;
3526 ch[i__ + (k * 5 + 5) * ch_dim1] = ti3 + ti4;
3527 ch[ic + (k * 5 + 4) * ch_dim1] = ti4 - ti3;
3528/* L102: */
3529 }
3530/* L103: */
3531 }
3532 return;
3533} /* radf5_ */
3534
3535/* Subroutine */ static void radfg(integer_t *ido, integer_t *ip, integer_t *l1, integer_t *
3536 idl1, real_t *cc, real_t *c1, real_t *c2, real_t *ch,
3537 real_t *ch2, real_t *wa)
3538{
3539 /* Initialized data */
3540
3541 static real_t tpi =
3542 REAL_CONSTANT(6.283185307179586476925286766559005768394338798750116419498891846);
3543
3544 /* System generated locals */
3545 integer_t ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
3546 c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset,
3547 i__1, i__2, i__3;
3548
3549 /* Local variables */
3550 integer_t i__, j, k, l, j2, ic, jc, lc, ik, is;
3551 real_t dc2, ai1, ai2, ar1, ar2, ds2;
3552 integer_t nbd;
3553 real_t dcp, arg, dsp, ar1h, ar2h;
3554 integer_t idp2, ipp2, idij, ipph;
3555
3556 /* Parameter adjustments */
3557 ch_dim1 = *ido;
3558 ch_dim2 = *l1;
3559 ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
3560 ch -= ch_offset;
3561 c1_dim1 = *ido;
3562 c1_dim2 = *l1;
3563 c1_offset = 1 + c1_dim1 * (1 + c1_dim2);
3564 c1 -= c1_offset;
3565 cc_dim1 = *ido;
3566 cc_dim2 = *ip;
3567 cc_offset = 1 + cc_dim1 * (1 + cc_dim2);
3568 cc -= cc_offset;
3569 ch2_dim1 = *idl1;
3570 ch2_offset = 1 + ch2_dim1;
3571 ch2 -= ch2_offset;
3572 c2_dim1 = *idl1;
3573 c2_offset = 1 + c2_dim1;
3574 c2 -= c2_offset;
3575 --wa;
3576
3577 /* Function Body */
3578 arg = tpi / (real_t) (*ip);
3579 dcp = cos(arg);
3580 dsp = sin(arg);
3581 ipph = (*ip + 1) / 2;
3582 ipp2 = *ip + 2;
3583 idp2 = *ido + 2;
3584 nbd = (*ido - 1) / 2;
3585 if (*ido == 1) {
3586 goto L119;
3587 }
3588 i__1 = *idl1;
3589 for (ik = 1; ik <= i__1; ++ik) {
3590 ch2[ik + ch2_dim1] = c2[ik + c2_dim1];
3591/* L101: */
3592 }
3593 i__1 = *ip;
3594 for (j = 2; j <= i__1; ++j) {
3595 i__2 = *l1;
3596 for (k = 1; k <= i__2; ++k) {
3597 ch[(k + j * ch_dim2) * ch_dim1 + 1] = c1[(k + j * c1_dim2) *
3598 c1_dim1 + 1];
3599/* L102: */
3600 }
3601/* L103: */
3602 }
3603 if (nbd > *l1) {
3604 goto L107;
3605 }
3606 is = -(*ido);
3607 i__1 = *ip;
3608 for (j = 2; j <= i__1; ++j) {
3609 is += *ido;
3610 idij = is;
3611 i__2 = *ido;
3612 for (i__ = 3; i__ <= i__2; i__ += 2) {
3613 idij += 2;
3614 i__3 = *l1;
3615 for (k = 1; k <= i__3; ++k) {
3616 ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1] = wa[idij - 1] * c1[
3617 i__ - 1 + (k + j * c1_dim2) * c1_dim1] + wa[idij] *
3618 c1[i__ + (k + j * c1_dim2) * c1_dim1];
3619 ch[i__ + (k + j * ch_dim2) * ch_dim1] = wa[idij - 1] * c1[i__
3620 + (k + j * c1_dim2) * c1_dim1] - wa[idij] * c1[i__ -
3621 1 + (k + j * c1_dim2) * c1_dim1];
3622/* L104: */
3623 }
3624/* L105: */
3625 }
3626/* L106: */
3627 }
3628 goto L111;
3629L107:
3630 is = -(*ido);
3631 i__1 = *ip;
3632 for (j = 2; j <= i__1; ++j) {
3633 is += *ido;
3634 i__2 = *l1;
3635 for (k = 1; k <= i__2; ++k) {
3636 idij = is;
3637 i__3 = *ido;
3638 for (i__ = 3; i__ <= i__3; i__ += 2) {
3639 idij += 2;
3640 ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1] = wa[idij - 1] * c1[
3641 i__ - 1 + (k + j * c1_dim2) * c1_dim1] + wa[idij] *
3642 c1[i__ + (k + j * c1_dim2) * c1_dim1];
3643 ch[i__ + (k + j * ch_dim2) * ch_dim1] = wa[idij - 1] * c1[i__
3644 + (k + j * c1_dim2) * c1_dim1] - wa[idij] * c1[i__ -
3645 1 + (k + j * c1_dim2) * c1_dim1];
3646/* L108: */
3647 }
3648/* L109: */
3649 }
3650/* L110: */
3651 }
3652L111:
3653 if (nbd < *l1) {
3654 goto L115;
3655 }
3656 i__1 = ipph;
3657 for (j = 2; j <= i__1; ++j) {
3658 jc = ipp2 - j;
3659 i__2 = *l1;
3660 for (k = 1; k <= i__2; ++k) {
3661 i__3 = *ido;
3662 for (i__ = 3; i__ <= i__3; i__ += 2) {
3663 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = ch[i__ - 1 + (k +
3664 j * ch_dim2) * ch_dim1] + ch[i__ - 1 + (k + jc *
3665 ch_dim2) * ch_dim1];
3666 c1[i__ - 1 + (k + jc * c1_dim2) * c1_dim1] = ch[i__ + (k + j *
3667 ch_dim2) * ch_dim1] - ch[i__ + (k + jc * ch_dim2) *
3668 ch_dim1];
3669 c1[i__ + (k + j * c1_dim2) * c1_dim1] = ch[i__ + (k + j *
3670 ch_dim2) * ch_dim1] + ch[i__ + (k + jc * ch_dim2) *
3671 ch_dim1];
3672 c1[i__ + (k + jc * c1_dim2) * c1_dim1] = ch[i__ - 1 + (k + jc
3673 * ch_dim2) * ch_dim1] - ch[i__ - 1 + (k + j * ch_dim2)
3674 * ch_dim1];
3675/* L112: */
3676 }
3677/* L113: */
3678 }
3679/* L114: */
3680 }
3681 goto L121;
3682L115:
3683 i__1 = ipph;
3684 for (j = 2; j <= i__1; ++j) {
3685 jc = ipp2 - j;
3686 i__2 = *ido;
3687 for (i__ = 3; i__ <= i__2; i__ += 2) {
3688 i__3 = *l1;
3689 for (k = 1; k <= i__3; ++k) {
3690 c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = ch[i__ - 1 + (k +
3691 j * ch_dim2) * ch_dim1] + ch[i__ - 1 + (k + jc *
3692 ch_dim2) * ch_dim1];
3693 c1[i__ - 1 + (k + jc * c1_dim2) * c1_dim1] = ch[i__ + (k + j *
3694 ch_dim2) * ch_dim1] - ch[i__ + (k + jc * ch_dim2) *
3695 ch_dim1];
3696 c1[i__ + (k + j * c1_dim2) * c1_dim1] = ch[i__ + (k + j *
3697 ch_dim2) * ch_dim1] + ch[i__ + (k + jc * ch_dim2) *
3698 ch_dim1];
3699 c1[i__ + (k + jc * c1_dim2) * c1_dim1] = ch[i__ - 1 + (k + jc
3700 * ch_dim2) * ch_dim1] - ch[i__ - 1 + (k + j * ch_dim2)
3701 * ch_dim1];
3702/* L116: */
3703 }
3704/* L117: */
3705 }
3706/* L118: */
3707 }
3708 goto L121;
3709L119:
3710 i__1 = *idl1;
3711 for (ik = 1; ik <= i__1; ++ik) {
3712 c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
3713/* L120: */
3714 }
3715L121:
3716 i__1 = ipph;
3717 for (j = 2; j <= i__1; ++j) {
3718 jc = ipp2 - j;
3719 i__2 = *l1;
3720 for (k = 1; k <= i__2; ++k) {
3721 c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) *
3722 ch_dim1 + 1] + ch[(k + jc * ch_dim2) * ch_dim1 + 1];
3723 c1[(k + jc * c1_dim2) * c1_dim1 + 1] = ch[(k + jc * ch_dim2) *
3724 ch_dim1 + 1] - ch[(k + j * ch_dim2) * ch_dim1 + 1];
3725/* L122: */
3726 }
3727/* L123: */
3728 }
3729
3730 ar1 = REAL_CONSTANT(1.0);
3731 ai1 = REAL_CONSTANT(0.0);
3732 i__1 = ipph;
3733 for (l = 2; l <= i__1; ++l) {
3734 lc = ipp2 - l;
3735 ar1h = dcp * ar1 - dsp * ai1;
3736 ai1 = dcp * ai1 + dsp * ar1;
3737 ar1 = ar1h;
3738 i__2 = *idl1;
3739 for (ik = 1; ik <= i__2; ++ik) {
3740 ch2[ik + l * ch2_dim1] = c2[ik + c2_dim1] + ar1 * c2[ik + (
3741 c2_dim1 << 1)];
3742 ch2[ik + lc * ch2_dim1] = ai1 * c2[ik + *ip * c2_dim1];
3743/* L124: */
3744 }
3745 dc2 = ar1;
3746 ds2 = ai1;
3747 ar2 = ar1;
3748 ai2 = ai1;
3749 i__2 = ipph;
3750 for (j = 3; j <= i__2; ++j) {
3751 jc = ipp2 - j;
3752 ar2h = dc2 * ar2 - ds2 * ai2;
3753 ai2 = dc2 * ai2 + ds2 * ar2;
3754 ar2 = ar2h;
3755 i__3 = *idl1;
3756 for (ik = 1; ik <= i__3; ++ik) {
3757 ch2[ik + l * ch2_dim1] += ar2 * c2[ik + j * c2_dim1];
3758 ch2[ik + lc * ch2_dim1] += ai2 * c2[ik + jc * c2_dim1];
3759/* L125: */
3760 }
3761/* L126: */
3762 }
3763/* L127: */
3764 }
3765 i__1 = ipph;
3766 for (j = 2; j <= i__1; ++j) {
3767 i__2 = *idl1;
3768 for (ik = 1; ik <= i__2; ++ik) {
3769 ch2[ik + ch2_dim1] += c2[ik + j * c2_dim1];
3770/* L128: */
3771 }
3772/* L129: */
3773 }
3774
3775 if (*ido < *l1) {
3776 goto L132;
3777 }
3778 i__1 = *l1;
3779 for (k = 1; k <= i__1; ++k) {
3780 i__2 = *ido;
3781 for (i__ = 1; i__ <= i__2; ++i__) {
3782 cc[i__ + (k * cc_dim2 + 1) * cc_dim1] = ch[i__ + (k + ch_dim2) *
3783 ch_dim1];
3784/* L130: */
3785 }
3786/* L131: */
3787 }
3788 goto L135;
3789L132:
3790 i__1 = *ido;
3791 for (i__ = 1; i__ <= i__1; ++i__) {
3792 i__2 = *l1;
3793 for (k = 1; k <= i__2; ++k) {
3794 cc[i__ + (k * cc_dim2 + 1) * cc_dim1] = ch[i__ + (k + ch_dim2) *
3795 ch_dim1];
3796/* L133: */
3797 }
3798/* L134: */
3799 }
3800L135:
3801 i__1 = ipph;
3802 for (j = 2; j <= i__1; ++j) {
3803 jc = ipp2 - j;
3804 j2 = j + j;
3805 i__2 = *l1;
3806 for (k = 1; k <= i__2; ++k) {
3807 cc[*ido + (j2 - 2 + k * cc_dim2) * cc_dim1] = ch[(k + j * ch_dim2)
3808 * ch_dim1 + 1];
3809 cc[(j2 - 1 + k * cc_dim2) * cc_dim1 + 1] = ch[(k + jc * ch_dim2) *
3810 ch_dim1 + 1];
3811/* L136: */
3812 }
3813/* L137: */
3814 }
3815 if (*ido == 1) {
3816 return;
3817 }
3818 if (nbd < *l1) {
3819 goto L141;
3820 }
3821 i__1 = ipph;
3822 for (j = 2; j <= i__1; ++j) {
3823 jc = ipp2 - j;
3824 j2 = j + j;
3825 i__2 = *l1;
3826 for (k = 1; k <= i__2; ++k) {
3827 i__3 = *ido;
3828 for (i__ = 3; i__ <= i__3; i__ += 2) {
3829 ic = idp2 - i__;
3830 cc[i__ - 1 + (j2 - 1 + k * cc_dim2) * cc_dim1] = ch[i__ - 1 +
3831 (k + j * ch_dim2) * ch_dim1] + ch[i__ - 1 + (k + jc *
3832 ch_dim2) * ch_dim1];
3833 cc[ic - 1 + (j2 - 2 + k * cc_dim2) * cc_dim1] = ch[i__ - 1 + (
3834 k + j * ch_dim2) * ch_dim1] - ch[i__ - 1 + (k + jc *
3835 ch_dim2) * ch_dim1];
3836 cc[i__ + (j2 - 1 + k * cc_dim2) * cc_dim1] = ch[i__ + (k + j *
3837 ch_dim2) * ch_dim1] + ch[i__ + (k + jc * ch_dim2) *
3838 ch_dim1];
3839 cc[ic + (j2 - 2 + k * cc_dim2) * cc_dim1] = ch[i__ + (k + jc *
3840 ch_dim2) * ch_dim1] - ch[i__ + (k + j * ch_dim2) *
3841 ch_dim1];
3842/* L138: */
3843 }
3844/* L139: */
3845 }
3846/* L140: */
3847 }
3848 return;
3849L141:
3850 i__1 = ipph;
3851 for (j = 2; j <= i__1; ++j) {
3852 jc = ipp2 - j;
3853 j2 = j + j;
3854 i__2 = *ido;
3855 for (i__ = 3; i__ <= i__2; i__ += 2) {
3856 ic = idp2 - i__;
3857 i__3 = *l1;
3858 for (k = 1; k <= i__3; ++k) {
3859 cc[i__ - 1 + (j2 - 1 + k * cc_dim2) * cc_dim1] = ch[i__ - 1 +
3860 (k + j * ch_dim2) * ch_dim1] + ch[i__ - 1 + (k + jc *
3861 ch_dim2) * ch_dim1];
3862 cc[ic - 1 + (j2 - 2 + k * cc_dim2) * cc_dim1] = ch[i__ - 1 + (
3863 k + j * ch_dim2) * ch_dim1] - ch[i__ - 1 + (k + jc *
3864 ch_dim2) * ch_dim1];
3865 cc[i__ + (j2 - 1 + k * cc_dim2) * cc_dim1] = ch[i__ + (k + j *
3866 ch_dim2) * ch_dim1] + ch[i__ + (k + jc * ch_dim2) *
3867 ch_dim1];
3868 cc[ic + (j2 - 2 + k * cc_dim2) * cc_dim1] = ch[i__ + (k + jc *
3869 ch_dim2) * ch_dim1] - ch[i__ + (k + j * ch_dim2) *
3870 ch_dim1];
3871/* L142: */
3872 }
3873/* L143: */
3874 }
3875/* L144: */
3876 }
3877 return;
3878} /* radfg_ */
3879
3880/* Subroutine */ void rfftb(integer_t *n, real_t *r__, real_t *wsave,
3881 integer_t *ifac)
3882{
3883 /* Parameter adjustments */
3884 --ifac;
3885 --wsave;
3886 --r__;
3887
3888 /* Function Body */
3889 if (*n == 1) {
3890 return;
3891 }
3892 rfftb1(n, &r__[1], &wsave[1], &wsave[*n + 1], &ifac[1]);
3893 return;
3894} /* rfftb_ */
3895
3896/* Subroutine */ static void rfftb1(integer_t *n, real_t *c__, real_t *ch,
3897 real_t *wa, integer_t *ifac)
3898{
3899 /* System generated locals */
3900 integer_t i__1;
3901
3902 /* Local variables */
3903 integer_t i__, k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
3904
3905 /* Parameter adjustments */
3906 --ifac;
3907 --wa;
3908 --ch;
3909 --c__;
3910
3911 /* Function Body */
3912 nf = ifac[2];
3913 na = 0;
3914 l1 = 1;
3915 iw = 1;
3916 i__1 = nf;
3917 for (k1 = 1; k1 <= i__1; ++k1) {
3918 ip = ifac[k1 + 2];
3919 l2 = ip * l1;
3920 ido = *n / l2;
3921 idl1 = ido * l1;
3922 if (ip != 4) {
3923 goto L103;
3924 }
3925 ix2 = iw + ido;
3926 ix3 = ix2 + ido;
3927 if (na != 0) {
3928 goto L101;
3929 }
3930 radb4(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
3931 goto L102;
3932L101:
3933 radb4(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3]);
3934L102:
3935 na = 1 - na;
3936 goto L115;
3937L103:
3938 if (ip != 2) {
3939 goto L106;
3940 }
3941 if (na != 0) {
3942 goto L104;
3943 }
3944 radb2(&ido, &l1, &c__[1], &ch[1], &wa[iw]);
3945 goto L105;
3946L104:
3947 radb2(&ido, &l1, &ch[1], &c__[1], &wa[iw]);
3948L105:
3949 na = 1 - na;
3950 goto L115;
3951L106:
3952 if (ip != 3) {
3953 goto L109;
3954 }
3955 ix2 = iw + ido;
3956 if (na != 0) {
3957 goto L107;
3958 }
3959 radb3(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2]);
3960 goto L108;
3961L107:
3962 radb3(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2]);
3963L108:
3964 na = 1 - na;
3965 goto L115;
3966L109:
3967 if (ip != 5) {
3968 goto L112;
3969 }
3970 ix2 = iw + ido;
3971 ix3 = ix2 + ido;
3972 ix4 = ix3 + ido;
3973 if (na != 0) {
3974 goto L110;
3975 }
3976 radb5(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
3977 ix4]);
3978 goto L111;
3979L110:
3980 radb5(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
3981 ix4]);
3982L111:
3983 na = 1 - na;
3984 goto L115;
3985L112:
3986 if (na != 0) {
3987 goto L113;
3988 }
3989 radbg(&ido, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[
3990 1], &wa[iw]);
3991 goto L114;
3992L113:
3993 radbg(&ido, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c__[1], &c__[1]
3994 , &wa[iw]);
3995L114:
3996 if (ido == 1) {
3997 na = 1 - na;
3998 }
3999L115:
4000 l1 = l2;
4001 iw += (ip - 1) * ido;
4002/* L116: */
4003 }
4004 if (na == 0) {
4005 return;
4006 }
4007 i__1 = *n;
4008 for (i__ = 1; i__ <= i__1; ++i__) {
4009 c__[i__] = ch[i__];
4010/* L117: */
4011 }
4012 return;
4013} /* rfftb1_ */
4014
4015/* Subroutine */ void rfftf(integer_t *n, real_t *r__, real_t *wsave,
4016 integer_t *ifac)
4017{
4018 /* Parameter adjustments */
4019 --ifac;
4020 --wsave;
4021 --r__;
4022
4023 /* Function Body */
4024 if (*n == 1) {
4025 return;
4026 }
4027 rfftf1(n, &r__[1], &wsave[1], &wsave[*n + 1], &ifac[1]);
4028 return;
4029} /* rfftf_ */
4030
4031/* Subroutine */ static void rfftf1(integer_t *n, real_t *c__, real_t *ch,
4032 real_t *wa, integer_t *ifac)
4033{
4034 /* System generated locals */
4035 integer_t i__1;
4036
4037 /* Local variables */
4038 integer_t i__, k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
4039
4040 /* Parameter adjustments */
4041 --ifac;
4042 --wa;
4043 --ch;
4044 --c__;
4045
4046 /* Function Body */
4047 nf = ifac[2];
4048 na = 1;
4049 l2 = *n;
4050 iw = *n;
4051 i__1 = nf;
4052 for (k1 = 1; k1 <= i__1; ++k1) {
4053 kh = nf - k1;
4054 ip = ifac[kh + 3];
4055 l1 = l2 / ip;
4056 ido = *n / l2;
4057 idl1 = ido * l1;
4058 iw -= (ip - 1) * ido;
4059 na = 1 - na;
4060 if (ip != 4) {
4061 goto L102;
4062 }
4063 ix2 = iw + ido;
4064 ix3 = ix2 + ido;
4065 if (na != 0) {
4066 goto L101;
4067 }
4068 radf4(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
4069 goto L110;
4070L101:
4071 radf4(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3]);
4072 goto L110;
4073L102:
4074 if (ip != 2) {
4075 goto L104;
4076 }
4077 if (na != 0) {
4078 goto L103;
4079 }
4080 radf2(&ido, &l1, &c__[1], &ch[1], &wa[iw]);
4081 goto L110;
4082L103:
4083 radf2(&ido, &l1, &ch[1], &c__[1], &wa[iw]);
4084 goto L110;
4085L104:
4086 if (ip != 3) {
4087 goto L106;
4088 }
4089 ix2 = iw + ido;
4090 if (na != 0) {
4091 goto L105;
4092 }
4093 radf3(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2]);
4094 goto L110;
4095L105:
4096 radf3(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2]);
4097 goto L110;
4098L106:
4099 if (ip != 5) {
4100 goto L108;
4101 }
4102 ix2 = iw + ido;
4103 ix3 = ix2 + ido;
4104 ix4 = ix3 + ido;
4105 if (na != 0) {
4106 goto L107;
4107 }
4108 radf5(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
4109 ix4]);
4110 goto L110;
4111L107:
4112 radf5(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
4113 ix4]);
4114 goto L110;
4115L108:
4116 if (ido == 1) {
4117 na = 1 - na;
4118 }
4119 if (na != 0) {
4120 goto L109;
4121 }
4122 radfg(&ido, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[
4123 1], &wa[iw]);
4124 na = 1;
4125 goto L110;
4126L109:
4127 radfg(&ido, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c__[1], &c__[1]
4128 , &wa[iw]);
4129 na = 0;
4130L110:
4131 l2 = l1;
4132/* L111: */
4133 }
4134 if (na == 1) {
4135 return;
4136 }
4137 i__1 = *n;
4138 for (i__ = 1; i__ <= i__1; ++i__) {
4139 c__[i__] = ch[i__];
4140/* L112: */
4141 }
4142 return;
4143} /* rfftf1_ */
4144
4145/* Subroutine */ void rffti(integer_t *n, real_t *wsave, integer_t *ifac)
4146{
4147 /* Parameter adjustments */
4148 --ifac;
4149 --wsave;
4150
4151 /* Function Body */
4152 if (*n == 1) {
4153 return;
4154 }
4155 rffti1(n, &wsave[*n + 1], &ifac[1]);
4156 return;
4157} /* rffti_ */
4158
4159/* Subroutine */ static void rffti1(integer_t *n, real_t *wa, integer_t *ifac)
4160{
4161 /* Initialized data */
4162
4163 static integer_t ntryh[4] = { 4,2,3,5 };
4164
4165 /* System generated locals */
4166 integer_t i__1, i__2, i__3;
4167
4168 /* Local variables */
4169 integer_t i__, j, k1, l1, l2, ib;
4170 real_t fi;
4171 integer_t ld, ii, nf, ip, nl, is, nq, nr;
4172 real_t arg;
4173 integer_t ido, ipm;
4174 real_t tpi;
4175 integer_t nfm1;
4176 real_t argh;
4177 integer_t ntry=0;
4178 real_t argld;
4179
4180 /* Parameter adjustments */
4181 --ifac;
4182 --wa;
4183
4184 /* Function Body */
4185 nl = *n;
4186 nf = 0;
4187 j = 0;
4188L101:
4189 ++j;
4190 if (j - 4 <= 0) {
4191 goto L102;
4192 } else {
4193 goto L103;
4194 }
4195L102:
4196 ntry = ntryh[j - 1];
4197 goto L104;
4198L103:
4199 ntry += 2;
4200L104:
4201 nq = nl / ntry;
4202 nr = nl - ntry * nq;
4203 if (nr != 0) {
4204 goto L101;
4205 } else {
4206 goto L105;
4207 }
4208L105:
4209 ++nf;
4210 ifac[nf + 2] = ntry;
4211 nl = nq;
4212 if (ntry != 2) {
4213 goto L107;
4214 }
4215 if (nf == 1) {
4216 goto L107;
4217 }
4218 i__1 = nf;
4219 for (i__ = 2; i__ <= i__1; ++i__) {
4220 ib = nf - i__ + 2;
4221 ifac[ib + 2] = ifac[ib + 1];
4222/* L106: */
4223 }
4224 ifac[3] = 2;
4225L107:
4226 if (nl != 1) {
4227 goto L104;
4228 }
4229 ifac[1] = *n;
4230 ifac[2] = nf;
4231 tpi = REAL_CONSTANT(6.283185307179586476925286766559005768394338798750211619498891846);
4232 argh = tpi / (real_t) (*n);
4233 is = 0;
4234 nfm1 = nf - 1;
4235 l1 = 1;
4236 if (nfm1 == 0) {
4237 return;
4238 }
4239 i__1 = nfm1;
4240 for (k1 = 1; k1 <= i__1; ++k1) {
4241 ip = ifac[k1 + 2];
4242 ld = 0;
4243 l2 = l1 * ip;
4244 ido = *n / l2;
4245 ipm = ip - 1;
4246 i__2 = ipm;
4247 for (j = 1; j <= i__2; ++j) {
4248 ld += l1;
4249 i__ = is;
4250 argld = (real_t) ld * argh;
4251 fi = REAL_CONSTANT(0.0);
4252 i__3 = ido;
4253 for (ii = 3; ii <= i__3; ii += 2) {
4254 i__ += 2;
4255 fi += REAL_CONSTANT(1.0);
4256 arg = fi * argld;
4257 wa[i__ - 1] = cos(arg);
4258 wa[i__] = sin(arg);
4259/* L108: */
4260 }
4261 is += ido;
4262/* L109: */
4263 }
4264 l1 = l2;
4265/* L110: */
4266 }
4267 return;
4268} /* rffti1_ */
4269
4270/* Subroutine */ void sinqb(integer_t *n, real_t *x, real_t *wsave,
4271 integer_t *ifac)
4272{
4273 /* System generated locals */
4274 integer_t i__1;
4275
4276 /* Local variables */
4277 integer_t k, kc, ns2;
4278 real_t xhold;
4279
4280 /* Parameter adjustments */
4281 --ifac;
4282 --wsave;
4283 --x;
4284
4285 /* Function Body */
4286 if (*n > 1) {
4287 goto L101;
4288 }
4289 x[1] *= REAL_CONSTANT(4.0);
4290 return;
4291L101:
4292 ns2 = *n / 2;
4293 i__1 = *n;
4294 for (k = 2; k <= i__1; k += 2) {
4295 x[k] = -x[k];
4296/* L102: */
4297 }
4298 cosqb(n, &x[1], &wsave[1], &ifac[1]);
4299 i__1 = ns2;
4300 for (k = 1; k <= i__1; ++k) {
4301 kc = *n - k;
4302 xhold = x[k];
4303 x[k] = x[kc + 1];
4304 x[kc + 1] = xhold;
4305/* L103: */
4306 }
4307 return;
4308} /* sinqb_ */
4309
4310/* Subroutine */ void sinqf(integer_t *n, real_t *x, real_t *wsave,
4311 integer_t *ifac)
4312{
4313 /* System generated locals */
4314 integer_t i__1;
4315
4316 /* Local variables */
4317 integer_t k, kc, ns2;
4318 real_t xhold;
4319
4320 /* Parameter adjustments */
4321 --ifac;
4322 --wsave;
4323 --x;
4324
4325 /* Function Body */
4326 if (*n == 1) {
4327 return;
4328 }
4329 ns2 = *n / 2;
4330 i__1 = ns2;
4331 for (k = 1; k <= i__1; ++k) {
4332 kc = *n - k;
4333 xhold = x[k];
4334 x[k] = x[kc + 1];
4335 x[kc + 1] = xhold;
4336/* L101: */
4337 }
4338 cosqf(n, &x[1], &wsave[1], &ifac[1]);
4339 i__1 = *n;
4340 for (k = 2; k <= i__1; k += 2) {
4341 x[k] = -x[k];
4342/* L102: */
4343 }
4344 return;
4345} /* sinqf_ */
4346
4347/* Subroutine */ void sinqi(integer_t *n, real_t *wsave, integer_t *ifac)
4348{
4349 /* Parameter adjustments */
4350 --ifac;
4351 --wsave;
4352
4353 /* Function Body */
4354 cosqi(n, &wsave[1], &ifac[1]);
4355 return;
4356} /* sinqi_ */
4357
4358/* Subroutine */ void sint(integer_t *n, real_t *x, real_t *wsave,
4359 integer_t *ifac)
4360{
4361 integer_t np1, iw1, iw2;
4362
4363 /* Parameter adjustments */
4364 --ifac;
4365 --wsave;
4366 --x;
4367
4368 /* Function Body */
4369 np1 = *n + 1;
4370 iw1 = *n / 2 + 1;
4371 iw2 = iw1 + np1;
4372 sint1(n, &x[1], &wsave[1], &wsave[iw1], &wsave[iw2], &ifac[1]);
4373 return;
4374} /* sint_ */
4375
4376/* Subroutine */ static void sint1(integer_t *n, real_t *war, real_t *was,
4377 real_t *xh, real_t *x, integer_t *ifac)
4378{
4379 /* Initialized data */
4380
4381 static real_t sqrt3 =
4382 REAL_CONSTANT(1.732050807568877293527446341505872366942805253803806280558069795);
4383
4384 /* System generated locals */
4385 integer_t i__1;
4386
4387 /* Local variables */
4388 integer_t i__, k;
4389 real_t t1, t2;
4390 integer_t kc, np1, ns2, modn;
4391 real_t xhold;
4392
4393 /* Parameter adjustments */
4394 --ifac;
4395 --x;
4396 --xh;
4397 --was;
4398 --war;
4399
4400 /* Function Body */
4401 i__1 = *n;
4402 for (i__ = 1; i__ <= i__1; ++i__) {
4403 xh[i__] = war[i__];
4404 war[i__] = x[i__];
4405/* L100: */
4406 }
4407 if ((i__1 = *n - 2) < 0) {
4408 goto L101;
4409 } else if (i__1 == 0) {
4410 goto L102;
4411 } else {
4412 goto L103;
4413 }
4414L101:
4415 xh[1] += xh[1];
4416 goto L106;
4417L102:
4418 xhold = sqrt3 * (xh[1] + xh[2]);
4419 xh[2] = sqrt3 * (xh[1] - xh[2]);
4420 xh[1] = xhold;
4421 goto L106;
4422L103:
4423 np1 = *n + 1;
4424 ns2 = *n / 2;
4425 x[1] = REAL_CONSTANT(0.0);
4426 i__1 = ns2;
4427 for (k = 1; k <= i__1; ++k) {
4428 kc = np1 - k;
4429 t1 = xh[k] - xh[kc];
4430 t2 = was[k] * (xh[k] + xh[kc]);
4431 x[k + 1] = t1 + t2;
4432 x[kc + 1] = t2 - t1;
4433/* L104: */
4434 }
4435 modn = *n % 2;
4436 if (modn != 0) {
4437 x[ns2 + 2] = xh[ns2 + 1] * REAL_CONSTANT(4.0);
4438 }
4439 rfftf1(&np1, &x[1], &xh[1], &war[1], &ifac[1]);
4440 xh[1] = x[1] * REAL_CONSTANT(0.5);
4441 i__1 = *n;
4442 for (i__ = 3; i__ <= i__1; i__ += 2) {
4443 xh[i__ - 1] = -x[i__];
4444 xh[i__] = xh[i__ - 2] + x[i__ - 1];
4445/* L105: */
4446 }
4447 if (modn != 0) {
4448 goto L106;
4449 }
4450 xh[*n] = -x[*n + 1];
4451L106:
4452 i__1 = *n;
4453 for (i__ = 1; i__ <= i__1; ++i__) {
4454 x[i__] = war[i__];
4455 war[i__] = xh[i__];
4456/* L107: */
4457 }
4458 return;
4459} /* sint1_ */
4460
4461/* Subroutine */ void sinti(integer_t *n, real_t *wsave, integer_t *ifac)
4462{
4463 /* Initialized data */
4464
4465 static real_t pi =
4466 REAL_CONSTANT(3.141592653589793238462643383279502884197169399375158209749445923);
4467
4468 /* System generated locals */
4469 integer_t i__1;
4470
4471 /* Local variables */
4472 integer_t k;
4473 real_t dt;
4474 integer_t np1, ns2;
4475
4476 /* Parameter adjustments */
4477 --ifac;
4478 --wsave;
4479
4480 /* Function Body */
4481 if (*n <= 1) {
4482 return;
4483 }
4484 ns2 = *n / 2;
4485 np1 = *n + 1;
4486 dt = pi / (real_t) np1;
4487 i__1 = ns2;
4488 for (k = 1; k <= i__1; ++k) {
4489 wsave[k] = sin(k * dt) * REAL_CONSTANT(2.0);
4490/* L101: */
4491 }
4492 rffti(&np1, &wsave[ns2 + 1], &ifac[1]);
4493 return;
4494} /* sinti_ */
4495
Generated by doxygen 1.13.2