/* fft.f -- translated by f2c (version 20030320). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Subroutine */ int ezffti_(integer *n, real *wsave) { static integer iw1, ns2; extern /* Subroutine */ int rffti_(integer *, real *); /* some minor modifications for ana interfaces, r. shine */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* version 2 february 1978 */ /* a package of fortran subprograms for the fast fourier */ /* transform of periodic and other symmetric sequences */ /* by */ /* paul n swarztrauber */ /* national center for atmospheric research boulder,colorado 80307 */ /* which is sponsored by the national science foundation */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* this package consists of programs which perform fast fourier */ /* transforms for both complex and real periodic sequences and */ /* certain other symmetric sequences that are listed below. */ /* 1. rffti initialize rfftf and rfftb */ /* 2. rfftf forward transform of a real periodic sequence */ /* 3. rfftb backward transform of a real coefficient array */ /* 4. ezfftf a simplified real periodic forward transform */ /* 5. ezfftb a simplified real periodic backward transform */ /* 6. sinti initialize sint */ /* 7. sint sine transform of a real odd sequence */ /* 8. costi initialize cost */ /* 9. cost cosine transform of a real even sequence */ /* 10. sinqi initialize sinqf and sinqb */ /* 11. sinqf forward sine transform with odd wave numbers */ /* 12. sinqb unnormalized inverse of sinqf */ /* 13. cosqi initialize cosqf and cosqb */ /* 14. cosqf forward cosine transform with odd wave numbers */ /* 15. cosqb unnormalized inverse of cosqf */ /* 16. cffti initialize cfftf and cfftb */ /* 17. cfftf forward transform of a complex periodic sequence */ /* 18. cfftb unnormalized inverse of cfftf */ /* --- 8/19/82 this version modified to work only with even n */ /* Parameter adjustments */ --wsave; /* Function Body */ if (*n <= 2) { return 0; } ns2 = *n / 2; /* modn = n-ns2-ns2 */ /* if (modn .ne. 0) go to 101 */ iw1 = *n + 1; rffti_(n, &wsave[iw1]); return 0; /* 101 iw1 = n+n+1 */ /* call cffti (n,wsave(iw1)) */ /* return */ } /* ezffti_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cffti_(integer *n, real *wsave) { static integer iw1, iw2; extern /* Subroutine */ int cffti1_(integer *, real *, integer *); /* subroutine cffti initializes the array wsave which is used in */ /* both cfftf and cfftb. the prime factorization of n together with */ /* a tabulation of the trigonometric functions are computed and */ /* stored in wsave. */ /* input parameter */ /* n the length of the sequence to be transformed */ /* output parameter */ /* wsave a work array which must be dimensioned at least 4*n+15 */ /* the same work array can be used for both cfftf and cfftb */ /* as long as n remains unchanged. different wsave arrays */ /* are required for different values of n. the contents of */ /* wsave must not be changed between calls of cfftf or cfftb. */ /* Parameter adjustments */ --wsave; /* Function Body */ if (*n == 1) { return 0; } iw1 = *n + *n + 1; iw2 = iw1 + *n + *n; cffti1_(n, &wsave[iw1], (integer *) &wsave[iw2]); return 0; } /* cffti_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cffti1_(integer *n, real *wa, integer *ifac) { /* Initialized data */ static integer ntryh[4] = { 3,4,2,5 }; /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer i__, j; static real dc; static integer ib, nf; static real ds; static integer nl, nq, nr, nt; static real tpi, arg1; static integer ntry; /* Parameter adjustments */ --ifac; --wa; /* Function Body */ nl = *n; nf = 0; j = 0; L101: ++j; if (j - 4 <= 0) { goto L102; } else { goto L103; } L102: ntry = ntryh[j - 1]; goto L104; L103: ntry += 2; L104: nq = nl / ntry; nr = nl - ntry * nq; if (nr != 0) { goto L101; } else { goto L105; } L105: ++nf; ifac[nf + 2] = ntry; nl = nq; if (ntry != 2) { goto L107; } if (nf == 1) { goto L107; } i__1 = nf; for (i__ = 2; i__ <= i__1; ++i__) { ib = nf - i__ + 2; ifac[ib + 2] = ifac[ib + 1]; /* L106: */ } ifac[3] = 2; L107: if (nl != 1) { goto L104; } ifac[1] = *n; ifac[2] = nf; tpi = atan(1.f) * 8.f; arg1 = tpi / (real) (*n); dc = cos(arg1); ds = sin(arg1); wa[1] = dc; wa[2] = ds; nt = *n + *n; i__1 = nt; for (i__ = 4; i__ <= i__1; i__ += 2) { wa[i__ - 1] = dc * wa[i__ - 3] - ds * wa[i__ - 2]; wa[i__] = ds * wa[i__ - 3] + dc * wa[i__ - 2]; /* L108: */ } return 0; } /* cffti1_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cfftb_(integer *n, real *c__, real *wsave) { static integer iw1, iw2, iw3; extern /* Subroutine */ int cfftb1_(integer *, real *, real *, real *, integer *); /* subroutine cfftb computes the backward complex discrete fourier */ /* transform (the fourier synthesis). equivalently , cfftb computes */ /* a complex periodic sequence from its fourier coefficients. */ /* the transform is defined below at output parameter c. */ /* a call of cfftf followed by a call of cfftb will multiply the */ /* sequence by n. */ /* the array wsave which is used by subroutine cfftb must be */ /* initialized by calling subroutine cffti(n,wsave). */ /* input parameters */ /* n the length of the complex sequence c. the method is */ /* more efficient when n is the product of small primes. */ /* c a complex array of length n which contains the sequence */ /* wsave a real work array which must be dimensioned at least 4n+15 */ /* in the program that calls cfftb. the wsave array must be */ /* initialized by calling subroutine cffti(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* the same wsave array can be used by cfftf and cfftb. */ /* output parameters */ /* c for j=1,...,n */ /* c(j)=the sum from k=1,...,n of */ /* c(k)*exp(i*j*k*2*pi/n) */ /* where i=sqrt(-1) */ /* wsave contains initialization calculations which must not be */ /* destroyed between calls of subroutine cfftf or cfftb */ /* Parameter adjustments */ --wsave; --c__; /* Function Body */ if (*n == 1) { return 0; } iw1 = *n + *n + 1; iw2 = iw1 + *n + *n; iw3 = iw2 + 1; cfftb1_(n, &c__[1], &wsave[1], &wsave[iw1], (integer *)&wsave[iw2]); return 0; } /* cfftb_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cfftb1_(integer *n, real *c__, real *ch, real *wa, integer *ifac) { /* System generated locals */ integer i__1; /* Local variables */ static integer k1, l1, l2, nf, ip, ix2, ix3, ido, idl1, idot; extern /* Subroutine */ int passb_(integer *, integer *, integer *, integer *, real *, real *, real *, real *, real *, real *), passb2_(integer *, integer *, integer *, real *, real *, real *, real *, real *, real *), passb3_(integer *, integer *, integer *, integer *, real *, real *, real *, real *, real *, real *, real *) , passb4_(integer *, integer *, integer *, integer *, integer *, real *, real *, real *, real *, real *, real *, real *, real *); /* Parameter adjustments */ --ifac; --wa; --ch; --c__; /* Function Body */ nf = ifac[2]; l1 = 1; i__1 = nf; for (k1 = 1; k1 <= i__1; ++k1) { ip = ifac[k1 + 2]; l2 = ip * l1; ido = *n / l2; idot = ido + ido; idl1 = idot * l1; if (ip != 4) { goto L101; } ix2 = l1 + l1; ix3 = ix2 + l1; passb4_(&idot, &l1, &idl1, &ix2, &ix3, &c__[1], &c__[1], &c__[1], &ch[ 1], &ch[1], &wa[1], &wa[1], &wa[1]); goto L104; L101: if (ip != 2) { goto L102; } passb2_(&idot, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[1], &wa[1]); goto L104; L102: if (ip != 3) { goto L103; } ix2 = l1 + l1; passb3_(&idot, &l1, &idl1, &ix2, &c__[1], &c__[1], &c__[1], &ch[1], & ch[1], &wa[1], &wa[1]); goto L104; L103: passb_(&idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[ 1], &wa[1]); L104: l1 = l2; /* L105: */ } return 0; } /* cfftb1_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passb2_(integer *ido, integer *l1, integer *idl1, real * cc, real *c1, real *c2, real *ch, real *ch2, real *wa1) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, i__1, i__2; /* Local variables */ static integer i__, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 3; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L107: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + (c1_dim2 << 1)) * c1_dim1 + 1] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1]; c1[(k + (c1_dim2 << 1)) * c1_dim1 + 2] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2]; /* L108: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L111; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L109: */ } /* L110: */ } return 0; L111: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; } /* passb2_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passb3_(integer *ido, integer *l1, integer *idl1, integer *ix2, real *cc, real *c1, real *c2, real *ch, real *ch2, real *wa1, real *wa2) { /* Initialized data */ static real taur = -.5f; static real taui = .866025403784439f; /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + (cc_dim1 << 2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 1)]; c2[ik + (c2_dim1 << 1)] = ch2[ik + ch2_dim1] + taur * ch2[ik + ( ch2_dim1 << 1)]; c2[ik + c2_dim1 * 3] = taui * ch2[ik + ch2_dim1 * 3]; /* L107: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik - 1 + (ch2_dim1 << 1)] = c2[ik - 1 + (c2_dim1 << 1)] - c2[ik + c2_dim1 * 3]; ch2[ik - 1 + ch2_dim1 * 3] = c2[ik - 1 + (c2_dim1 << 1)] + c2[ik + c2_dim1 * 3]; /* L108: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik + (ch2_dim1 << 1)] = c2[ik + (c2_dim1 << 1)] + c2[ik - 1 + c2_dim1 * 3]; ch2[ik + ch2_dim1 * 3] = c2[ik + (c2_dim1 << 1)] - c2[ik - 1 + c2_dim1 * 3]; /* L109: */ } for (j = 2; j <= 3; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L110: */ } /* L111: */ } if (*ido == 2) { return 0; } if (idot - 1 < *l1) { goto L114; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; L114: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L115: */ } /* L116: */ } return 0; } /* passb3_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passb4_(integer *ido, integer *l1, integer *idl1, integer *ix2, integer *ix3, real *cc, real *c1, real *c2, real *ch, real *ch2, real *wa1, real *wa2, real *wa3) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, wa3_dim1, wa3_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 5; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; wa3_dim1 = *ix3; wa3_offset = 1 + wa3_dim1; wa3 -= wa3_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L106; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 4) * cc_dim1] - cc[i__ + ((k << 2) + 2) * cc_dim1]; /* L101: */ } i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 2) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 4) * cc_dim1]; /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L104: */ } /* L105: */ } goto L111; L106: i__1 = *ido; for (i__ = 2; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 4) * cc_dim1] - cc[i__ + ((k << 2) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 2) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 4) * cc_dim1]; /* L107: */ } /* L108: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L109: */ } /* L110: */ } L111: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + (ch2_dim1 << 1)] + ch2[ik + ch2_dim1 * 3]; /* L112: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + ch2_dim1 * 3] = ch2[ik + (ch2_dim1 << 1)] - ch2[ik + ch2_dim1 * 3]; /* L113: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 1)] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 2)]; /* L114: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 2)] = ch2[ik + ch2_dim1] - ch2[ik + (ch2_dim1 << 2)]; /* L115: */ } for (j = 2; j <= 4; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L116: */ } /* L117: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L120; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L118: */ } /* L119: */ } return 0; L120: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L121: */ } /* L122: */ } return 0; } /* passb4_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passb_(integer *ido, integer *ip, integer *l1, integer * idl1, real *cc, real *c1, real *c2, real *ch, real *ch2, real *wa) { /* System generated locals */ integer ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, i__1, i__2, i__3; /* Local variables */ static integer i__, j, k, l, jc, lc, ik, nt, l1t, idj, idl; static real wai, war; static integer ipp2, idij, idlj, idot, ipph; /* Parameter adjustments */ ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_dim2 = *ip; cc_offset = 1 + cc_dim1 * (1 + cc_dim2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; --wa; /* Function Body */ idot = *ido / 2; nt = *ip * *idl1; ipp2 = *ip + 2; ipph = (*ip + 1) / 2; l1t = *l1 + *l1; if (*ido < *l1) { goto L106; } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *l1; for (k = 1; k <= i__2; ++k) { i__3 = *ido; for (i__ = 1; i__ <= i__3; ++i__) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L101: */ } /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L104: */ } /* L105: */ } goto L112; L106: i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = *l1; for (k = 1; k <= i__3; ++k) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L107: */ } /* L108: */ } /* L109: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L110: */ } /* L111: */ } L112: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L113: */ } idj = 0; i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; idj += *idl1; i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + j * c2_dim1] = ch2[ik + ch2_dim1] + wa[idj - 1] * ch2[ik + (ch2_dim1 << 1)]; c2[ik + jc * c2_dim1] = wa[idj] * ch2[ik + *ip * ch2_dim1]; /* L114: */ } /* L115: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + c2_dim1] += ch2[ik + j * ch2_dim1]; /* L116: */ } /* L117: */ } idl = 0; i__1 = ipph; for (l = 2; l <= i__1; ++l) { lc = ipp2 - l; idl += *idl1; idlj = idl; i__2 = ipph; for (j = 3; j <= i__2; ++j) { jc = ipp2 - j; idlj += idl; if (idlj > nt) { idlj -= nt; } war = wa[idlj - 1]; wai = wa[idlj]; i__3 = *idl1; for (ik = 1; ik <= i__3; ++ik) { c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1]; c2[ik + lc * c2_dim1] += wai * ch2[ik + jc * ch2_dim1]; /* L118: */ } /* L119: */ } /* L120: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik + jc * c2_dim1]; ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik + jc * c2_dim1]; /* L121: */ } /* L122: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc * c2_dim1]; ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc * c2_dim1]; /* L123: */ } /* L124: */ } i__1 = *ip; for (j = 2; j <= i__1; ++j) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L125: */ } /* L126: */ } if (*ido == 2) { return 0; } idj = 0; if (idot > *l1) { goto L130; } i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; idij = 0; i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { idij += idj; i__3 = *l1; for (k = 1; k <= i__3; ++k) { c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L127: */ } /* L128: */ } /* L129: */ } return 0; L130: i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; i__2 = *l1; for (k = 1; k <= i__2; ++k) { idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; /* L131: */ } idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L132: */ } /* L133: */ } /* L134: */ } return 0; } /* passb_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cfftf_(integer *n, real *c__, real *wsave) { static integer iw1, iw2; extern /* Subroutine */ int cfftf1_(integer *, real *, real *, real *, integer *); /* subroutine cfftf computes the forward complex discrete fourier */ /* transform (the fourier analysis). equivalently , cfftf computes */ /* the fourier coefficients of a complex periodic sequence. */ /* the transform is defined below at output parameter c. */ /* the transform is not normalized. to obtain a normalized transform */ /* the output must be divided by n. otherwise a call of cfftf */ /* followed by a call of cfftb will multiply the sequence by n. */ /* the array wsave which is used by subroutine cfftf must be */ /* initialized by calling subroutine cffti(n,wsave). */ /* input parameters */ /* n the length of the complex sequence c. the method is */ /* more efficient when n is the product of small primes. n */ /* c a complex array of length n which contains the sequence */ /* wsave a real work array which must be dimensioned at least 4n+15 */ /* in the program that calls cfftf. the wsave array must be */ /* initialized by calling subroutine cffti(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* the same wsave array can be used by cfftf and cfftb. */ /* output parameters */ /* c for j=1,...,n */ /* c(j)=the sum from k=1,...,n of */ /* c(k)*exp(-i*j*k*2*pi/n) */ /* where i=sqrt(-1) */ /* wsave contains initialization calculations which must not be */ /* destroyed between calls of subroutine cfftf or cfftb */ /* Parameter adjustments */ --wsave; --c__; /* Function Body */ if (*n == 1) { return 0; } iw1 = *n + *n + 1; iw2 = iw1 + *n + *n; cfftf1_(n, &c__[1], &wsave[1], &wsave[iw1], (integer *) &wsave[iw2]); return 0; } /* cfftf_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cfftf1_(integer *n, real *c__, real *ch, real *wa, integer *ifac) { /* System generated locals */ integer i__1; /* Local variables */ static integer k1, l1, l2, nf, ip, ix2, ix3, ido, idl1, idot; extern /* Subroutine */ int passf_(integer *, integer *, integer *, integer *, real *, real *, real *, real *, real *, real *), passf2_(integer *, integer *, integer *, real *, real *, real *, real *, real *, real *), passf3_(integer *, integer *, integer *, integer *, real *, real *, real *, real *, real *, real *, real *) , passf4_(integer *, integer *, integer *, integer *, integer *, real *, real *, real *, real *, real *, real *, real *, real *); /* Parameter adjustments */ --ifac; --wa; --ch; --c__; /* Function Body */ nf = ifac[2]; l1 = 1; i__1 = nf; for (k1 = 1; k1 <= i__1; ++k1) { ip = ifac[k1 + 2]; l2 = ip * l1; ido = *n / l2; idot = ido + ido; idl1 = idot * l1; if (ip != 4) { goto L101; } ix2 = l1 + l1; ix3 = ix2 + l1; passf4_(&idot, &l1, &idl1, &ix2, &ix3, &c__[1], &c__[1], &c__[1], &ch[ 1], &ch[1], &wa[1], &wa[1], &wa[1]); goto L104; L101: if (ip != 2) { goto L102; } passf2_(&idot, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[1], &wa[1]); goto L104; L102: if (ip != 3) { goto L103; } ix2 = l1 + l1; passf3_(&idot, &l1, &idl1, &ix2, &c__[1], &c__[1], &c__[1], &ch[1], & ch[1], &wa[1], &wa[1]); goto L104; L103: passf_(&idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[ 1], &wa[1]); L104: l1 = l2; /* L105: */ } return 0; } /* cfftf1_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passf2_(integer *ido, integer *l1, integer *idl1, real * cc, real *c1, real *c2, real *ch, real *ch2, real *wa1) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, i__1, i__2; /* Local variables */ static integer i__, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 3; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L107: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + (c1_dim2 << 1)) * c1_dim1 + 1] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1]; c1[(k + (c1_dim2 << 1)) * c1_dim1 + 2] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2]; /* L108: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L111; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L109: */ } /* L110: */ } return 0; L111: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; } /* passf2_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passf3_(integer *ido, integer *l1, integer *idl1, integer *ix2, real *cc, real *c1, real *c2, real *ch, real *ch2, real *wa1, real *wa2) { /* Initialized data */ static real taur = -.5f; static real taui = -.866025403784439f; /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + (cc_dim1 << 2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 1)]; c2[ik + (c2_dim1 << 1)] = ch2[ik + ch2_dim1] + taur * ch2[ik + ( ch2_dim1 << 1)]; c2[ik + c2_dim1 * 3] = taui * ch2[ik + ch2_dim1 * 3]; /* L107: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik - 1 + (ch2_dim1 << 1)] = c2[ik - 1 + (c2_dim1 << 1)] - c2[ik + c2_dim1 * 3]; ch2[ik - 1 + ch2_dim1 * 3] = c2[ik - 1 + (c2_dim1 << 1)] + c2[ik + c2_dim1 * 3]; /* L108: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik + (ch2_dim1 << 1)] = c2[ik + (c2_dim1 << 1)] + c2[ik - 1 + c2_dim1 * 3]; ch2[ik + ch2_dim1 * 3] = c2[ik + (c2_dim1 << 1)] - c2[ik - 1 + c2_dim1 * 3]; /* L109: */ } for (j = 2; j <= 3; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L110: */ } /* L111: */ } if (*ido == 2) { return 0; } if (idot - 1 < *l1) { goto L114; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; L114: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L115: */ } /* L116: */ } return 0; } /* passf3_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passf4_(integer *ido, integer *l1, integer *idl1, integer *ix2, integer *ix3, real *cc, real *c1, real *c2, real *ch, real *ch2, real *wa1, real *wa2, real *wa3) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, wa3_dim1, wa3_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 5; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; wa3_dim1 = *ix3; wa3_offset = 1 + wa3_dim1; wa3 -= wa3_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L106; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] - cc[i__ + ((k << 2) + 4) * cc_dim1]; /* L101: */ } i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 4) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 2) * cc_dim1]; /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L104: */ } /* L105: */ } goto L111; L106: i__1 = *ido; for (i__ = 2; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 4) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 2) * cc_dim1]; ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] - cc[i__ + ((k << 2) + 4) * cc_dim1]; /* L107: */ } /* L108: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L109: */ } /* L110: */ } L111: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + (ch2_dim1 << 1)] + ch2[ik + ch2_dim1 * 3]; /* L112: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + ch2_dim1 * 3] = ch2[ik + (ch2_dim1 << 1)] - ch2[ik + ch2_dim1 * 3]; /* L113: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 1)] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 2)]; /* L114: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 2)] = ch2[ik + ch2_dim1] - ch2[ik + (ch2_dim1 << 2)]; /* L115: */ } for (j = 2; j <= 4; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L116: */ } /* L117: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L120; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L118: */ } /* L119: */ } return 0; L120: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L121: */ } /* L122: */ } return 0; } /* passf4_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int passf_(integer *ido, integer *ip, integer *l1, integer * idl1, real *cc, real *c1, real *c2, real *ch, real *ch2, real *wa) { /* System generated locals */ integer ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, i__1, i__2, i__3; /* Local variables */ static integer i__, j, k, l, jc, lc, ik, nt, l1t, idj, idl; static real wai, war; static integer ipp2, idij, idlj, idot, ipph; /* Parameter adjustments */ ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_dim2 = *ip; cc_offset = 1 + cc_dim1 * (1 + cc_dim2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; --wa; /* Function Body */ idot = *ido / 2; nt = *ip * *idl1; ipp2 = *ip + 2; ipph = (*ip + 1) / 2; l1t = *l1 + *l1; if (*ido < *l1) { goto L106; } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *l1; for (k = 1; k <= i__2; ++k) { i__3 = *ido; for (i__ = 1; i__ <= i__3; ++i__) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L101: */ } /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L104: */ } /* L105: */ } goto L112; L106: i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = *l1; for (k = 1; k <= i__3; ++k) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L107: */ } /* L108: */ } /* L109: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L110: */ } /* L111: */ } L112: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L113: */ } idj = 0; i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; idj += *idl1; i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + j * c2_dim1] = ch2[ik + ch2_dim1] + wa[idj - 1] * ch2[ik + (ch2_dim1 << 1)]; c2[ik + jc * c2_dim1] = -wa[idj] * ch2[ik + *ip * ch2_dim1]; /* L114: */ } /* L115: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + c2_dim1] += ch2[ik + j * ch2_dim1]; /* L116: */ } /* L117: */ } idl = 0; i__1 = ipph; for (l = 2; l <= i__1; ++l) { lc = ipp2 - l; idl += *idl1; idlj = idl; i__2 = ipph; for (j = 3; j <= i__2; ++j) { jc = ipp2 - j; idlj += idl; if (idlj > nt) { idlj -= nt; } war = wa[idlj - 1]; wai = wa[idlj]; i__3 = *idl1; for (ik = 1; ik <= i__3; ++ik) { c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1]; c2[ik + lc * c2_dim1] -= wai * ch2[ik + jc * ch2_dim1]; /* L118: */ } /* L119: */ } /* L120: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik + jc * c2_dim1]; ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik + jc * c2_dim1]; /* L121: */ } /* L122: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc * c2_dim1]; ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc * c2_dim1]; /* L123: */ } /* L124: */ } i__1 = *ip; for (j = 2; j <= i__1; ++j) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L125: */ } /* L126: */ } if (*ido == 2) { return 0; } idj = 0; if (idot > *l1) { goto L130; } i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; idij = 0; i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { idij += idj; i__3 = *l1; for (k = 1; k <= i__3; ++k) { c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L127: */ } /* L128: */ } /* L129: */ } return 0; L130: i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; i__2 = *l1; for (k = 1; k <= i__2; ++k) { idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; /* L131: */ } idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L132: */ } /* L133: */ } /* L134: */ } return 0; } /* passf_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cosqi_(integer *n, real *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static real dc, ds, dt; static integer iw1; static real pih, wsk; extern /* Subroutine */ int rffti_(integer *, real *); /* subroutine cosqi initializes the array wsave which is used in */ /* both cosqf and cosqb. the prime factorization of n together with */ /* a tabulation of the trigonometric functions are computed and */ /* stored in wsave. */ /* input parameter */ /* n the length of the sequence to be transformed. n must be */ /* even */ /* output parameter */ /* wsave a work array which must be dimensioned at least 3.5*n+15. */ /* the same work array can be used for both cosqf and cosqb */ /* as long as n remains unchanged. different wsave arrays */ /* are required for different values of n. the contents of */ /* wsave must not be changed between calls of cosqf or cosqb. */ /* Parameter adjustments */ --wsave; /* Function Body */ iw1 = *n + 1; pih = atan(1.f) * 2.f; dt = pih / (real) (*n); dc = cos(dt); ds = sin(dt); wsave[1] = dc; wsk = ds; i__1 = *n; for (k = 2; k <= i__1; ++k) { wsave[k] = dc * wsave[k - 1] - ds * wsk; wsk = ds * wsave[k - 1] + dc * wsk; /* L101: */ } rffti_(n, &wsave[iw1]); return 0; } /* cosqi_ */ /* ------------------------------------------------------------------------ */ /* Subroutine */ int cosqf_(integer *n, real *x, real *wsave) { /* Initialized data */ static real tsq2 = 2.82842712474619f; static real tx; static integer iw1, iw2; static real tsqx; extern /* Subroutine */ int cosqf1_(integer *, real *, real *, real *, real *); /* subroutine cosqf computes the fast fourier transform of quarter */ /* wave data. that is , cosqf computes the coefficients in a cosine */ /* series representation with only odd wave numbers. the transform */ /* is defined below at output parameter x */ /* cosqf is the unnormalized inverse of cosqb since a call of cosqf */ /* followed by a call of cosqb will multiply the input sequence x */ /* by 8*n. */ /* the array wsave which is used by subroutine cosqf must be */ /* initialized by calling subroutine cosqi(n,wsave). */ /* input parameters */ /* n the length of the array x. n must be even and the method */ /* is most efficient when n is a product of small primes. */ /* x an array which contains the sequence to be transformed */ /* wsave a work array which must be dimensioned at least 3.5n+15 */ /* in the program that calls cosqf. the wsave array must be */ /* initialized by calling subroutine cosqi(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* output parameters */ /* x for i=1,...,n */ /* x(i)=2*x(1) plus the sum from k=2 to k=n of */ /* 4.*x(k)*cos((2*i-1)*(k-1)*pi/(2*n)) */ /* a call of cosqf followed by a call of */ /* cosqb will multiply the sequence x by 8*n */ /* therefore cosqb is the unnormalized inverse */ /* of cosqf. */ /* wsave contains initialization calculations which must not */ /* be destroyed between calls of cosqf or cosqb. */ /* Parameter adjustments */ --wsave; --x; /* Function Body */ /* type 2,n, x(1) */ /* L2: */ if (*n > 2) { goto L101; } tx = x[1] + x[1]; tsqx = tsq2 * x[2]; x[1] = tx + tsqx; x[2] = tx - tsqx; return 0; L101: iw1 = *n + 1; iw2 = *n / 2 + iw1; cosqf1_(n, &x[1], &wsave[1], &wsave[iw1], &wsave[iw2]); return 0; } /* cosqf_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cosqf1_(integer *n, real *x, real *w, real *w1, real *xh) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k, kc; static real xn; static integer nm1, np2, ns2; extern /* Subroutine */ int rfftf_(integer *, real *, real *); /* Parameter adjustments */ --xh; --w1; --w; --x; /* Function Body */ ns2 = *n / 2; nm1 = *n - 1; np2 = *n + 2; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; xh[k] = x[k] + x[kc]; xh[kc] = x[k] - x[kc]; /* L101: */ } xh[ns2 + 1] = x[ns2 + 1] + x[ns2 + 1]; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; x[k] = w[k - 1] * xh[kc] + w[kc - 1] * xh[k]; x[kc] = w[k - 1] * xh[k] - w[kc - 1] * xh[kc]; /* L102: */ } x[ns2 + 1] = w[ns2] * xh[ns2 + 1]; rfftf_(n, &x[1], &w1[1]); xn = x[2]; i__1 = nm1; for (i__ = 3; i__ <= i__1; i__ += 2) { x[i__ - 1] = x[i__] + x[i__ + 1]; x[i__] -= x[i__ + 1]; /* L103: */ } x[*n] = xn; return 0; } /* cosqf1_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cosqb_(integer *n, real *x, real *wsave) { /* Initialized data */ static real tsq2 = 2.82842712474619f; static real x1; static integer iw1, iw2; extern /* Subroutine */ int cosqb1_(integer *, real *, real *, real *, real *); /* subroutine cosqb computes the fast fourier transform of quarter */ /* wave data. that is , cosqb computes a sequence from its */ /* representation in terms of a cosine series with odd wave numbers. */ /* the transform is defined below at output parameter x. */ /* cosqb is the unnormalized inverse of cosqf since a call of cosqb */ /* followed by a call of cosqf will multiply the input sequence x */ /* by 8*n. */ /* the array wsave which is used by subroutine cosqb must be */ /* initialized by calling subroutine cosqi(n,wsave). */ /* input parameters */ /* n the length of the array x. n must be even and the method */ /* is most efficient when n is a product of small primes. */ /* x an array which contains the sequence to be transformed */ /* wsave a work array which must be dimensioned at least 3.5n+15 */ /* in the program that calls cosqb. the wsave array must be */ /* initialized by calling subroutine cosqi(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* output parameters */ /* x for i=1,...,n */ /* x(i)= the sum from k=1 to k=n of */ /* 4*x(k)*cos((2*k-1)*(i-1)*pi/(2*n)) */ /* a call of cosqb followed by a call of */ /* cosqf will multiply the sequence x by 8*n */ /* therefore cosqf is the unnormalized inverse */ /* of cosqb. */ /* wsave contains initialization calculations which must not */ /* be destroyed between calls of cosqb or cosqf. */ /* Parameter adjustments */ --wsave; --x; /* Function Body */ if (*n > 2) { goto L101; } x1 = (x[1] + x[2]) * 4.f; x[2] = tsq2 * (x[1] - x[2]); x[1] = x1; return 0; L101: iw1 = *n + 1; iw2 = *n / 2 + iw1; cosqb1_(n, &x[1], &wsave[1], &wsave[iw1], &wsave[iw2]); return 0; } /* cosqb_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cosqb1_(integer *n, real *x, real *w, real *w1, real *xh) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k, kc, ii; static real xn; static integer nm1, np2, ns2; extern /* Subroutine */ int rfftb_(integer *, real *, real *); /* Parameter adjustments */ --xh; --w1; --w; --x; /* Function Body */ ns2 = *n / 2; nm1 = *n - 1; np2 = *n + 2; xn = x[*n]; i__1 = nm1; for (ii = 3; ii <= i__1; ii += 2) { i__ = np2 - ii; x[i__ + 1] = x[i__ - 1] - x[i__]; x[i__] += x[i__ - 1]; /* L101: */ } x[1] += x[1]; x[2] = xn + xn; rfftb_(n, &x[1], &w1[1]); i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; xh[k] = w[k - 1] * x[kc] + w[kc - 1] * x[k]; xh[kc] = w[k - 1] * x[k] - w[kc - 1] * x[kc]; /* L102: */ } x[ns2 + 1] = w[ns2] * (x[ns2 + 1] + x[ns2 + 1]); i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; x[k] = xh[k] + xh[kc]; x[kc] = xh[k] - xh[kc]; /* L103: */ } x[1] += x[1]; return 0; } /* cosqb1_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int costi_(integer *n, real *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static real dt, pi; static integer nm1, iw1, ns2; static real dcs, wck, dss; static integer ns2m; extern /* Subroutine */ int rffti_(integer *, real *); /* subroutine costi initializes the array wsave which is used in */ /* subroutine cost. the prime factorization of n together with */ /* a tabulation of the trigonometric functions are computed and */ /* stored in wsave. */ /* input parameter */ /* n the length of the sequence to be transformed. n must be */ /* odd and greater than 1. */ /* output parameter */ /* wsave a work array which must be dimensioned at least 3*n+15. */ /* different wsave arrays are required for different values */ /* of n. the contents of wsave must not be changed between */ /* calls of cost. */ /* Parameter adjustments */ --wsave; /* Function Body */ if (*n <= 3) { return 0; } nm1 = *n - 1; ns2 = nm1 / 2; ns2m = ns2 - 1; iw1 = ns2 + 1; pi = atan(1.f) * 4.f; dt = pi / (real) nm1; dcs = cos(dt); dss = sin(dt); wsave[1] = dss + dss; wck = dcs + dcs; if (ns2m < 2) { goto L102; } i__1 = ns2m; for (k = 2; k <= i__1; ++k) { wsave[k] = dss * wck + dcs * wsave[k - 1]; wck = dcs * wck - dss * wsave[k - 1]; /* L101: */ } L102: rffti_(&nm1, &wsave[iw1]); return 0; } /* costi_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int cost_(integer *n, real *x, real *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k; static real c1, t1, t2; static integer kc; static real cn; static integer ks, nm1, np1; static real fx2; static integer iw1, ns2; static real tx13, x1px3; extern /* Subroutine */ int rfftf_(integer *, real *, real *); /* subroutine cost computes the discrete fourier cosine transform */ /* of an even sequence x(i). the transform is defined below at output */ /* parameter x. */ /* cost is the unnormalized inverse of itself since a call of cost */ /* followed by another call of cost will multiply the input sequence */ /* x by 8*(n+1). the transform is defined below at output parameter x */ /* the array wsave which is used by subroutine cost must be */ /* initialized by calling subroutine costi(n,wsave). */ /* input parameters */ /* n the length of the sequence x. n must be odd and greater */ /* than 1. the method is most efficient when n-1 is a product */ /* of small primes. */ /* x an array which contains the sequence to be transformed */ /* wsave a work array which must be dimensioned at least 3*n+15 */ /* in the program that calls cost. the wsave array must be */ /* initialized by calling subroutine costi(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* output parameters */ /* x for i=1,...,n */ /* x(i)= 2*x(1)+2*(-1)**(i+1)*x(n) */ /* + the sum from k=2 to k=n-1 */ /* 4*x(k)*cos((k-1)*(i-1)*pi/(n-1)) */ /* a call of cost followed by another call of */ /* cost will multiply the sequence x by 8(n+1) */ /* hence cost is the unnormalized inverse */ /* of itself. */ /* wsave contains initialization calculations which must not be */ /* destroyed between calls of cost. */ /* Parameter adjustments */ --wsave; --x; /* Function Body */ nm1 = *n - 1; np1 = *n + 1; ns2 = nm1 / 2; iw1 = ns2 + 1; if (*n > 3) { goto L101; } x1px3 = x[1] + x[3]; tx13 = x1px3 + x1px3; fx2 = x[2] * 4.f; x[2] = (x[1] - x[3]) * 2.f; x[1] = tx13 + fx2; x[3] = tx13 - fx2; return 0; L101: c1 = x[1] - x[*n]; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = nm1 - k; ks = ns2 - k; c1 += wsave[ks + 1] * (x[k] - x[kc + 2]); /* L102: */ } c1 += c1; x[1] += x[*n]; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np1 - k; t1 = x[k] + x[kc]; t2 = wsave[k - 1] * (x[k] - x[kc]); x[k] = t1 - t2; x[kc] = t1 + t2; /* L103: */ } x[ns2 + 1] += x[ns2 + 1]; rfftf_(&nm1, &x[1], &wsave[iw1]); cn = x[2]; x[2] = c1; i__1 = nm1; for (i__ = 4; i__ <= i__1; i__ += 2) { x[i__] += x[i__ - 2]; /* L104: */ } x[*n] = cn; return 0; } /* cost_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int rffti_(integer *n, real *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static real dc; static integer kc; static real ds, dt; static integer iw1, ns2, nqm; static real tpi; extern /* Subroutine */ int cffti_(integer *, real *); /* subroutine rffti initializes the array wsave which is used in */ /* both rfftf and rfftb. the prime factorization of n together with */ /* a tabulation of the trigonometric functions are computed and */ /* stored in wsave. */ /* input parameter */ /* n the length of the sequence to be transformed. n must be */ /* even */ /* output parameter */ /* wsave a work array which must be dimensioned at least 2.5*n+15. */ /* the same work array can be used for both rfftf and rfftb */ /* as long as n remains unchanged. different wsave arrays */ /* are required for different values of n. the contents of */ /* wsave must not be changed between calls of rfftf or rfftb. */ /* Parameter adjustments */ --wsave; /* Function Body */ ns2 = *n / 2; nqm = (ns2 - 1) / 2; tpi = atan(1.f) * 8.f; dt = tpi / (real) (*n); dc = cos(dt); ds = sin(dt); wsave[1] = dc; wsave[ns2 - 1] = ds; if (nqm < 2) { goto L102; } i__1 = nqm; for (k = 2; k <= i__1; ++k) { kc = ns2 - k; wsave[k] = dc * wsave[k - 1] - ds * wsave[kc + 1]; wsave[kc] = ds * wsave[k - 1] + dc * wsave[kc + 1]; /* L101: */ } L102: iw1 = ns2 + 1; cffti_(&ns2, &wsave[iw1]); return 0; } /* rffti_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int rfftf_(integer *n, real *r__, real *wsave) { static real r1; static integer iw1; extern /* Subroutine */ int rfftf1_(integer *, real *, real *, real *); /* subroutine rfftf computes the fourier coefficients of a real */ /* perodic sequence (fourier analysis). the transform is defined */ /* below at output parameter r. */ /* input parameters */ /* n the length of the array r. n must be even and the method */ /* is most efficient when n is a product of small primes. */ /* n may change so long as different work arrays are provided */ /* r a real array of length n which contains the sequence */ /* to be transformed */ /* wsave a work array which must be dimensioned at least 2.5*n+15 */ /* in the program that calls rfftf. the wsave array must be */ /* initialized by calling subroutine rffti(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* the same wsave array can be used by rfftf and rfftb. */ /* output parameters */ /* r for k=2,...,n/2 */ /* r(2*k-1)= the sum from i=1 to i=n of */ /* 2.*r(i)*cos((k-1)*(i-1)*2*pi/n) */ /* r(2*k)= the sum from i=1 to i=n of */ /* 2.*r(i)*sin((k-1)*(i-1)*2*pi/n) */ /* also */ /* r(1)= the sum from i=1 to i=n of 2.*r(i) */ /* r(2)= the sum from i=1 to i=n of 2.*(-1)**(i+1)*r(i) */ /* ***** note */ /* this transform is unnormalized since a call of rfftf */ /* followed by a call of rfftb will multiply the input */ /* sequence by 2*n. */ /* wsave contains results which must not be destroyed between */ /* calls of rfftf or rfftb. */ /* type *,'n, r(1,1) =', n, r(1,1) */ /* Parameter adjustments */ --wsave; r__ -= 3; /* Function Body */ if (*n > 2) { goto L101; } r1 = (r__[3] + r__[4]) * 2.f; r__[4] = (r__[3] - r__[4]) * 2.f; r__[3] = r1; return 0; L101: iw1 = *n / 2 + 1; rfftf1_(n, &r__[3], &wsave[iw1], &wsave[1]); return 0; } /* rfftf_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int rfftf1_(integer *n, real *x, real *xh, real *w) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, nq, ns2, nqm, nqp, ipar, ns2p2; extern /* Subroutine */ int cfftf_(integer *, real *, real *); static real xhold1, xhold2; /* Parameter adjustments */ --w; xh -= 3; x -= 3; /* Function Body */ ns2 = *n / 2; ns2p2 = ns2 + 2; nq = ns2 / 2; ipar = ns2 - nq - nq; nqm = nq; if (ipar == 0) { --nqm; } nqp = nqm + 1; cfftf_(&ns2, &x[3], &xh[3]); if (nqp < 2) { goto L107; } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = x[(k << 1) + 2] + x[(kc << 1) + 2]; xh[(kc << 1) + 2] = x[(kc << 1) + 1] - x[(k << 1) + 1]; /* L101: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(k << 1) + 1] = x[(k << 1) + 1] + x[(kc << 1) + 1]; xh[(k << 1) + 2] = x[(k << 1) + 2] - x[(kc << 1) + 2]; /* L102: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = w[k - 1] * xh[(kc << 1) + 1] + w[kc - 1] * xh[(kc << 1) + 2]; x[(kc << 1) + 2] = w[k - 1] * xh[(kc << 1) + 2] - w[kc - 1] * xh[(kc << 1) + 1]; /* L103: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = x[(kc << 1) + 1]; xh[(kc << 1) + 2] = x[(kc << 1) + 2]; /* L104: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = xh[(k << 1) + 1] - xh[(kc << 1) + 1]; x[(kc << 1) + 2] = xh[(k << 1) + 2] - xh[(kc << 1) + 2]; /* L105: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(k << 1) + 1] = xh[(k << 1) + 1] + xh[(kc << 1) + 1]; x[(k << 1) + 2] = -xh[(k << 1) + 2] - xh[(kc << 1) + 2]; /* L106: */ } if (ipar != 0) { goto L108; } L107: x[(nqp + 1 << 1) + 1] += x[(nqp + 1 << 1) + 1]; x[(nqp + 1 << 1) + 2] += x[(nqp + 1 << 1) + 2]; L108: xhold1 = x[4] + x[3]; xhold2 = x[3] - x[4]; x[4] = xhold2 + xhold2; x[3] = xhold1 + xhold1; return 0; } /* rfftf1_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int rfftb_(integer *n, real *r__, real *wsave) { static real r1; static integer iw1; extern /* Subroutine */ int rfftb1_(integer *, real *, real *, real *); /* subroutine rfftb computes the real perodic sequence from its */ /* fourier coefficients (fourier synthesis). the transform is defined */ /* below at output parameter r. */ /* input parameters */ /* n the length of the array r. n must be even and the method */ /* is most efficient when n is a product of small primes. */ /* n may change so long as different work arrays are provided */ /* r a real array of length n which contains the sequence */ /* to be transformed */ /* wsave a work array which must be dimensioned at least 2.5*n+15 */ /* in the program that calls rfftb. the wsave array must be */ /* initialized by calling subroutine rffti(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* the same wsave array can be used by rfftf and rfftb. */ /* output parameters */ /* r for i=1,...,n */ /* r(i)=x(1)+(-1)**(i+1)*x(2) */ /* plus the sum from k=2 to k=n/2 of */ /* 2*r(2k-1)*cos((k-1)*(i-1)*2*pi/n) */ /* +2*r(2k)*sin((k-1)*(i-1)*2*pi/n) */ /* ***** note */ /* this transform is unnormalized since a call of rfftf */ /* followed by a call of rfftb will multiply the input */ /* sequence by 2*n. */ /* wsave contains results which must not be destroyed between */ /* calls of rfftb or rfftf. */ /* Parameter adjustments */ --wsave; r__ -= 3; /* Function Body */ if (*n > 2) { goto L101; } r1 = r__[3] + r__[4]; r__[4] = r__[3] - r__[4]; r__[3] = r1; return 0; L101: iw1 = *n / 2 + 1; rfftb1_(n, &r__[3], &wsave[iw1], &wsave[1]); return 0; } /* rfftb_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int rfftb1_(integer *n, real *x, real *xh, real *w) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, nq, ns2, nqm, nqp, ipar, ns2p2; extern /* Subroutine */ int cfftb_(integer *, real *, real *); static real xhold1; /* Parameter adjustments */ --w; xh -= 3; x -= 3; /* Function Body */ ns2 = *n / 2; ns2p2 = ns2 + 2; nq = ns2 / 2; ipar = ns2 - nq - nq; nqm = nq; if (ipar == 0) { --nqm; } nqp = nqm + 1; xhold1 = x[3] - x[4]; x[3] = x[4] + x[3]; x[4] = xhold1; if (ipar != 0) { goto L101; } x[(nqp + 1 << 1) + 1] += x[(nqp + 1 << 1) + 1]; x[(nqp + 1 << 1) + 2] += x[(nqp + 1 << 1) + 2]; L101: if (nqp < 2) { goto L108; } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(k << 1) + 1] = x[(k << 1) + 1] + x[(kc << 1) + 1]; xh[(k << 1) + 2] = x[(kc << 1) + 2] - x[(k << 1) + 2]; /* L102: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = x[(k << 1) + 1] - x[(kc << 1) + 1]; xh[(kc << 1) + 2] = -x[(k << 1) + 2] - x[(kc << 1) + 2]; /* L103: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = xh[(kc << 1) + 1]; x[(kc << 1) + 2] = xh[(kc << 1) + 2]; /* L104: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = w[k - 1] * x[(kc << 1) + 1] - w[kc - 1] * x[(kc << 1) + 2]; xh[(kc << 1) + 2] = w[k - 1] * x[(kc << 1) + 2] + w[kc - 1] * x[(kc << 1) + 1]; /* L105: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(k << 1) + 1] = xh[(k << 1) + 1] - xh[(kc << 1) + 2]; x[(k << 1) + 2] = xh[(k << 1) + 2] + xh[(kc << 1) + 1]; /* L106: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = xh[(k << 1) + 1] + xh[(kc << 1) + 2]; x[(kc << 1) + 2] = xh[(kc << 1) + 1] - xh[(k << 1) + 2]; /* L107: */ } L108: cfftb_(&ns2, &x[3], &xh[3]); return 0; } /* rfftb1_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int sinqi_(integer *n, real *wsave) { extern /* Subroutine */ int cosqi_(integer *, real *); /* subroutine sinqi initializes the array wsave which is used in */ /* both sinqf and sinqb. the prime factorization of n together with */ /* a tabulation of the trigonometric functions are computed and */ /* stored in wsave. */ /* input parameter */ /* n the length of the sequence to be transformed. n must be */ /* even */ /* output parameter */ /* wsave a work array which must be dimensioned at least 3.5*n+15. */ /* the same work array can be used for both sinqf and sinqb */ /* as long as n remains unchanged. different wsave arrays */ /* are required for different values of n. the contents of */ /* wsave must not be changed between calls of sinqf or sinqb. */ /* Parameter adjustments */ --wsave; /* Function Body */ cosqi_(n, &wsave[1]); return 0; } /* sinqi_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int sinqf_(integer *n, real *x, real *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, ns2; extern /* Subroutine */ int cosqf_(integer *, real *, real *); static real xhold; /* subroutine sinqf computes the fast fourier transform of quarter */ /* wave data. that is , sinqf computes the coefficients in a sine */ /* series representation with only odd wave numbers. the transform */ /* is defined below at output parameter x. */ /* sinqb is the unnormalized inverse of sinqf since a call of sinqf */ /* followed by a call of sinqb will multiply the input sequence x */ /* by 8*n. */ /* the array wsave which is used by subroutine sinqf must be */ /* initialized by calling subroutine sinqi(n,wsave). */ /* input parameters */ /* n the length of the array x. n must be even and the method */ /* is most efficient when n is a product of small primes. */ /* x an array which contains the sequence to be transformed */ /* wsave a work array which must be dimensioned at least 3.5n+15 */ /* in the program that calls sinqf. the wsave array must be */ /* initialized by calling subroutine sinqi(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* output parameters */ /* x for i=1,...,n */ /* x(i)=2*(-1)**(i+1)*x(n) */ /* + the sum from k=1 to k=n-1 of */ /* 4*x(k)*sin((2i-1)*k*pi/(2*n)) */ /* a call of sinqf followed by a call of */ /* sinqb will multiply the sequence x by 8*n */ /* therefore sinqb is the unnormalized inverse */ /* of sinqf. */ /* wsave contains initialization calculations which must not */ /* be destroyed between calls of sinqf or sinqb. */ /* Parameter adjustments */ --wsave; --x; /* Function Body */ ns2 = *n / 2; i__1 = ns2; for (k = 1; k <= i__1; ++k) { kc = *n - k; xhold = x[k]; x[k] = x[kc + 1]; x[kc + 1] = xhold; /* L101: */ } cosqf_(n, &x[1], &wsave[1]); i__1 = *n; for (k = 2; k <= i__1; k += 2) { x[k] = -x[k]; /* L102: */ } return 0; } /* sinqf_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int sinqb_(integer *n, real *x, real *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, ns2; extern /* Subroutine */ int cosqb_(integer *, real *, real *); static real xhold; /* subroutine sinqb computes the fast fourier transform of quarter */ /* wave data. that is , sinqb computes a sequence from its */ /* representation in terms of a sine series with odd wave numbers. */ /* the transform is defined below at output parameter x. */ /* sinqf is the unnormalized inverse of sinqb since a call of sinqb */ /* followed by a call of sinqf will multiply the input sequence x */ /* by 8*n. */ /* the array wsave which is used by subroutine sinqb must be */ /* initialized by calling subroutine sinqi(n,wsave). */ /* input parameters */ /* n the length of the array x. n must be even and the method */ /* is most efficient when n is a product of small primes. */ /* x an array which contains the sequence to be transformed */ /* wsave a work array which must be dimensioned at least 3.5n+15 */ /* in the program that calls sinqb. the wsave array must be */ /* initialized by calling subroutine sinqi(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* output parameters */ /* x for i=1,...,n */ /* x(i)= the sum from k=1 to k=n of */ /* 4*x(k)*sin((2k-1)*i*pi/(2*n)) */ /* a call of sinqb followed by a call of */ /* sinqf will multiply the sequence x by 8*n */ /* therefore sinqf is the unnormalized inverse */ /* of sinqb. */ /* wsave contains initialization calculations which must not */ /* be destroyed between calls of sinqb or sinqf. */ /* Parameter adjustments */ --wsave; --x; /* Function Body */ ns2 = *n / 2; i__1 = *n; for (k = 2; k <= i__1; k += 2) { x[k] = -x[k]; /* L101: */ } cosqb_(n, &x[1], &wsave[1]); i__1 = ns2; for (k = 1; k <= i__1; ++k) { kc = *n - k; xhold = x[k]; x[k] = x[kc + 1]; x[kc + 1] = xhold; /* L102: */ } return 0; } /* sinqb_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int sinti_(integer *n, real *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static real dt, pi; static integer np1, iw1, ns2; static real dcs, wck, dss; static integer ns2m; extern /* Subroutine */ int rffti_(integer *, real *); /* subroutine sinti initializes the array wsave which is used in */ /* subroutine sint. the prime factorization of n together with */ /* a tabulation of the trigonometric functions are computed and */ /* stored in wsave. */ /* input parameter */ /* n the length of the sequence to be transformed. n must be */ /* odd. */ /* output parameter */ /* wsave a work array which must be dimensioned at least 3*n+18. */ /* different wsave arrays are required for different values */ /* of n. the contents of wsave must not be changed between */ /* calls of sint. */ /* Parameter adjustments */ --wsave; /* Function Body */ if (*n <= 1) { return 0; } np1 = *n + 1; ns2 = np1 / 2; ns2m = ns2 - 1; iw1 = ns2 + 1; pi = atan(1.f) * 4.f; dt = pi / (real) np1; dcs = cos(dt); dss = sin(dt); wsave[1] = dss + dss; wck = dcs + dcs; if (ns2m < 2) { goto L102; } i__1 = ns2m; for (k = 2; k <= i__1; ++k) { wsave[k] = dss * wck + dcs * wsave[k - 1]; wck = dcs * wck - dss * wsave[k - 1]; /* L101: */ } L102: rffti_(&np1, &wsave[iw1]); return 0; } /* sinti_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int sint_(integer *n, real *x, real *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k; static real t1, t2, x1; static integer kc, np1, iw1, ns2, ns2m; extern /* Subroutine */ int rfftf_(integer *, real *, real *); /* subroutine sint computes the discrete fourier sine transform */ /* of an odd sequence x(i). the transform is defined below at */ /* output parameter x. */ /* sint is the unnormalized inverse of itself since a call of sint */ /* followed by another call of sint will multiply the input sequence */ /* x by 8*(n+1). */ /* the array wsave which is used by subroutine sint must be */ /* initialized by calling subroutine sinti(n,wsave). */ /* input parameters */ /* n n is the length of x. n must be odd and the method is */ /* most efficient when n+1 is a product of small primes */ /* x an array which contains the sequence to be transformed */ /* ************important************* */ /* x must be dimensioned at least n+1 */ /* wsave a work array which must be dimensioned at least 3n+18 */ /* in the program that calls sint. the wsave array must be */ /* initialized by calling subroutine sinti(n,wsave) and a */ /* different wsave array must be used for each different */ /* value of n. this initialization does not have to be */ /* repeated so long as n remains unchanged thus subsequent */ /* transforms can be obtained faster than the first. */ /* output parameters */ /* x for i=1,...,n */ /* x(i)= the sum from k=1 to k=n */ /* 4*x(k)*sin(k*i*pi/(n+1)) */ /* a call of sint followed by another call of */ /* sint will multiply the sequence x by 8(n+1) */ /* hence sint is the unnormalized inverse */ /* of itself. */ /* wsave contains initialization calculations which must not be */ /* destroyed between calls of sint. */ /* Parameter adjustments */ --wsave; --x; /* Function Body */ if (*n > 1) { goto L101; } x[1] *= 4.f; return 0; L101: np1 = *n + 1; ns2 = np1 / 2; ns2m = ns2 - 1; iw1 = ns2 + 1; x1 = x[1]; x[1] = 0.f; i__1 = ns2m; for (k = 1; k <= i__1; ++k) { kc = np1 - k; t1 = x1 - x[kc]; t2 = wsave[k] * (x1 + x[kc]); x1 = x[k + 1]; x[k + 1] = t1 + t2; x[kc + 1] = t2 - t1; /* L102: */ } x[ns2 + 1] = x1 * 4.f; rfftf_(&np1, &x[1], &wsave[iw1]); x[1] *= .5f; i__1 = *n; for (i__ = 3; i__ <= i__1; i__ += 2) { x[i__ - 1] = x[i__ + 1]; x[i__] += x[i__ - 2]; /* L103: */ } return 0; } /* sint_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezpower_(integer *n, real *r__, real *p, real *wsave) { /* System generated locals */ integer i__1; real r__1, r__2; /* Local variables */ static integer i__; static real cf, cf2; static integer iw1, ns2, ns2m; extern /* Subroutine */ int rfftf_(integer *, real *, real *); /* modified from ezfftf to just return power */ /* Parameter adjustments */ --wsave; --p; --r__; /* Function Body */ if (*n <= 2) { /* type 2,n */ /* L2: */ return 0; } /* L102: */ ns2 = *n / 2; /* modn = n-ns2-ns2 */ /* if (modn .ne. 0) go to 105 */ iw1 = *n + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { wsave[i__] = r__[i__]; /* L103: */ } rfftf_(n, &wsave[1], &wsave[iw1]); cf = 1.f / (real) (*n); cf2 = cf * cf; /* azero = .5*cf*wsave(1) */ p[1] = cf2 * .25f * wsave[1] * wsave[1]; /* a(ns2) = .5*cf*wsave(2) */ /* b(ns2) = 0. */ p[ns2 + 1] = cf2 * .25f * wsave[2] * wsave[2]; ns2m = ns2 - 1; i__1 = ns2m; for (i__ = 1; i__ <= i__1; ++i__) { /* a(i) = cf*wsave(2*i+1) */ /* b(i) = cf*wsave(2*i+2) */ /* Computing 2nd power */ r__1 = wsave[(i__ << 1) + 2]; /* Computing 2nd power */ r__2 = wsave[(i__ << 1) + 1]; p[i__ + 1] = cf2 * (r__1 * r__1 + r__2 * r__2); /* L104: */ } return 0; /* 105 iw1 = n+n+1 */ /* do 106 i=1,n */ /* wsave(2*i-1) = r(i) */ /* wsave(2*i) = 0. */ /* 106 continue */ /* call cfftf (n,wsave,wsave(iw1)) */ /* cf = 2./float(n) */ /* cf2=cf*cf */ /* c azero = .5*cf*wsave(1) */ /* nm1s2 = (n-1)/2 */ /* do 107 i=1,nm1s2 */ /* c a(i) = cf*wsave(2*i+1) */ /* c b(i) = -cf*wsave(2*i+2) */ /* p(i)=cf2*(wsave(2*i+2)**2+wsave(2*i+1)**2) */ /* 107 continue */ /* return */ } /* ezpower_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezsc_(integer *n, real *r__, real *b, real *a, real * wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__; static real cf; static integer iw1, ns2; extern /* Subroutine */ int rfftf_(integer *, real *, real *); /* -- modified from ezfftf to return sine and cosine series for a real */ /* - variable in a format suitable for ana, this means that first */ /* cosine is set to mean (azero here) and first and last sine are 0 */ /* the sine corresponds to b and the cosine to a */ /* we also assume that n is always even (checked by calling routine */ /* which is sc in fun2.mar */ /* - also note that call to ezffti to set up wsave is made in fun2 */ /* Parameter adjustments */ --wsave; --a; --b; --r__; /* Function Body */ if (*n <= 2) { /* type 2,n */ /* L2: */ return 0; } /* L102: */ ns2 = *n / 2; iw1 = *n + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { wsave[i__] = r__[i__]; /* L103: */ } rfftf_(n, &wsave[1], &wsave[iw1]); cf = 1.f / (real) (*n); a[1] = cf * .5f * wsave[1]; /* this was azero */ b[1] = 0.f; /* first and last sine terms are 0 */ a[ns2 + 1] = cf * .5f * wsave[2]; b[ns2 + 1] = 0.f; i__1 = ns2; for (i__ = 2; i__ <= i__1; ++i__) { a[i__] = cf * wsave[(i__ << 1) - 1]; b[i__] = cf * wsave[i__ * 2]; /* L104: */ } return 0; } /* ezsc_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezdsc_(integer *n, real *r__, real *b, real *a, real * wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, iw1, ns2; extern /* Subroutine */ int rfftf_(integer *, real *, real *); /* -- modified from ezfftf to return sine and cosine series for a real */ /* - variable, this is different and more sensible than the sc */ /* - routine for multi-dimensional extensions */ /* - the lengths are each n/2 and the first value in the sine */ /* - array is actually the highest cosine coefficient */ /* we also assume that n is always even (checked by calling routine */ /* which is dsc in fun2.mar */ /* - result is unnormalized to save a little time, divide by */ /* - n to convert most to the actual sine/cosine coefficients */ /* - also note that call to ezffti to set up wsave is made in fun2 */ /* Parameter adjustments */ --wsave; --a; --b; --r__; /* Function Body */ if (*n <= 2) { /* type 2,n */ /* L2: */ return 0; } /* L102: */ ns2 = *n / 2; iw1 = *n + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { wsave[i__] = r__[i__]; /* L103: */ } rfftf_(n, &wsave[1], &wsave[iw1]); i__1 = ns2; for (i__ = 1; i__ <= i__1; ++i__) { a[i__] = wsave[(i__ << 1) - 1]; b[i__] = wsave[i__ * 2]; /* L104: */ } return 0; } /* ezdsc_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezscb_(integer *n, real *r__, real *a, real *b, real * wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, iw1, ns2; extern /* Subroutine */ int rfftb_(integer *, real *, real *); /* -- this is a modified version of ezfftb for ana */ /* it is similar but the sine and cosine series are 1 element longer */ /* with c(1)=azero and s(1)=0 always assumed */ /* also we assume n is even */ /* Parameter adjustments */ --r__; --a; --b; --wsave; /* Function Body */ if (*n > 1) { goto L101; } r__[1] = a[1]; return 0; L101: if (*n > 2) { goto L102; } r__[1] = a[1] + a[2]; r__[2] = a[1] - a[2]; return 0; L102: ns2 = *n / 2; iw1 = *n + 1; i__1 = ns2; for (i__ = 2; i__ <= i__1; ++i__) { r__[(i__ << 1) - 1] = a[i__] * .5f; r__[i__ * 2] = b[i__] * .5f; /* L103: */ } r__[1] = a[1]; r__[2] = a[ns2 + 1]; rfftb_(n, &r__[1], &wsave[iw1]); return 0; } /* ezscb_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezfftb_(integer *n, real *r__, real *azero, real *a, real *b, real *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, ic, iw1, ns2, ns2m, modn, nm1s2; extern /* Subroutine */ int cfftb_(integer *, real *, real *), rfftb_( integer *, real *, real *), ezffti_(integer *, real *); /* subroutine ezfftb computes a real perodic sequence from its */ /* fourier coefficients (fourier synthesis). the transform is */ /* defined below at output parameter r. ezfftb is a simplified */ /* version of rfftb. it is not as fast as rfftb since scaling and */ /* initialization are computed for each transform. the repeated */ /* initialization can be suppressed by removeing the statment */ /* ( call ezffti(n,wsave) ) from both ezfftf and ezfftb and inserting */ /* it at the appropriate place in your program. */ /* input parameters */ /* n the length of the array r. ezfftb is about twice as fast */ /* for even n as it is for odd n. also ezfftb is more */ /* efficient when n is a product of small primes. */ /* azero the constant fourier coefficient */ /* a,b arrays which contain the remaining fourier coefficients */ /* these arrays are not destroyed. */ /* the length of these arrays depends on whether n is even or */ /* odd. */ /* if n is even n/2 locations are required */ /* if n is odd (n-1)/2 locations are required */ /* wsave a work array whose length depends on whether n is even or */ /* odd. */ /* if n is even 3.5*n+15 locations are required */ /* if n is odd 6*n+15 locations are required */ /* output parameters */ /* r if n is even define kmax=n/2 */ /* if n is odd define kmax=(n-1)/2 */ /* then for i=1,...,n */ /* r(i)=azero plus the sum from k=1 to k=kmax of */ /* a(k)*cos(k*(i-1)*2*pi/n)+b(k)*sin(k*(i-1)*2*pi/n) */ /* ********************* complex notation ************************** */ /* for j=1,...,n */ /* r(j) equals the sum from k=-kmax to k=kmax of */ /* c(k)*exp(i*k*(j-1)*2*pi/n) */ /* where */ /* c(k) = .5*cmplx(a(k),-b(k)) for k=1,...,kmax */ /* c(-k) = conjg(c(k)) */ /* c(0) = azero */ /* and i=sqrt(-1) */ /* *************** amplitude - phase notation *********************** */ /* for i=1,...,n */ /* r(i) equals azero plus the sum from k=1 to k=kmax of */ /* alpha(k)*cos(k*(i-1)*2*pi/n+beta(k)) */ /* where */ /* alpha(k) = sqrt(a(k)*a(k)+b(k)*b(k)) */ /* cos(beta(k))=a(k)/alpha(k) */ /* sin(beta(k))=-b(k)/alpha(k) */ /* Parameter adjustments */ --wsave; --b; --a; --r__; /* Function Body */ if (*n > 1) { goto L101; } r__[1] = *azero; return 0; L101: if (*n > 2) { goto L102; } r__[1] = *azero + a[1]; r__[2] = *azero - a[1]; return 0; L102: ns2 = *n / 2; /* to supress repeated initialization remove the following statment */ /* ( call ezffti(n,wsave)) from both this program and ezfftf and */ /* insert it at the appropriate place in your program. */ ezffti_(n, &wsave[1]); modn = *n - ns2 - ns2; if (modn != 0) { goto L104; } iw1 = *n + 1; ns2m = ns2 - 1; i__1 = ns2m; for (i__ = 1; i__ <= i__1; ++i__) { r__[(i__ << 1) + 1] = a[i__] * .5f; r__[(i__ << 1) + 2] = b[i__] * .5f; /* L103: */ } r__[1] = *azero; r__[2] = a[ns2]; rfftb_(n, &r__[1], &wsave[iw1]); return 0; L104: iw1 = *n + *n + 1; nm1s2 = (*n - 1) / 2; i__1 = nm1s2; for (i__ = 1; i__ <= i__1; ++i__) { ic = *n - i__; wsave[(i__ << 1) + 1] = a[i__]; wsave[(i__ << 1) + 2] = -b[i__]; wsave[(ic << 1) + 1] = a[i__]; wsave[(ic << 1) + 2] = b[i__]; /* L105: */ } wsave[1] = *azero + *azero; wsave[2] = 0.f; cfftb_(n, &wsave[1], &wsave[iw1]); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { r__[i__] = wsave[(i__ << 1) - 1] * .5f; /* L106: */ } return 0; } /* ezfftb_ */ /* ----------------------------------------------------------------------------- */ /* ----------------------------------------------------------------------------- */ /* the dp versions follow */ /* Subroutine */ int ezfftid_(integer *n, doublereal *wsave) { static integer iw1, ns2; extern /* Subroutine */ int rfftid_(integer *, doublereal *); /* --- 8/19/82 this version modified to work only with even n */ /* Parameter adjustments */ --wsave; /* Function Body */ if (*n <= 2) { return 0; } ns2 = *n / 2; iw1 = *n + 1; rfftid_(n, &wsave[iw1]); return 0; } /* ezfftid_ */ /* Subroutine */ int cfftid_(integer *n, doublereal *wsave) { static integer iw1, iw2; extern /* Subroutine */ int cffti1d_(integer *, doublereal *, integer * ); /* Parameter adjustments */ --wsave; /* Function Body */ if (*n == 1) { return 0; } iw1 = *n + *n + 1; iw2 = iw1 + *n + *n; cffti1d_(n, &wsave[iw1], (integer *) &wsave[iw2]); return 0; } /* cfftid_ */ /* Subroutine */ int cffti1d_(integer *n, doublereal *wa, integer *ifac) { /* Initialized data */ static integer ntryh[4] = { 3,4,2,5 }; /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer i__, j; static doublereal dc; static integer ib, nf; static doublereal ds; static integer nl, nq, nr, nt; static doublereal tpi, arg1; static integer ntry; /* Parameter adjustments */ --ifac; --wa; /* Function Body */ nl = *n; nf = 0; j = 0; L101: ++j; if (j - 4 <= 0) { goto L102; } else { goto L103; } L102: ntry = ntryh[j - 1]; goto L104; L103: ntry += 2; L104: nq = nl / ntry; nr = nl - ntry * nq; if (nr != 0) { goto L101; } else { goto L105; } L105: ++nf; ifac[nf + 2] = ntry; nl = nq; if (ntry != 2) { goto L107; } if (nf == 1) { goto L107; } i__1 = nf; for (i__ = 2; i__ <= i__1; ++i__) { ib = nf - i__ + 2; ifac[ib + 2] = ifac[ib + 1]; /* L106: */ } ifac[3] = 2; L107: if (nl != 1) { goto L104; } ifac[1] = *n; ifac[2] = nf; tpi = atan(1.f) * 8.f; arg1 = tpi / (real) (*n); dc = cos(arg1); ds = sin(arg1); wa[1] = dc; wa[2] = ds; nt = *n + *n; i__1 = nt; for (i__ = 4; i__ <= i__1; i__ += 2) { wa[i__ - 1] = dc * wa[i__ - 3] - ds * wa[i__ - 2]; wa[i__] = ds * wa[i__ - 3] + dc * wa[i__ - 2]; /* L108: */ } return 0; } /* cffti1d_ */ /* Subroutine */ int cfftbd_(integer *n, doublereal *c__, doublereal *wsave) { static integer iw1, iw2, iw3; extern /* Subroutine */ int cfftb1d_(integer *, doublereal *, doublereal * , doublereal *, integer *); /* Parameter adjustments */ --wsave; --c__; /* Function Body */ if (*n == 1) { return 0; } iw1 = *n + *n + 1; iw2 = iw1 + *n + *n; iw3 = iw2 + 1; cfftb1d_(n, &c__[1], &wsave[1], &wsave[iw1], (integer *) &wsave[iw2]); return 0; } /* cfftbd_ */ /* Subroutine */ int cfftb1d_(integer *n, doublereal *c__, doublereal *ch, doublereal *wa, integer *ifac) { /* System generated locals */ integer i__1; /* Local variables */ static integer k1, l1, l2, nf, ip, ix2, ix3, ido, idl1, idot; extern /* Subroutine */ int passbd_(integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), passb2d_(integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), passb3d_(integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), passb4d_( integer *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); /* Parameter adjustments */ --ifac; --wa; --ch; --c__; /* Function Body */ nf = ifac[2]; l1 = 1; i__1 = nf; for (k1 = 1; k1 <= i__1; ++k1) { ip = ifac[k1 + 2]; l2 = ip * l1; ido = *n / l2; idot = ido + ido; idl1 = idot * l1; if (ip != 4) { goto L101; } ix2 = l1 + l1; ix3 = ix2 + l1; passb4d_(&idot, &l1, &idl1, &ix2, &ix3, &c__[1], &c__[1], &c__[1], & ch[1], &ch[1], &wa[1], &wa[1], &wa[1]); goto L104; L101: if (ip != 2) { goto L102; } passb2d_(&idot, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[1], &wa[1]); goto L104; L102: if (ip != 3) { goto L103; } ix2 = l1 + l1; passb3d_(&idot, &l1, &idl1, &ix2, &c__[1], &c__[1], &c__[1], &ch[1], & ch[1], &wa[1], &wa[1]); goto L104; L103: passbd_(&idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], & ch[1], &wa[1]); L104: l1 = l2; /* L105: */ } return 0; } /* cfftb1d_ */ /* Subroutine */ int passb2d_(integer *ido, integer *l1, integer *idl1, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa1) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, i__1, i__2; /* Local variables */ static integer i__, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 3; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L107: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + (c1_dim2 << 1)) * c1_dim1 + 1] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1]; c1[(k + (c1_dim2 << 1)) * c1_dim1 + 2] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2]; /* L108: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L111; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L109: */ } /* L110: */ } return 0; L111: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; } /* passb2d_ */ /* Subroutine */ int passb3d_(integer *ido, integer *l1, integer *idl1, integer *ix2, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa1, doublereal *wa2) { /* Initialized data */ static doublereal taur = -.5; static doublereal taui = .866025403784439; /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + (cc_dim1 << 2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 1)]; c2[ik + (c2_dim1 << 1)] = ch2[ik + ch2_dim1] + taur * ch2[ik + ( ch2_dim1 << 1)]; c2[ik + c2_dim1 * 3] = taui * ch2[ik + ch2_dim1 * 3]; /* L107: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik - 1 + (ch2_dim1 << 1)] = c2[ik - 1 + (c2_dim1 << 1)] - c2[ik + c2_dim1 * 3]; ch2[ik - 1 + ch2_dim1 * 3] = c2[ik - 1 + (c2_dim1 << 1)] + c2[ik + c2_dim1 * 3]; /* L108: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik + (ch2_dim1 << 1)] = c2[ik + (c2_dim1 << 1)] + c2[ik - 1 + c2_dim1 * 3]; ch2[ik + ch2_dim1 * 3] = c2[ik + (c2_dim1 << 1)] - c2[ik - 1 + c2_dim1 * 3]; /* L109: */ } for (j = 2; j <= 3; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L110: */ } /* L111: */ } if (*ido == 2) { return 0; } if (idot - 1 < *l1) { goto L114; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; L114: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L115: */ } /* L116: */ } return 0; } /* passb3d_ */ /* Subroutine */ int passb4d_(integer *ido, integer *l1, integer *idl1, integer *ix2, integer *ix3, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa1, doublereal *wa2, doublereal *wa3) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, wa3_dim1, wa3_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 5; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; wa3_dim1 = *ix3; wa3_offset = 1 + wa3_dim1; wa3 -= wa3_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L106; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 4) * cc_dim1] - cc[i__ + ((k << 2) + 2) * cc_dim1]; /* L101: */ } i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 2) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 4) * cc_dim1]; /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L104: */ } /* L105: */ } goto L111; L106: i__1 = *ido; for (i__ = 2; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 4) * cc_dim1] - cc[i__ + ((k << 2) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 2) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 4) * cc_dim1]; /* L107: */ } /* L108: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L109: */ } /* L110: */ } L111: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + (ch2_dim1 << 1)] + ch2[ik + ch2_dim1 * 3]; /* L112: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + ch2_dim1 * 3] = ch2[ik + (ch2_dim1 << 1)] - ch2[ik + ch2_dim1 * 3]; /* L113: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 1)] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 2)]; /* L114: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 2)] = ch2[ik + ch2_dim1] - ch2[ik + (ch2_dim1 << 2)]; /* L115: */ } for (j = 2; j <= 4; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L116: */ } /* L117: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L120; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L118: */ } /* L119: */ } return 0; L120: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] + wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] - wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L121: */ } /* L122: */ } return 0; } /* passb4d_ */ /* Subroutine */ int passbd_(integer *ido, integer *ip, integer *l1, integer * idl1, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa) { /* System generated locals */ integer ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, i__1, i__2, i__3; /* Local variables */ static integer i__, j, k, l, jc, lc, ik, nt, l1t, idj, idl; static doublereal wai, war; static integer ipp2, idij, idlj, idot, ipph; /* Parameter adjustments */ ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_dim2 = *ip; cc_offset = 1 + cc_dim1 * (1 + cc_dim2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; --wa; /* Function Body */ idot = *ido / 2; nt = *ip * *idl1; ipp2 = *ip + 2; ipph = (*ip + 1) / 2; l1t = *l1 + *l1; if (*ido < *l1) { goto L106; } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *l1; for (k = 1; k <= i__2; ++k) { i__3 = *ido; for (i__ = 1; i__ <= i__3; ++i__) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L101: */ } /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L104: */ } /* L105: */ } goto L112; L106: i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = *l1; for (k = 1; k <= i__3; ++k) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L107: */ } /* L108: */ } /* L109: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L110: */ } /* L111: */ } L112: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L113: */ } idj = 0; i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; idj += *idl1; i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + j * c2_dim1] = ch2[ik + ch2_dim1] + wa[idj - 1] * ch2[ik + (ch2_dim1 << 1)]; c2[ik + jc * c2_dim1] = wa[idj] * ch2[ik + *ip * ch2_dim1]; /* L114: */ } /* L115: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + c2_dim1] += ch2[ik + j * ch2_dim1]; /* L116: */ } /* L117: */ } idl = 0; i__1 = ipph; for (l = 2; l <= i__1; ++l) { lc = ipp2 - l; idl += *idl1; idlj = idl; i__2 = ipph; for (j = 3; j <= i__2; ++j) { jc = ipp2 - j; idlj += idl; if (idlj > nt) { idlj -= nt; } war = wa[idlj - 1]; wai = wa[idlj]; i__3 = *idl1; for (ik = 1; ik <= i__3; ++ik) { c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1]; c2[ik + lc * c2_dim1] += wai * ch2[ik + jc * ch2_dim1]; /* L118: */ } /* L119: */ } /* L120: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik + jc * c2_dim1]; ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik + jc * c2_dim1]; /* L121: */ } /* L122: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc * c2_dim1]; ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc * c2_dim1]; /* L123: */ } /* L124: */ } i__1 = *ip; for (j = 2; j <= i__1; ++j) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L125: */ } /* L126: */ } if (*ido == 2) { return 0; } idj = 0; if (idot > *l1) { goto L130; } i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; idij = 0; i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { idij += idj; i__3 = *l1; for (k = 1; k <= i__3; ++k) { c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L127: */ } /* L128: */ } /* L129: */ } return 0; L130: i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; i__2 = *l1; for (k = 1; k <= i__2; ++k) { idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; /* L131: */ } idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L132: */ } /* L133: */ } /* L134: */ } return 0; } /* passbd_ */ /* Subroutine */ int cfftfd_(integer *n, doublereal *c__, doublereal *wsave) { static integer iw1, iw2; extern /* Subroutine */ int cfftf1d_(integer *, doublereal *, doublereal * , doublereal *, integer *); /* Parameter adjustments */ --wsave; --c__; /* Function Body */ if (*n == 1) { return 0; } iw1 = *n + *n + 1; iw2 = iw1 + *n + *n; cfftf1d_(n, &c__[1], &wsave[1], &wsave[iw1], (integer *) &wsave[iw2]); return 0; } /* cfftfd_ */ /* Subroutine */ int cfftf1d_(integer *n, doublereal *c__, doublereal *ch, doublereal *wa, integer *ifac) { /* System generated locals */ integer i__1; /* Local variables */ static integer k1, l1, l2, nf, ip, ix2, ix3, ido, idl1, idot; extern /* Subroutine */ int passfd_(integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), passf2d_(integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), passf3d_(integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), passf4d_( integer *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *); /* Parameter adjustments */ --ifac; --wa; --ch; --c__; /* Function Body */ nf = ifac[2]; l1 = 1; i__1 = nf; for (k1 = 1; k1 <= i__1; ++k1) { ip = ifac[k1 + 2]; l2 = ip * l1; ido = *n / l2; idot = ido + ido; idl1 = idot * l1; if (ip != 4) { goto L101; } ix2 = l1 + l1; ix3 = ix2 + l1; passf4d_(&idot, &l1, &idl1, &ix2, &ix3, &c__[1], &c__[1], &c__[1], & ch[1], &ch[1], &wa[1], &wa[1], &wa[1]); goto L104; L101: if (ip != 2) { goto L102; } passf2d_(&idot, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[1], &wa[1]); goto L104; L102: if (ip != 3) { goto L103; } ix2 = l1 + l1; passf3d_(&idot, &l1, &idl1, &ix2, &c__[1], &c__[1], &c__[1], &ch[1], & ch[1], &wa[1], &wa[1]); goto L104; L103: passfd_(&idot, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], & ch[1], &wa[1]); L104: l1 = l2; /* L105: */ } return 0; } /* cfftf1d_ */ /* Subroutine */ int passf2d_(integer *ido, integer *l1, integer *idl1, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa1) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, i__1, i__2; /* Local variables */ static integer i__, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 3; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[i__ + ((k << 1) + 2) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 1) + 1) * cc_dim1] - cc[i__ + ((k << 1) + 2) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L107: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + (c1_dim2 << 1)) * c1_dim1 + 1] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1]; c1[(k + (c1_dim2 << 1)) * c1_dim1 + 2] = ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2]; /* L108: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L111; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L109: */ } /* L110: */ } return 0; L111: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; } /* passf2d_ */ /* Subroutine */ int passf3d_(integer *ido, integer *l1, integer *idl1, integer *ix2, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa1, doublereal *wa2) { /* Initialized data */ static doublereal taur = -.5; static doublereal taui = -.866025403784439; /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + (cc_dim1 << 2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L103; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L101: */ } /* L102: */ } goto L106; L103: i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * 3 + 1) * cc_dim1]; ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] + cc[i__ + (k * 3 + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + (k * 3 + 2) * cc_dim1] - cc[i__ + (k * 3 + 3) * cc_dim1]; /* L104: */ } /* L105: */ } L106: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 1)]; c2[ik + (c2_dim1 << 1)] = ch2[ik + ch2_dim1] + taur * ch2[ik + ( ch2_dim1 << 1)]; c2[ik + c2_dim1 * 3] = taui * ch2[ik + ch2_dim1 * 3]; /* L107: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik - 1 + (ch2_dim1 << 1)] = c2[ik - 1 + (c2_dim1 << 1)] - c2[ik + c2_dim1 * 3]; ch2[ik - 1 + ch2_dim1 * 3] = c2[ik - 1 + (c2_dim1 << 1)] + c2[ik + c2_dim1 * 3]; /* L108: */ } i__1 = *idl1; for (ik = 2; ik <= i__1; ik += 2) { ch2[ik + (ch2_dim1 << 1)] = c2[ik + (c2_dim1 << 1)] + c2[ik - 1 + c2_dim1 * 3]; ch2[ik + ch2_dim1 * 3] = c2[ik + (c2_dim1 << 1)] - c2[ik - 1 + c2_dim1 * 3]; /* L109: */ } for (j = 2; j <= 3; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L110: */ } /* L111: */ } if (*ido == 2) { return 0; } if (idot - 1 < *l1) { goto L114; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L112: */ } /* L113: */ } return 0; L114: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; /* L115: */ } /* L116: */ } return 0; } /* passf3d_ */ /* Subroutine */ int passf4d_(integer *ido, integer *l1, integer *idl1, integer *ix2, integer *ix3, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa1, doublereal *wa2, doublereal *wa3) { /* System generated locals */ integer cc_dim1, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch_dim1, ch_dim2, ch_offset, ch2_dim1, ch2_offset, wa1_dim1, wa1_offset, wa2_dim1, wa2_offset, wa3_dim1, wa3_offset, i__1, i__2; /* Local variables */ static integer i__, j, k, ik, idot; /* Parameter adjustments */ wa1_dim1 = *l1; wa1_offset = 1 + wa1_dim1; wa1 -= wa1_offset; ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_offset = 1 + cc_dim1 * 5; cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; wa2_dim1 = *ix2; wa2_offset = 1 + wa2_dim1; wa2 -= wa2_offset; wa3_dim1 = *ix3; wa3_offset = 1 + wa3_dim1; wa3 -= wa3_offset; /* Function Body */ idot = *ido / 2; if (*ido < *l1) { goto L106; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] - cc[i__ + ((k << 2) + 4) * cc_dim1]; /* L101: */ } i__2 = *ido; for (i__ = 2; i__ <= i__2; i__ += 2) { ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 4) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 2) * cc_dim1]; /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L104: */ } /* L105: */ } goto L111; L106: i__1 = *ido; for (i__ = 2; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ - 1 + ((k << 2) + 4) * cc_dim1] - cc[i__ - 1 + ((k << 2) + 2) * cc_dim1]; ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] - cc[i__ + ((k << 2) + 4) * cc_dim1]; /* L107: */ } /* L108: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] + cc[i__ + ((k << 2) + 3) * cc_dim1]; ch[i__ + (k + ch_dim2 * 3) * ch_dim1] = cc[i__ + ((k << 2) + 2) * cc_dim1] + cc[i__ + ((k << 2) + 4) * cc_dim1]; ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 2) + 1) * cc_dim1] - cc[i__ + ((k << 2) + 3) * cc_dim1]; /* L109: */ } /* L110: */ } L111: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + (ch2_dim1 << 1)] + ch2[ik + ch2_dim1 * 3]; /* L112: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + ch2_dim1 * 3] = ch2[ik + (ch2_dim1 << 1)] - ch2[ik + ch2_dim1 * 3]; /* L113: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 1)] = ch2[ik + ch2_dim1] + ch2[ik + (ch2_dim1 << 2)]; /* L114: */ } i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { ch2[ik + (ch2_dim1 << 2)] = ch2[ik + ch2_dim1] - ch2[ik + (ch2_dim1 << 2)]; /* L115: */ } for (j = 2; j <= 4; ++j) { i__1 = *l1; for (k = 1; k <= i__1; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L116: */ } /* L117: */ } if (*ido == 2) { return 0; } if (idot < *l1) { goto L120; } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L118: */ } /* L119: */ } return 0; L120: i__1 = *ido; for (i__ = 4; i__ <= i__1; i__ += 2) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[i__ + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] - wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 1)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 1)) * c1_dim1] = wa1[*l1 - 1 + (i__ - 2) * wa1_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] + wa1[*l1 + (i__ - 2) * wa1_dim1] * ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1]; c1[i__ + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1] - wa2[* ix2 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ - 1 + (k + c1_dim2 * 3) * c1_dim1] = wa2[*ix2 - 1 + (i__ - 2) * wa2_dim1] * ch[i__ - 1 + (k + ch_dim2 * 3) * ch_dim1] + wa2[*ix2 + (i__ - 2) * wa2_dim1] * ch[i__ + (k + ch_dim2 * 3) * ch_dim1]; c1[i__ + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + (i__ - 2) * wa3_dim1] * ch[i__ + (k + (ch_dim2 << 2)) * ch_dim1] - wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + ( ch_dim2 << 2)) * ch_dim1]; c1[i__ - 1 + (k + (c1_dim2 << 2)) * c1_dim1] = wa3[*ix3 - 1 + ( i__ - 2) * wa3_dim1] * ch[i__ - 1 + (k + (ch_dim2 << 2)) * ch_dim1] + wa3[*ix3 + (i__ - 2) * wa3_dim1] * ch[i__ + ( k + (ch_dim2 << 2)) * ch_dim1]; /* L121: */ } /* L122: */ } return 0; } /* passf4d_ */ /* Subroutine */ int passfd_(integer *ido, integer *ip, integer *l1, integer * idl1, doublereal *cc, doublereal *c1, doublereal *c2, doublereal *ch, doublereal *ch2, doublereal *wa) { /* System generated locals */ integer ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1, c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, i__1, i__2, i__3; /* Local variables */ static integer i__, j, k, l, jc, lc, ik, nt, l1t, idj, idl; static doublereal wai, war; static integer ipp2, idij, idlj, idot, ipph; /* Parameter adjustments */ ch_dim1 = *ido; ch_dim2 = *l1; ch_offset = 1 + ch_dim1 * (1 + ch_dim2); ch -= ch_offset; c1_dim1 = *ido; c1_dim2 = *l1; c1_offset = 1 + c1_dim1 * (1 + c1_dim2); c1 -= c1_offset; cc_dim1 = *ido; cc_dim2 = *ip; cc_offset = 1 + cc_dim1 * (1 + cc_dim2); cc -= cc_offset; ch2_dim1 = *idl1; ch2_offset = 1 + ch2_dim1; ch2 -= ch2_offset; c2_dim1 = *idl1; c2_offset = 1 + c2_dim1; c2 -= c2_offset; --wa; /* Function Body */ idot = *ido / 2; nt = *ip * *idl1; ipp2 = *ip + 2; ipph = (*ip + 1) / 2; l1t = *l1 + *l1; if (*ido < *l1) { goto L106; } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *l1; for (k = 1; k <= i__2; ++k) { i__3 = *ido; for (i__ = 1; i__ <= i__3; ++i__) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L101: */ } /* L102: */ } /* L103: */ } i__1 = *l1; for (k = 1; k <= i__1; ++k) { i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L104: */ } /* L105: */ } goto L112; L106: i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *ido; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = *l1; for (k = 1; k <= i__3; ++k) { ch[i__ + (k + j * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] + cc[i__ + (jc + k * cc_dim2) * cc_dim1]; ch[i__ + (k + jc * ch_dim2) * ch_dim1] = cc[i__ + (j + k * cc_dim2) * cc_dim1] - cc[i__ + (jc + k * cc_dim2) * cc_dim1]; /* L107: */ } /* L108: */ } /* L109: */ } i__1 = *ido; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + (k * cc_dim2 + 1) * cc_dim1]; /* L110: */ } /* L111: */ } L112: i__1 = *idl1; for (ik = 1; ik <= i__1; ++ik) { c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; /* L113: */ } idj = 0; i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; idj += *idl1; i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + j * c2_dim1] = ch2[ik + ch2_dim1] + wa[idj - 1] * ch2[ik + (ch2_dim1 << 1)]; c2[ik + jc * c2_dim1] = -wa[idj] * ch2[ik + *ip * ch2_dim1]; /* L114: */ } /* L115: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { i__2 = *idl1; for (ik = 1; ik <= i__2; ++ik) { c2[ik + c2_dim1] += ch2[ik + j * ch2_dim1]; /* L116: */ } /* L117: */ } idl = 0; i__1 = ipph; for (l = 2; l <= i__1; ++l) { lc = ipp2 - l; idl += *idl1; idlj = idl; i__2 = ipph; for (j = 3; j <= i__2; ++j) { jc = ipp2 - j; idlj += idl; if (idlj > nt) { idlj -= nt; } war = wa[idlj - 1]; wai = wa[idlj]; i__3 = *idl1; for (ik = 1; ik <= i__3; ++ik) { c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1]; c2[ik + lc * c2_dim1] -= wai * ch2[ik + jc * ch2_dim1]; /* L118: */ } /* L119: */ } /* L120: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik + jc * c2_dim1]; ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik + jc * c2_dim1]; /* L121: */ } /* L122: */ } i__1 = ipph; for (j = 2; j <= i__1; ++j) { jc = ipp2 - j; i__2 = *idl1; for (ik = 2; ik <= i__2; ik += 2) { ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc * c2_dim1]; ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc * c2_dim1]; /* L123: */ } /* L124: */ } i__1 = *ip; for (j = 2; j <= i__1; ++j) { i__2 = *l1; for (k = 1; k <= i__2; ++k) { c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * ch_dim1 + 1]; c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * ch_dim1 + 2]; /* L125: */ } /* L126: */ } if (*ido == 2) { return 0; } idj = 0; if (idot > *l1) { goto L130; } i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; idij = 0; i__2 = *ido; for (i__ = 4; i__ <= i__2; i__ += 2) { idij += idj; i__3 = *l1; for (k = 1; k <= i__3; ++k) { c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L127: */ } /* L128: */ } /* L129: */ } return 0; L130: i__1 = *ip; for (j = 2; j <= i__1; ++j) { idj += l1t; i__2 = *l1; for (k = 1; k <= i__2; ++k) { idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[ i__ - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i__ + (k + j * ch_dim2) * ch_dim1]; /* L131: */ } idij = 0; i__3 = *ido; for (i__ = 4; i__ <= i__3; i__ += 2) { idij += idj; c1[i__ + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i__ + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i__ - 1 + (k + j * ch_dim2) * ch_dim1]; /* L132: */ } /* L133: */ } /* L134: */ } return 0; } /* passfd_ */ /* Subroutine */ int cosqid_(integer *n, doublereal *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static doublereal dc, ds, dt; static integer iw1; static doublereal pih, wsk; extern /* Subroutine */ int rfftid_(integer *, doublereal *); /* Parameter adjustments */ --wsave; /* Function Body */ iw1 = *n + 1; pih = atan(1.f) * 2.f; dt = pih / (real) (*n); dc = cos(dt); ds = sin(dt); wsave[1] = dc; wsk = ds; i__1 = *n; for (k = 2; k <= i__1; ++k) { wsave[k] = dc * wsave[k - 1] - ds * wsk; wsk = ds * wsave[k - 1] + dc * wsk; /* L101: */ } rfftid_(n, &wsave[iw1]); return 0; } /* cosqid_ */ /* Subroutine */ int cosqfd_(integer *n, doublereal *x, doublereal *wsave) { /* Initialized data */ static doublereal tsq2 = 2.82842712474619; static doublereal tx; static integer iw1, iw2; static doublereal tsqx; extern /* Subroutine */ int cosqf1d_(integer *, doublereal *, doublereal * , doublereal *, doublereal *); /* Parameter adjustments */ --wsave; --x; /* Function Body */ if (*n > 2) { goto L101; } tx = x[1] + x[1]; tsqx = tsq2 * x[2]; x[1] = tx + tsqx; x[2] = tx - tsqx; return 0; L101: iw1 = *n + 1; iw2 = *n / 2 + iw1; cosqf1d_(n, &x[1], &wsave[1], &wsave[iw1], &wsave[iw2]); return 0; } /* cosqfd_ */ /* Subroutine */ int cosqf1d_(integer *n, doublereal *x, doublereal *w, doublereal *w1, doublereal *xh) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k, kc; static doublereal xn; static integer nm1, np2, ns2; extern /* Subroutine */ int rfftfd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --xh; --w1; --w; --x; /* Function Body */ ns2 = *n / 2; nm1 = *n - 1; np2 = *n + 2; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; xh[k] = x[k] + x[kc]; xh[kc] = x[k] - x[kc]; /* L101: */ } xh[ns2 + 1] = x[ns2 + 1] + x[ns2 + 1]; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; x[k] = w[k - 1] * xh[kc] + w[kc - 1] * xh[k]; x[kc] = w[k - 1] * xh[k] - w[kc - 1] * xh[kc]; /* L102: */ } x[ns2 + 1] = w[ns2] * xh[ns2 + 1]; rfftfd_(n, &x[1], &w1[1]); xn = x[2]; i__1 = nm1; for (i__ = 3; i__ <= i__1; i__ += 2) { x[i__ - 1] = x[i__] + x[i__ + 1]; x[i__] -= x[i__ + 1]; /* L103: */ } x[*n] = xn; return 0; } /* cosqf1d_ */ /* Subroutine */ int cosqbd_(integer *n, doublereal *x, doublereal *wsave) { /* Initialized data */ static doublereal tsq2 = 2.82842712474619; static doublereal x1; static integer iw1, iw2; extern /* Subroutine */ int cosqb1d_(integer *, doublereal *, doublereal * , doublereal *, doublereal *); /* Parameter adjustments */ --wsave; --x; /* Function Body */ if (*n > 2) { goto L101; } x1 = (x[1] + x[2]) * 4.f; x[2] = tsq2 * (x[1] - x[2]); x[1] = x1; return 0; L101: iw1 = *n + 1; iw2 = *n / 2 + iw1; cosqb1d_(n, &x[1], &wsave[1], &wsave[iw1], &wsave[iw2]); return 0; } /* cosqbd_ */ /* Subroutine */ int cosqb1d_(integer *n, doublereal *x, doublereal *w, doublereal *w1, doublereal *xh) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k, kc, ii; static doublereal xn; static integer nm1, np2, ns2; extern /* Subroutine */ int rfftbd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --xh; --w1; --w; --x; /* Function Body */ ns2 = *n / 2; nm1 = *n - 1; np2 = *n + 2; xn = x[*n]; i__1 = nm1; for (ii = 3; ii <= i__1; ii += 2) { i__ = np2 - ii; x[i__ + 1] = x[i__ - 1] - x[i__]; x[i__] += x[i__ - 1]; /* L101: */ } x[1] += x[1]; x[2] = xn + xn; rfftbd_(n, &x[1], &w1[1]); i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; xh[k] = w[k - 1] * x[kc] + w[kc - 1] * x[k]; xh[kc] = w[k - 1] * x[k] - w[kc - 1] * x[kc]; /* L102: */ } x[ns2 + 1] = w[ns2] * (x[ns2 + 1] + x[ns2 + 1]); i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np2 - k; x[k] = xh[k] + xh[kc]; x[kc] = xh[k] - xh[kc]; /* L103: */ } x[1] += x[1]; return 0; } /* cosqb1d_ */ /* Subroutine */ int costid_(integer *n, doublereal *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static doublereal dt, pi; static integer nm1, iw1, ns2; static doublereal dcs, wck, dss; static integer ns2m; extern /* Subroutine */ int rfftid_(integer *, doublereal *); /* Parameter adjustments */ --wsave; /* Function Body */ if (*n <= 3) { return 0; } nm1 = *n - 1; ns2 = nm1 / 2; ns2m = ns2 - 1; iw1 = ns2 + 1; pi = atan(1.f) * 4.f; dt = pi / (real) nm1; dcs = cos(dt); dss = sin(dt); wsave[1] = dss + dss; wck = dcs + dcs; if (ns2m < 2) { goto L102; } i__1 = ns2m; for (k = 2; k <= i__1; ++k) { wsave[k] = dss * wck + dcs * wsave[k - 1]; wck = dcs * wck - dss * wsave[k - 1]; /* L101: */ } L102: rfftid_(&nm1, &wsave[iw1]); return 0; } /* costid_ */ /* Subroutine */ int costd_(integer *n, doublereal *x, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k; static doublereal c1, t1, t2; static integer kc; static doublereal cn; static integer ks, nm1, np1; static doublereal fx2; static integer iw1, ns2; static doublereal tx13, x1px3; extern /* Subroutine */ int rfftfd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --wsave; --x; /* Function Body */ nm1 = *n - 1; np1 = *n + 1; ns2 = nm1 / 2; iw1 = ns2 + 1; if (*n > 3) { goto L101; } x1px3 = x[1] + x[3]; tx13 = x1px3 + x1px3; fx2 = x[2] * 4.f; x[2] = (x[1] - x[3]) * 2.f; x[1] = tx13 + fx2; x[3] = tx13 - fx2; return 0; L101: c1 = x[1] - x[*n]; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = nm1 - k; ks = ns2 - k; c1 += wsave[ks + 1] * (x[k] - x[kc + 2]); /* L102: */ } c1 += c1; x[1] += x[*n]; i__1 = ns2; for (k = 2; k <= i__1; ++k) { kc = np1 - k; t1 = x[k] + x[kc]; t2 = wsave[k - 1] * (x[k] - x[kc]); x[k] = t1 - t2; x[kc] = t1 + t2; /* L103: */ } x[ns2 + 1] += x[ns2 + 1]; rfftfd_(&nm1, &x[1], &wsave[iw1]); cn = x[2]; x[2] = c1; i__1 = nm1; for (i__ = 4; i__ <= i__1; i__ += 2) { x[i__] += x[i__ - 2]; /* L104: */ } x[*n] = cn; return 0; } /* costd_ */ /* Subroutine */ int rfftid_(integer *n, doublereal *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static doublereal dc; static integer kc; static doublereal ds, dt; static integer iw1, ns2, nqm; static doublereal tpi; extern /* Subroutine */ int cfftid_(integer *, doublereal *); /* Parameter adjustments */ --wsave; /* Function Body */ ns2 = *n / 2; nqm = (ns2 - 1) / 2; tpi = atan(1.f) * 8.f; dt = tpi / (real) (*n); dc = cos(dt); ds = sin(dt); wsave[1] = dc; wsave[ns2 - 1] = ds; if (nqm < 2) { goto L102; } i__1 = nqm; for (k = 2; k <= i__1; ++k) { kc = ns2 - k; wsave[k] = dc * wsave[k - 1] - ds * wsave[kc + 1]; wsave[kc] = ds * wsave[k - 1] + dc * wsave[kc + 1]; /* L101: */ } L102: iw1 = ns2 + 1; cfftid_(&ns2, &wsave[iw1]); return 0; } /* rfftid_ */ /* Subroutine */ int rfftfd_(integer *n, doublereal *r__, doublereal *wsave) { static doublereal r1; static integer iw1; extern /* Subroutine */ int rfftf1d_(integer *, doublereal *, doublereal * , doublereal *); /* Parameter adjustments */ --wsave; r__ -= 3; /* Function Body */ if (*n > 2) { goto L101; } r1 = (r__[3] + r__[4]) * 2.f; r__[4] = (r__[3] - r__[4]) * 2.f; r__[3] = r1; return 0; L101: iw1 = *n / 2 + 1; rfftf1d_(n, &r__[3], &wsave[iw1], &wsave[1]); return 0; } /* rfftfd_ */ /* Subroutine */ int rfftf1d_(integer *n, doublereal *x, doublereal *xh, doublereal *w) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, nq, ns2, nqm, nqp, ipar, ns2p2; static doublereal xhold1, xhold2; extern /* Subroutine */ int cfftfd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --w; xh -= 3; x -= 3; /* Function Body */ ns2 = *n / 2; ns2p2 = ns2 + 2; nq = ns2 / 2; ipar = ns2 - nq - nq; nqm = nq; if (ipar == 0) { --nqm; } nqp = nqm + 1; cfftfd_(&ns2, &x[3], &xh[3]); if (nqp < 2) { goto L107; } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = x[(k << 1) + 2] + x[(kc << 1) + 2]; xh[(kc << 1) + 2] = x[(kc << 1) + 1] - x[(k << 1) + 1]; /* L101: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(k << 1) + 1] = x[(k << 1) + 1] + x[(kc << 1) + 1]; xh[(k << 1) + 2] = x[(k << 1) + 2] - x[(kc << 1) + 2]; /* L102: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = w[k - 1] * xh[(kc << 1) + 1] + w[kc - 1] * xh[(kc << 1) + 2]; x[(kc << 1) + 2] = w[k - 1] * xh[(kc << 1) + 2] - w[kc - 1] * xh[(kc << 1) + 1]; /* L103: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = x[(kc << 1) + 1]; xh[(kc << 1) + 2] = x[(kc << 1) + 2]; /* L104: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = xh[(k << 1) + 1] - xh[(kc << 1) + 1]; x[(kc << 1) + 2] = xh[(k << 1) + 2] - xh[(kc << 1) + 2]; /* L105: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(k << 1) + 1] = xh[(k << 1) + 1] + xh[(kc << 1) + 1]; x[(k << 1) + 2] = -xh[(k << 1) + 2] - xh[(kc << 1) + 2]; /* L106: */ } if (ipar != 0) { goto L108; } L107: x[(nqp + 1 << 1) + 1] += x[(nqp + 1 << 1) + 1]; x[(nqp + 1 << 1) + 2] += x[(nqp + 1 << 1) + 2]; L108: xhold1 = x[4] + x[3]; xhold2 = x[3] - x[4]; x[4] = xhold2 + xhold2; x[3] = xhold1 + xhold1; return 0; } /* rfftf1d_ */ /* Subroutine */ int rfftbd_(integer *n, doublereal *r__, doublereal *wsave) { static doublereal r1; static integer iw1; extern /* Subroutine */ int rfftb1d_(integer *, doublereal *, doublereal * , doublereal *); /* Parameter adjustments */ --wsave; r__ -= 3; /* Function Body */ if (*n > 2) { goto L101; } r1 = r__[3] + r__[4]; r__[4] = r__[3] - r__[4]; r__[3] = r1; return 0; L101: iw1 = *n / 2 + 1; rfftb1d_(n, &r__[3], &wsave[iw1], &wsave[1]); return 0; } /* rfftbd_ */ /* Subroutine */ int rfftb1d_(integer *n, doublereal *x, doublereal *xh, doublereal *w) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, nq, ns2, nqm, nqp, ipar, ns2p2; static doublereal xhold1; extern /* Subroutine */ int cfftbd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --w; xh -= 3; x -= 3; /* Function Body */ ns2 = *n / 2; ns2p2 = ns2 + 2; nq = ns2 / 2; ipar = ns2 - nq - nq; nqm = nq; if (ipar == 0) { --nqm; } nqp = nqm + 1; xhold1 = x[3] - x[4]; x[3] = x[4] + x[3]; x[4] = xhold1; if (ipar != 0) { goto L101; } x[(nqp + 1 << 1) + 1] += x[(nqp + 1 << 1) + 1]; x[(nqp + 1 << 1) + 2] += x[(nqp + 1 << 1) + 2]; L101: if (nqp < 2) { goto L108; } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(k << 1) + 1] = x[(k << 1) + 1] + x[(kc << 1) + 1]; xh[(k << 1) + 2] = x[(kc << 1) + 2] - x[(k << 1) + 2]; /* L102: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = x[(k << 1) + 1] - x[(kc << 1) + 1]; xh[(kc << 1) + 2] = -x[(k << 1) + 2] - x[(kc << 1) + 2]; /* L103: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = xh[(kc << 1) + 1]; x[(kc << 1) + 2] = xh[(kc << 1) + 2]; /* L104: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; xh[(kc << 1) + 1] = w[k - 1] * x[(kc << 1) + 1] - w[kc - 1] * x[(kc << 1) + 2]; xh[(kc << 1) + 2] = w[k - 1] * x[(kc << 1) + 2] + w[kc - 1] * x[(kc << 1) + 1]; /* L105: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(k << 1) + 1] = xh[(k << 1) + 1] - xh[(kc << 1) + 2]; x[(k << 1) + 2] = xh[(k << 1) + 2] + xh[(kc << 1) + 1]; /* L106: */ } i__1 = nqp; for (k = 2; k <= i__1; ++k) { kc = ns2p2 - k; x[(kc << 1) + 1] = xh[(k << 1) + 1] + xh[(kc << 1) + 2]; x[(kc << 1) + 2] = xh[(kc << 1) + 1] - xh[(k << 1) + 2]; /* L107: */ } L108: cfftbd_(&ns2, &x[3], &xh[3]); return 0; } /* rfftb1d_ */ /* Subroutine */ int sinqid_(integer *n, doublereal *wsave) { extern /* Subroutine */ int cosqid_(integer *, doublereal *); /* Parameter adjustments */ --wsave; /* Function Body */ cosqid_(n, &wsave[1]); return 0; } /* sinqid_ */ /* Subroutine */ int sinqfd_(integer *n, doublereal *x, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, ns2; static doublereal xhold; extern /* Subroutine */ int cosqfd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --wsave; --x; /* Function Body */ ns2 = *n / 2; i__1 = ns2; for (k = 1; k <= i__1; ++k) { kc = *n - k; xhold = x[k]; x[k] = x[kc + 1]; x[kc + 1] = xhold; /* L101: */ } cosqfd_(n, &x[1], &wsave[1]); i__1 = *n; for (k = 2; k <= i__1; k += 2) { x[k] = -x[k]; /* L102: */ } return 0; } /* sinqfd_ */ /* Subroutine */ int sinqbd_(integer *n, doublereal *x, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer k, kc, ns2; static doublereal xhold; extern /* Subroutine */ int cosqbd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --wsave; --x; /* Function Body */ ns2 = *n / 2; i__1 = *n; for (k = 2; k <= i__1; k += 2) { x[k] = -x[k]; /* L101: */ } cosqbd_(n, &x[1], &wsave[1]); i__1 = ns2; for (k = 1; k <= i__1; ++k) { kc = *n - k; xhold = x[k]; x[k] = x[kc + 1]; x[kc + 1] = xhold; /* L102: */ } return 0; } /* sinqbd_ */ /* Subroutine */ int sintid_(integer *n, doublereal *wsave) { /* System generated locals */ integer i__1; /* Builtin functions */ double atan(doublereal), cos(doublereal), sin(doublereal); /* Local variables */ static integer k; static doublereal dt, pi; static integer np1, iw1, ns2; static doublereal dcs, wck, dss; static integer ns2m; extern /* Subroutine */ int rfftid_(integer *, doublereal *); /* Parameter adjustments */ --wsave; /* Function Body */ if (*n <= 1) { return 0; } np1 = *n + 1; ns2 = np1 / 2; ns2m = ns2 - 1; iw1 = ns2 + 1; pi = atan(1.f) * 4.f; dt = pi / (real) np1; dcs = cos(dt); dss = sin(dt); wsave[1] = dss + dss; wck = dcs + dcs; if (ns2m < 2) { goto L102; } i__1 = ns2m; for (k = 2; k <= i__1; ++k) { wsave[k] = dss * wck + dcs * wsave[k - 1]; wck = dcs * wck - dss * wsave[k - 1]; /* L101: */ } L102: rfftid_(&np1, &wsave[iw1]); return 0; } /* sintid_ */ /* Subroutine */ int sintd_(integer *n, doublereal *x, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k; static doublereal t1, t2, x1; static integer kc, np1, iw1, ns2, ns2m; extern /* Subroutine */ int rfftfd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --wsave; --x; /* Function Body */ if (*n > 1) { goto L101; } x[1] *= 4.f; return 0; L101: np1 = *n + 1; ns2 = np1 / 2; ns2m = ns2 - 1; iw1 = ns2 + 1; x1 = x[1]; x[1] = 0.f; i__1 = ns2m; for (k = 1; k <= i__1; ++k) { kc = np1 - k; t1 = x1 - x[kc]; t2 = wsave[k] * (x1 + x[kc]); x1 = x[k + 1]; x[k + 1] = t1 + t2; x[kc + 1] = t2 - t1; /* L102: */ } x[ns2 + 1] = x1 * 4.f; rfftfd_(&np1, &x[1], &wsave[iw1]); x[1] *= .5f; i__1 = *n; for (i__ = 3; i__ <= i__1; i__ += 2) { x[i__ - 1] = x[i__ + 1]; x[i__] += x[i__ - 2]; /* L103: */ } return 0; } /* sintd_ */ /* Subroutine */ int ezpowerd_(integer *n, doublereal *r__, doublereal *p, doublereal *wsave) { /* System generated locals */ integer i__1; doublereal d__1, d__2; /* Local variables */ static integer i__; static doublereal cf, cf2; static integer iw1, ns2, ns2m; extern /* Subroutine */ int rfftfd_(integer *, doublereal *, doublereal *) ; /* ****************************************************************** */ /* modified from ezfftf to just return power */ /* Parameter adjustments */ --wsave; --p; --r__; /* Function Body */ if (*n <= 2) { /* type 2,n */ /* L2: */ return 0; } /* L102: */ ns2 = *n / 2; /* modn = n-ns2-ns2 */ /* if (modn .ne. 0) go to 105 */ iw1 = *n + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { wsave[i__] = r__[i__]; /* L103: */ } rfftfd_(n, &wsave[1], &wsave[iw1]); cf = 1.f / (real) (*n); cf2 = cf * cf; /* azero = .5*cf*wsave(1) */ p[1] = cf2 * .25f * wsave[1] * wsave[1]; /* a(ns2) = .5*cf*wsave(2) */ /* b(ns2) = 0. */ p[ns2 + 1] = cf2 * .25f * wsave[2] * wsave[2]; ns2m = ns2 - 1; i__1 = ns2m; for (i__ = 1; i__ <= i__1; ++i__) { /* a(i) = cf*wsave(2*i+1) */ /* b(i) = cf*wsave(2*i+2) */ /* Computing 2nd power */ d__1 = wsave[(i__ << 1) + 2]; /* Computing 2nd power */ d__2 = wsave[(i__ << 1) + 1]; p[i__ + 1] = cf2 * (d__1 * d__1 + d__2 * d__2); /* L104: */ } return 0; /* 105 iw1 = n+n+1 */ /* do 106 i=1,n */ /* wsave(2*i-1) = r(i) */ /* wsave(2*i) = 0. */ /* 106 continue */ /* call cfftfd (n,wsave,wsave(iw1)) */ /* cf = 2./float(n) */ /* cf2=cf*cf */ /* c azero = .5*cf*wsave(1) */ /* nm1s2 = (n-1)/2 */ /* do 107 i=1,nm1s2 */ /* c a(i) = cf*wsave(2*i+1) */ /* c b(i) = -cf*wsave(2*i+2) */ /* p(i)=cf2*(wsave(2*i+2)**2+wsave(2*i+1)**2) */ /* 107 continue */ /* return */ } /* ezpowerd_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezscd_(integer *n, doublereal *r__, doublereal *b, doublereal *a, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__; static doublereal cf, cf2; static integer iw1, ns2; extern /* Subroutine */ int rfftfd_(integer *, doublereal *, doublereal *) ; /* ****************************************************************** */ /* -- modified from ezfftf to return sine and cosine series for a real */ /* - variable in a format suitable for ana, this means that first */ /* cosine is set to mean (azero here) and first and last sine are 0 */ /* the sine corresponds to b and the cosine to a */ /* we also assume that n is always even (checked by calling routine */ /* which is sc in fun2.mar */ /* - also note that call to ezffti to set up wsave is made in fun2 */ /* Parameter adjustments */ --wsave; --a; --b; --r__; /* Function Body */ if (*n <= 2) { /* type 2,n */ /* L2: */ return 0; } /* L102: */ ns2 = *n / 2; iw1 = *n + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { wsave[i__] = r__[i__]; /* L103: */ } rfftfd_(n, &wsave[1], &wsave[iw1]); cf = 1.f / (real) (*n); cf2 = cf * cf; a[1] = cf * .5f * wsave[1]; /* this was azero */ b[1] = 0.f; /* first and last sine terms are 0 */ a[ns2 + 1] = cf * .5f * wsave[2]; b[ns2 + 1] = 0.f; i__1 = ns2; for (i__ = 2; i__ <= i__1; ++i__) { a[i__] = cf * wsave[(i__ << 1) - 1]; b[i__] = cf * wsave[i__ * 2]; /* L104: */ } return 0; } /* ezscd_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezdscd_(integer *n, doublereal *r__, doublereal *b, doublereal *a, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, iw1, ns2; extern /* Subroutine */ int rfftfd_(integer *, doublereal *, doublereal *) ; /* Parameter adjustments */ --wsave; --a; --b; --r__; /* Function Body */ if (*n <= 2) { /* type 2,n */ /* L2: */ return 0; } /* L102: */ ns2 = *n / 2; iw1 = *n + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { wsave[i__] = r__[i__]; /* L103: */ } rfftfd_(n, &wsave[1], &wsave[iw1]); i__1 = ns2; for (i__ = 1; i__ <= i__1; ++i__) { a[i__] = wsave[(i__ << 1) - 1]; b[i__] = wsave[i__ * 2]; /* L104: */ } return 0; } /* ezdscd_ */ /* ---------------------------------------------------------------------- */ /* Subroutine */ int ezscbd_(integer *n, doublereal *r__, doublereal *a, doublereal *b, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, iw1, ns2; extern /* Subroutine */ int rfftbd_(integer *, doublereal *, doublereal *) ; /* -- this is a modified version of ezfftb for ana */ /* it is similar but the sine and cosine series are 1 element longer */ /* with c(1)=azero and s(1)=0 always assumed */ /* also we assume n is even */ /* Parameter adjustments */ --r__; --a; --b; --wsave; /* Function Body */ if (*n > 1) { goto L101; } r__[1] = a[1]; return 0; L101: if (*n > 2) { goto L102; } r__[1] = a[1] + a[2]; r__[2] = a[1] - a[2]; return 0; L102: ns2 = *n / 2; iw1 = *n + 1; i__1 = ns2; for (i__ = 2; i__ <= i__1; ++i__) { r__[(i__ << 1) - 1] = a[i__] * .5f; r__[i__ * 2] = b[i__] * .5f; /* L103: */ } r__[1] = a[1]; r__[2] = a[ns2 + 1]; rfftbd_(n, &r__[1], &wsave[iw1]); return 0; } /* ezscbd_ */ /* Subroutine */ int ezfftbd_(integer *n, doublereal *r__, doublereal *azero, doublereal *a, doublereal *b, doublereal *wsave) { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, ic, iw1, ns2, ns2m, modn, nm1s2; extern /* Subroutine */ int cfftbd_(integer *, doublereal *, doublereal *) , rfftbd_(integer *, doublereal *, doublereal *), ezfftid_( integer *, doublereal *); /* Parameter adjustments */ --wsave; --b; --a; --r__; /* Function Body */ if (*n > 1) { goto L101; } r__[1] = *azero; return 0; L101: if (*n > 2) { goto L102; } r__[1] = *azero + a[1]; r__[2] = *azero - a[1]; return 0; L102: ns2 = *n / 2; /* to supress repeated initialization remove the following statment */ /* ( call ezffti(n,wsave)) from both this program and ezfftf and */ /* insert it at the appropriate place in your program. */ ezfftid_(n, &wsave[1]); modn = *n - ns2 - ns2; if (modn != 0) { goto L104; } iw1 = *n + 1; ns2m = ns2 - 1; i__1 = ns2m; for (i__ = 1; i__ <= i__1; ++i__) { r__[(i__ << 1) + 1] = a[i__] * .5f; r__[(i__ << 1) + 2] = b[i__] * .5f; /* L103: */ } r__[1] = *azero; r__[2] = a[ns2]; rfftbd_(n, &r__[1], &wsave[iw1]); return 0; L104: iw1 = *n + *n + 1; nm1s2 = (*n - 1) / 2; i__1 = nm1s2; for (i__ = 1; i__ <= i__1; ++i__) { ic = *n - i__; wsave[(i__ << 1) + 1] = a[i__]; wsave[(i__ << 1) + 2] = -b[i__]; wsave[(ic << 1) + 1] = a[i__]; wsave[(ic << 1) + 2] = b[i__]; /* L105: */ } wsave[1] = *azero + *azero; wsave[2] = 0.f; cfftbd_(n, &wsave[1], &wsave[iw1]); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { r__[i__] = wsave[(i__ << 1) - 1] * .5f; /* L106: */ } return 0; } /* ezfftbd_ */