pub3dfft.C

Go to the documentation of this file.
00001 
00007 //
00008 // XXX static and global variables are unsafe for shared memory builds.
00009 // This FFT module cannot safely be used by multithreaded builds of NAMD.
00010 //
00011 
00012 #include <math.h>
00013 #include "pub3dfft.h"
00014 #ifndef M_PI
00015 #define M_PI 3.14159265358979323846
00016 #endif
00017 
00018 /* Subroutine */ int passb(int *nac, int *ido, int *ip, int *
00019         l1, int *idl1, double *cc, double *c1, double *c2, 
00020         double *ch, double *ch2, double *wa)
00021 {
00022     /* System generated locals */
00023     int ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
00024              c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, 
00025             i_1, i_2, i_3;
00026 
00027     /* Local variables */
00028     static int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, 
00029             idl, inc, idp;
00030     static double wai, war;
00031     static int ipp2;
00032 
00033     /* Parameter adjustments */
00034     cc_dim1 = *ido;
00035     cc_dim2 = *ip;
00036     cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
00037     cc -= cc_offset;
00038     c1_dim1 = *ido;
00039     c1_dim2 = *l1;
00040     c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
00041     c1 -= c1_offset;
00042     c2_dim1 = *idl1;
00043     c2_offset = c2_dim1 + 1;
00044     c2 -= c2_offset;
00045     ch_dim1 = *ido;
00046     ch_dim2 = *l1;
00047     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00048     ch -= ch_offset;
00049     ch2_dim1 = *idl1;
00050     ch2_offset = ch2_dim1 + 1;
00051     ch2 -= ch2_offset;
00052     --wa;
00053 
00054     /* Function Body */
00055     idot = *ido / 2;
00056     nt = *ip * *idl1;
00057     ipp2 = *ip + 2;
00058     ipph = (*ip + 1) / 2;
00059     idp = *ip * *ido;
00060 
00061     if (*ido < *l1) {
00062         goto L106;
00063     }
00064     i_1 = ipph;
00065     for (j = 2; j <= i_1; ++j) {
00066         jc = ipp2 - j;
00067         i_2 = *l1;
00068         for (k = 1; k <= i_2; ++k) {
00069             i_3 = *ido;
00070             for (i = 1; i <= i_3; ++i) {
00071                 ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
00072                          * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
00073                 ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k * 
00074                         cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) * 
00075                         cc_dim1];
00076 /* L101: */
00077             }
00078 /* L102: */
00079         }
00080 /* L103: */
00081     }
00082     i_1 = *l1;
00083     for (k = 1; k <= i_1; ++k) {
00084         i_2 = *ido;
00085         for (i = 1; i <= i_2; ++i) {
00086             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * 
00087                     cc_dim1];
00088 /* L104: */
00089         }
00090 /* L105: */
00091     }
00092     goto L112;
00093 L106:
00094     i_1 = ipph;
00095     for (j = 2; j <= i_1; ++j) {
00096         jc = ipp2 - j;
00097         i_2 = *ido;
00098         for (i = 1; i <= i_2; ++i) {
00099             i_3 = *l1;
00100             for (k = 1; k <= i_3; ++k) {
00101                 ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
00102                          * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
00103                 ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k * 
00104                         cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) * 
00105                         cc_dim1];
00106 /* L107: */
00107             }
00108 /* L108: */
00109         }
00110 /* L109: */
00111     }
00112     i_1 = *ido;
00113     for (i = 1; i <= i_1; ++i) {
00114         i_2 = *l1;
00115         for (k = 1; k <= i_2; ++k) {
00116             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * 
00117                     cc_dim1];
00118 /* L110: */
00119         }
00120 /* L111: */
00121     }
00122 L112:
00123     idl = 2 - *ido;
00124     inc = 0;
00125     i_1 = ipph;
00126     for (l = 2; l <= i_1; ++l) {
00127         lc = ipp2 - l;
00128         idl += *ido;
00129         i_2 = *idl1;
00130         for (ik = 1; ik <= i_2; ++ik) {
00131             c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] * ch2[ik 
00132                     + (ch2_dim1 << 1)];
00133             c2[ik + lc * c2_dim1] = wa[idl] * ch2[ik + *ip * ch2_dim1];
00134 /* L113: */
00135         }
00136         idlj = idl;
00137         inc += *ido;
00138         i_2 = ipph;
00139         for (j = 3; j <= i_2; ++j) {
00140             jc = ipp2 - j;
00141             idlj += inc;
00142             if (idlj > idp) {
00143                 idlj -= idp;
00144             }
00145             war = wa[idlj - 1];
00146             wai = wa[idlj];
00147             i_3 = *idl1;
00148             for (ik = 1; ik <= i_3; ++ik) {
00149                 c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
00150                 c2[ik + lc * c2_dim1] += wai * ch2[ik + jc * ch2_dim1];
00151 /* L114: */
00152             }
00153 /* L115: */
00154         }
00155 /* L116: */
00156     }
00157     i_1 = ipph;
00158     for (j = 2; j <= i_1; ++j) {
00159         i_2 = *idl1;
00160         for (ik = 1; ik <= i_2; ++ik) {
00161             ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
00162 /* L117: */
00163         }
00164 /* L118: */
00165     }
00166     i_1 = ipph;
00167     for (j = 2; j <= i_1; ++j) {
00168         jc = ipp2 - j;
00169         i_2 = *idl1;
00170         for (ik = 2; ik <= i_2; ik += 2) {
00171             ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik + 
00172                     jc * c2_dim1];
00173             ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik + 
00174                     jc * c2_dim1];
00175             ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc * 
00176                     c2_dim1];
00177             ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc * 
00178                     c2_dim1];
00179 /* L119: */
00180         }
00181 /* L120: */
00182     }
00183     *nac = 1;
00184     if (*ido == 2) {
00185         return 0;
00186     }
00187     *nac = 0;
00188     i_1 = *idl1;
00189     for (ik = 1; ik <= i_1; ++ik) {
00190         c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
00191 /* L121: */
00192     }
00193     i_1 = *ip;
00194     for (j = 2; j <= i_1; ++j) {
00195         i_2 = *l1;
00196         for (k = 1; k <= i_2; ++k) {
00197             c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * 
00198                     ch_dim1 + 1];
00199             c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * 
00200                     ch_dim1 + 2];
00201 /* L122: */
00202         }
00203 /* L123: */
00204     }
00205     if (idot > *l1) {
00206         goto L127;
00207     }
00208     idij = 0;
00209     i_1 = *ip;
00210     for (j = 2; j <= i_1; ++j) {
00211         idij += 2;
00212         i_2 = *ido;
00213         for (i = 4; i <= i_2; i += 2) {
00214             idij += 2;
00215             i_3 = *l1;
00216             for (k = 1; k <= i_3; ++k) {
00217                 c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i 
00218                         - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i 
00219                         + (k + j * ch_dim2) * ch_dim1];
00220                 c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
00221                         k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i - 1 + (
00222                         k + j * ch_dim2) * ch_dim1];
00223 /* L124: */
00224             }
00225 /* L125: */
00226         }
00227 /* L126: */
00228     }
00229     return 0;
00230 L127:
00231     idj = 2 - *ido;
00232     i_1 = *ip;
00233     for (j = 2; j <= i_1; ++j) {
00234         idj += *ido;
00235         i_2 = *l1;
00236         for (k = 1; k <= i_2; ++k) {
00237             idij = idj;
00238             i_3 = *ido;
00239             for (i = 4; i <= i_3; i += 2) {
00240                 idij += 2;
00241                 c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i 
00242                         - 1 + (k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i 
00243                         + (k + j * ch_dim2) * ch_dim1];
00244                 c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
00245                         k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i - 1 + (
00246                         k + j * ch_dim2) * ch_dim1];
00247 /* L128: */
00248             }
00249 /* L129: */
00250         }
00251 /* L130: */
00252     }
00253     return 0;
00254 } /* passb_ */
00255 
00256 
00257 /* Subroutine */ int passb2(int *ido, int *l1, double *cc, 
00258         double *ch, double *wa1)
00259 {
00260     /* System generated locals */
00261     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
00262 
00263     /* Local variables */
00264     static int i, k;
00265     static double ti2, tr2;
00266 
00267     /* Parameter adjustments */
00268     cc_dim1 = *ido;
00269     cc_offset = cc_dim1 * 3 + 1;
00270     cc -= cc_offset;
00271     ch_dim1 = *ido;
00272     ch_dim2 = *l1;
00273     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00274     ch -= ch_offset;
00275     --wa1;
00276 
00277     /* Function Body */
00278     if (*ido > 2) {
00279         goto L102;
00280     }
00281     i_1 = *l1;
00282     for (k = 1; k <= i_1; ++k) {
00283         ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] + 
00284                 cc[((k << 1) + 2) * cc_dim1 + 1];
00285         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 
00286                 + 1] - cc[((k << 1) + 2) * cc_dim1 + 1];
00287         ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] + 
00288                 cc[((k << 1) + 2) * cc_dim1 + 2];
00289         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 
00290                 + 2] - cc[((k << 1) + 2) * cc_dim1 + 2];
00291 /* L101: */
00292     }
00293     return 0;
00294 L102:
00295     i_1 = *l1;
00296     for (k = 1; k <= i_1; ++k) {
00297         i_2 = *ido;
00298         for (i = 2; i <= i_2; i += 2) {
00299             ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + ((k << 1) + 1) * 
00300                     cc_dim1] + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
00301             tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1] - cc[i - 1 + ((k << 1) 
00302                     + 2) * cc_dim1];
00303             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1]
00304                      + cc[i + ((k << 1) + 2) * cc_dim1];
00305             ti2 = cc[i + ((k << 1) + 1) * cc_dim1] - cc[i + ((k << 1) + 2) * 
00306                     cc_dim1];
00307             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2 + wa1[i]
00308                      * tr2;
00309             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2 - 
00310                     wa1[i] * ti2;
00311 /* L103: */
00312         }
00313 /* L104: */
00314     }
00315     return 0;
00316 } /* passb2_ */
00317 
00318 
00319 /* Subroutine */ int passb3(int *ido, int *l1, double *cc, 
00320         double *ch, double *wa1, double *wa2)
00321 {
00322     /* Initialized data */
00323 
00324     static double taur = -.5;
00325     static double taui = .866025403784439;
00326 
00327     /* System generated locals */
00328     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
00329 
00330     /* Local variables */
00331     static int i, k;
00332     static double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
00333 
00334     /* Parameter adjustments */
00335     cc_dim1 = *ido;
00336     cc_offset = (cc_dim1 << 2) + 1;
00337     cc -= cc_offset;
00338     ch_dim1 = *ido;
00339     ch_dim2 = *l1;
00340     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00341     ch -= ch_offset;
00342     --wa1;
00343     --wa2;
00344 
00345     /* Function Body */
00346     if (*ido != 2) {
00347         goto L102;
00348     }
00349     i_1 = *l1;
00350     for (k = 1; k <= i_1; ++k) {
00351         tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
00352         cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
00353         ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
00354 
00355         ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
00356         ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
00357         ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
00358 
00359         cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] - cc[(k * 3 + 3) * 
00360                 cc_dim1 + 1]);
00361         ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] - cc[(k * 3 + 3) * 
00362                 cc_dim1 + 2]);
00363         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
00364         ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
00365         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
00366         ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
00367 /* L101: */
00368     }
00369     return 0;
00370 L102:
00371     i_1 = *l1;
00372     for (k = 1; k <= i_1; ++k) {
00373         i_2 = *ido;
00374         for (i = 2; i <= i_2; i += 2) {
00375             tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1] + cc[i - 1 + (k * 3 + 3) *
00376                      cc_dim1];
00377             cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
00378             ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1) * 
00379                     cc_dim1] + tr2;
00380             ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * 
00381                     cc_dim1];
00382             ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
00383             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + 
00384                     ti2;
00385             cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1] - cc[i - 1 + (k * 
00386                     3 + 3) * cc_dim1]);
00387             ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1] - cc[i + (k * 3 + 3) *
00388                      cc_dim1]);
00389             dr2 = cr2 - ci3;
00390             dr3 = cr2 + ci3;
00391             di2 = ci2 + cr3;
00392             di3 = ci2 - cr3;
00393             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 + wa1[i]
00394                      * dr2;
00395             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 - 
00396                     wa1[i] * di2;
00397             ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 + wa2[i] * 
00398                     dr3;
00399             ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 - wa2[
00400                     i] * di3;
00401 /* L103: */
00402         }
00403 /* L104: */
00404     }
00405     return 0;
00406 } /* passb3_ */
00407 
00408 /* Subroutine */ int passb4(int *ido, int *l1, double *cc, 
00409         double *ch, double *wa1, double *wa2, double *wa3)
00410 {
00411     /* System generated locals */
00412     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
00413 
00414     /* Local variables */
00415     static int i, k;
00416     static double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, 
00417             tr2, tr3, tr4;
00418 
00419     /* Parameter adjustments */
00420     cc_dim1 = *ido;
00421     cc_offset = cc_dim1 * 5 + 1;
00422     cc -= cc_offset;
00423     ch_dim1 = *ido;
00424     ch_dim2 = *l1;
00425     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00426     ch -= ch_offset;
00427     --wa1;
00428     --wa2;
00429     --wa3;
00430 
00431     /* Function Body */
00432     if (*ido != 2) {
00433         goto L102;
00434     }
00435     i_1 = *l1;
00436     for (k = 1; k <= i_1; ++k) {
00437         ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] - cc[((k << 2) + 3) * cc_dim1 
00438                 + 2];
00439         ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] + cc[((k << 2) + 3) * cc_dim1 
00440                 + 2];
00441         tr4 = cc[((k << 2) + 4) * cc_dim1 + 2] - cc[((k << 2) + 2) * cc_dim1 
00442                 + 2];
00443         ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] + cc[((k << 2) + 4) * cc_dim1 
00444                 + 2];
00445         tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - cc[((k << 2) + 3) * cc_dim1 
00446                 + 1];
00447         tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + cc[((k << 2) + 3) * cc_dim1 
00448                 + 1];
00449         ti4 = cc[((k << 2) + 2) * cc_dim1 + 1] - cc[((k << 2) + 4) * cc_dim1 
00450                 + 1];
00451         tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] + cc[((k << 2) + 4) * cc_dim1 
00452                 + 1];
00453         ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
00454         ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
00455         ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
00456         ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
00457         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
00458         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
00459         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
00460         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
00461 /* L101: */
00462     }
00463     return 0;
00464 L102:
00465     i_1 = *l1;
00466     for (k = 1; k <= i_1; ++k) {
00467         i_2 = *ido;
00468         for (i = 2; i <= i_2; i += 2) {
00469             ti1 = cc[i + ((k << 2) + 1) * cc_dim1] - cc[i + ((k << 2) + 3) * 
00470                     cc_dim1];
00471             ti2 = cc[i + ((k << 2) + 1) * cc_dim1] + cc[i + ((k << 2) + 3) * 
00472                     cc_dim1];
00473             ti3 = cc[i + ((k << 2) + 2) * cc_dim1] + cc[i + ((k << 2) + 4) * 
00474                     cc_dim1];
00475             tr4 = cc[i + ((k << 2) + 4) * cc_dim1] - cc[i + ((k << 2) + 2) * 
00476                     cc_dim1];
00477             tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] - cc[i - 1 + ((k << 2) 
00478                     + 3) * cc_dim1];
00479             tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] + cc[i - 1 + ((k << 2) 
00480                     + 3) * cc_dim1];
00481             ti4 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] - cc[i - 1 + ((k << 2) 
00482                     + 4) * cc_dim1];
00483             tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] + cc[i - 1 + ((k << 2) 
00484                     + 4) * cc_dim1];
00485             ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
00486             cr3 = tr2 - tr3;
00487             ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
00488             ci3 = ti2 - ti3;
00489             cr2 = tr1 + tr4;
00490             cr4 = tr1 - tr4;
00491             ci2 = ti1 + ti4;
00492             ci4 = ti1 - ti4;
00493             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2 - 
00494                     wa1[i] * ci2;
00495             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2 + wa1[i]
00496                      * cr2;
00497             ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3 - wa2[
00498                     i] * ci3;
00499             ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3 + wa2[i] * 
00500                     cr3;
00501             ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4 - 
00502                     wa3[i] * ci4;
00503             ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4 + wa3[i]
00504                      * cr4;
00505 /* L103: */
00506         }
00507 /* L104: */
00508     }
00509     return 0;
00510 } /* passb4_ */
00511 
00512 /* Subroutine */ int passb5(int *ido, int *l1, double *cc, 
00513         double *ch, double *wa1, double *wa2, double *wa3, 
00514         double *wa4)
00515 {
00516     /* Initialized data */
00517 
00518     static double tr11 = .309016994374947;
00519     static double ti11 = .951056516295154;
00520     static double tr12 = -.809016994374947;
00521     static double ti12 = .587785252292473;
00522 
00523     /* System generated locals */
00524     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
00525 
00526     /* Local variables */
00527     static int i, k;
00528     static double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, 
00529             cr4, ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
00530 
00531     /* Parameter adjustments */
00532     cc_dim1 = *ido;
00533     cc_offset = cc_dim1 * 6 + 1;
00534     cc -= cc_offset;
00535     ch_dim1 = *ido;
00536     ch_dim2 = *l1;
00537     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00538     ch -= ch_offset;
00539     --wa1;
00540     --wa2;
00541     --wa3;
00542     --wa4;
00543 
00544     /* Function Body */
00545     if (*ido != 2) {
00546         goto L102;
00547     }
00548     i_1 = *l1;
00549     for (k = 1; k <= i_1; ++k) {
00550         ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
00551         ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
00552         ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
00553         ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
00554         tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
00555         tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
00556         tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
00557         tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
00558         ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + tr2 
00559                 + tr3;
00560         ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] + ti2 
00561                 + ti3;
00562         cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
00563         ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
00564         cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
00565         ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
00566         cr5 = ti11 * tr5 + ti12 * tr4;
00567         ci5 = ti11 * ti5 + ti12 * ti4;
00568         cr4 = ti12 * tr5 - ti11 * tr4;
00569         ci4 = ti12 * ti5 - ti11 * ti4;
00570         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
00571         ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
00572         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
00573         ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
00574         ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
00575         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
00576         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
00577         ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
00578 /* L101: */
00579     }
00580     return 0;
00581 L102:
00582     i_1 = *l1;
00583     for (k = 1; k <= i_1; ++k) {
00584         i_2 = *ido;
00585         for (i = 2; i <= i_2; i += 2) {
00586             ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * 
00587                     cc_dim1];
00588             ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * 
00589                     cc_dim1];
00590             ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * 
00591                     cc_dim1];
00592             ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * 
00593                     cc_dim1];
00594             tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1] - cc[i - 1 + (k * 5 + 5) *
00595                      cc_dim1];
00596             tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1] + cc[i - 1 + (k * 5 + 5) *
00597                      cc_dim1];
00598             tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1] - cc[i - 1 + (k * 5 + 4) *
00599                      cc_dim1];
00600             tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1] + cc[i - 1 + (k * 5 + 4) *
00601                      cc_dim1];
00602             ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1) * 
00603                     cc_dim1] + tr2 + tr3;
00604             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1) * cc_dim1] + 
00605                     ti2 + ti3;
00606             cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
00607 
00608             ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
00609             cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
00610 
00611             ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
00612             cr5 = ti11 * tr5 + ti12 * tr4;
00613             ci5 = ti11 * ti5 + ti12 * ti4;
00614             cr4 = ti12 * tr5 - ti11 * tr4;
00615             ci4 = ti12 * ti5 - ti11 * ti4;
00616             dr3 = cr3 - ci4;
00617             dr4 = cr3 + ci4;
00618             di3 = ci3 + cr4;
00619             di4 = ci3 - cr4;
00620             dr5 = cr2 + ci5;
00621             dr2 = cr2 - ci5;
00622             di5 = ci2 - cr5;
00623             di2 = ci2 + cr5;
00624             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 - 
00625                     wa1[i] * di2;
00626             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 + wa1[i]
00627                      * dr2;
00628             ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 - wa2[
00629                     i] * di3;
00630             ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 + wa2[i] * 
00631                     dr3;
00632             ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4 - 
00633                     wa3[i] * di4;
00634             ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4 + wa3[i]
00635                      * dr4;
00636             ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5 - wa4[
00637                     i] * di5;
00638             ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5 + wa4[i] * 
00639                     dr5;
00640 /* L103: */
00641         }
00642 /* L104: */
00643     }
00644     return 0;
00645 } /* passb5_ */
00646 
00647 /* Subroutine */ int passf(int *nac, int *ido, int *ip, int *
00648         l1, int *idl1, double *cc, double *c1, double *c2, 
00649         double *ch, double *ch2, double *wa)
00650 {
00651     /* System generated locals */
00652     int ch_dim1, ch_dim2, ch_offset, cc_dim1, cc_dim2, cc_offset, c1_dim1,
00653              c1_dim2, c1_offset, c2_dim1, c2_offset, ch2_dim1, ch2_offset, 
00654             i_1, i_2, i_3;
00655 
00656     /* Local variables */
00657     static int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, 
00658             idl, inc, idp;
00659     static double wai, war;
00660     static int ipp2;
00661 
00662     /* Parameter adjustments */
00663     cc_dim1 = *ido;
00664     cc_dim2 = *ip;
00665     cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
00666     cc -= cc_offset;
00667     c1_dim1 = *ido;
00668     c1_dim2 = *l1;
00669     c1_offset = c1_dim1 * (c1_dim2 + 1) + 1;
00670     c1 -= c1_offset;
00671     c2_dim1 = *idl1;
00672     c2_offset = c2_dim1 + 1;
00673     c2 -= c2_offset;
00674     ch_dim1 = *ido;
00675     ch_dim2 = *l1;
00676     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00677     ch -= ch_offset;
00678     ch2_dim1 = *idl1;
00679     ch2_offset = ch2_dim1 + 1;
00680     ch2 -= ch2_offset;
00681     --wa;
00682 
00683     /* Function Body */
00684     idot = *ido / 2;
00685     nt = *ip * *idl1;
00686     ipp2 = *ip + 2;
00687     ipph = (*ip + 1) / 2;
00688     idp = *ip * *ido;
00689 
00690     if (*ido < *l1) {
00691         goto L106;
00692     }
00693     i_1 = ipph;
00694     for (j = 2; j <= i_1; ++j) {
00695         jc = ipp2 - j;
00696         i_2 = *l1;
00697         for (k = 1; k <= i_2; ++k) {
00698             i_3 = *ido;
00699             for (i = 1; i <= i_3; ++i) {
00700                 ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
00701                          * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
00702                 ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k * 
00703                         cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) * 
00704                         cc_dim1];
00705 /* L101: */
00706             }
00707 /* L102: */
00708         }
00709 /* L103: */
00710     }
00711     i_1 = *l1;
00712     for (k = 1; k <= i_1; ++k) {
00713         i_2 = *ido;
00714         for (i = 1; i <= i_2; ++i) {
00715             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * 
00716                     cc_dim1];
00717 /* L104: */
00718         }
00719 /* L105: */
00720     }
00721     goto L112;
00722 L106:
00723     i_1 = ipph;
00724     for (j = 2; j <= i_1; ++j) {
00725         jc = ipp2 - j;
00726         i_2 = *ido;
00727         for (i = 1; i <= i_2; ++i) {
00728             i_3 = *l1;
00729             for (k = 1; k <= i_3; ++k) {
00730                 ch[i + (k + j * ch_dim2) * ch_dim1] = cc[i + (j + k * cc_dim2)
00731                          * cc_dim1] + cc[i + (jc + k * cc_dim2) * cc_dim1];
00732                 ch[i + (k + jc * ch_dim2) * ch_dim1] = cc[i + (j + k * 
00733                         cc_dim2) * cc_dim1] - cc[i + (jc + k * cc_dim2) * 
00734                         cc_dim1];
00735 /* L107: */
00736             }
00737 /* L108: */
00738         }
00739 /* L109: */
00740     }
00741     i_1 = *ido;
00742     for (i = 1; i <= i_1; ++i) {
00743         i_2 = *l1;
00744         for (k = 1; k <= i_2; ++k) {
00745             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * 
00746                     cc_dim1];
00747 /* L110: */
00748         }
00749 /* L111: */
00750     }
00751 L112:
00752     idl = 2 - *ido;
00753     inc = 0;
00754     i_1 = ipph;
00755     for (l = 2; l <= i_1; ++l) {
00756         lc = ipp2 - l;
00757         idl += *ido;
00758         i_2 = *idl1;
00759         for (ik = 1; ik <= i_2; ++ik) {
00760             c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] * ch2[ik 
00761                     + (ch2_dim1 << 1)];
00762             c2[ik + lc * c2_dim1] = -wa[idl] * ch2[ik + *ip * ch2_dim1];
00763 /* L113: */
00764         }
00765         idlj = idl;
00766         inc += *ido;
00767         i_2 = ipph;
00768         for (j = 3; j <= i_2; ++j) {
00769             jc = ipp2 - j;
00770             idlj += inc;
00771             if (idlj > idp) {
00772                 idlj -= idp;
00773             }
00774             war = wa[idlj - 1];
00775             wai = wa[idlj];
00776             i_3 = *idl1;
00777             for (ik = 1; ik <= i_3; ++ik) {
00778                 c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1];
00779                 c2[ik + lc * c2_dim1] -= wai * ch2[ik + jc * ch2_dim1];
00780 /* L114: */
00781             }
00782 /* L115: */
00783         }
00784 /* L116: */
00785     }
00786     i_1 = ipph;
00787     for (j = 2; j <= i_1; ++j) {
00788         i_2 = *idl1;
00789         for (ik = 1; ik <= i_2; ++ik) {
00790             ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1];
00791 /* L117: */
00792         }
00793 /* L118: */
00794     }
00795     i_1 = ipph;
00796     for (j = 2; j <= i_1; ++j) {
00797         jc = ipp2 - j;
00798         i_2 = *idl1;
00799         for (ik = 2; ik <= i_2; ik += 2) {
00800             ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - c2[ik + 
00801                     jc * c2_dim1];
00802             ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] + c2[ik + 
00803                     jc * c2_dim1];
00804             ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] + c2[ik - 1 + jc * 
00805                     c2_dim1];
00806             ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - c2[ik - 1 + jc * 
00807                     c2_dim1];
00808 /* L119: */
00809         }
00810 /* L120: */
00811     }
00812     *nac = 1;
00813     if (*ido == 2) {
00814         return 0;
00815     }
00816     *nac = 0;
00817     i_1 = *idl1;
00818     for (ik = 1; ik <= i_1; ++ik) {
00819         c2[ik + c2_dim1] = ch2[ik + ch2_dim1];
00820 /* L121: */
00821     }
00822     i_1 = *ip;
00823     for (j = 2; j <= i_1; ++j) {
00824         i_2 = *l1;
00825         for (k = 1; k <= i_2; ++k) {
00826             c1[(k + j * c1_dim2) * c1_dim1 + 1] = ch[(k + j * ch_dim2) * 
00827                     ch_dim1 + 1];
00828             c1[(k + j * c1_dim2) * c1_dim1 + 2] = ch[(k + j * ch_dim2) * 
00829                     ch_dim1 + 2];
00830 /* L122: */
00831         }
00832 /* L123: */
00833     }
00834     if (idot > *l1) {
00835         goto L127;
00836     }
00837     idij = 0;
00838     i_1 = *ip;
00839     for (j = 2; j <= i_1; ++j) {
00840         idij += 2;
00841         i_2 = *ido;
00842         for (i = 4; i <= i_2; i += 2) {
00843             idij += 2;
00844             i_3 = *l1;
00845             for (k = 1; k <= i_3; ++k) {
00846                 c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i 
00847                         - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i 
00848                         + (k + j * ch_dim2) * ch_dim1];
00849                 c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
00850                         k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i - 1 + (
00851                         k + j * ch_dim2) * ch_dim1];
00852 /* L124: */
00853             }
00854 /* L125: */
00855         }
00856 /* L126: */
00857     }
00858     return 0;
00859 L127:
00860     idj = 2 - *ido;
00861     i_1 = *ip;
00862     for (j = 2; j <= i_1; ++j) {
00863         idj += *ido;
00864         i_2 = *l1;
00865         for (k = 1; k <= i_2; ++k) {
00866             idij = idj;
00867             i_3 = *ido;
00868             for (i = 4; i <= i_3; i += 2) {
00869                 idij += 2;
00870                 c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i 
00871                         - 1 + (k + j * ch_dim2) * ch_dim1] + wa[idij] * ch[i 
00872                         + (k + j * ch_dim2) * ch_dim1];
00873                 c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] * ch[i + (
00874                         k + j * ch_dim2) * ch_dim1] - wa[idij] * ch[i - 1 + (
00875                         k + j * ch_dim2) * ch_dim1];
00876 /* L128: */
00877             }
00878 /* L129: */
00879         }
00880 /* L130: */
00881     }
00882     return 0;
00883 } /* passf_ */
00884 
00885 /* Subroutine */ int passf2(int *ido, int *l1, double *cc, 
00886         double *ch, double *wa1)
00887 {
00888     /* System generated locals */
00889     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
00890 
00891     /* Local variables */
00892     static int i, k;
00893     static double ti2, tr2;
00894 
00895     /* Parameter adjustments */
00896     cc_dim1 = *ido;
00897     cc_offset = cc_dim1 * 3 + 1;
00898     cc -= cc_offset;
00899     ch_dim1 = *ido;
00900     ch_dim2 = *l1;
00901     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00902     ch -= ch_offset;
00903     --wa1;
00904 
00905     /* Function Body */
00906     if (*ido > 2) {
00907         goto L102;
00908     }
00909     i_1 = *l1;
00910     for (k = 1; k <= i_1; ++k) {
00911         ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] + 
00912                 cc[((k << 1) + 2) * cc_dim1 + 1];
00913         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 
00914                 + 1] - cc[((k << 1) + 2) * cc_dim1 + 1];
00915         ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] + 
00916                 cc[((k << 1) + 2) * cc_dim1 + 2];
00917         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 
00918                 + 2] - cc[((k << 1) + 2) * cc_dim1 + 2];
00919 /* L101: */
00920     }
00921     return 0;
00922 L102:
00923     i_1 = *l1;
00924     for (k = 1; k <= i_1; ++k) {
00925         i_2 = *ido;
00926         for (i = 2; i <= i_2; i += 2) {
00927             ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + ((k << 1) + 1) * 
00928                     cc_dim1] + cc[i - 1 + ((k << 1) + 2) * cc_dim1];
00929             tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1] - cc[i - 1 + ((k << 1) 
00930                     + 2) * cc_dim1];
00931             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1]
00932                      + cc[i + ((k << 1) + 2) * cc_dim1];
00933             ti2 = cc[i + ((k << 1) + 1) * cc_dim1] - cc[i + ((k << 1) + 2) * 
00934                     cc_dim1];
00935             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2 - wa1[i]
00936                      * tr2;
00937             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2 + 
00938                     wa1[i] * ti2;
00939 /* L103: */
00940         }
00941 /* L104: */
00942     }
00943     return 0;
00944 } /* passf2_ */
00945 
00946 /* Subroutine */ int passf3(int *ido, int *l1, double *cc, 
00947         double *ch, double *wa1, double *wa2)
00948 {
00949     /* Initialized data */
00950 
00951     static double taur = -.5;
00952     static double taui = -.866025403784439;
00953 
00954     /* System generated locals */
00955     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
00956 
00957     /* Local variables */
00958     static int i, k;
00959     static double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
00960 
00961     /* Parameter adjustments */
00962     cc_dim1 = *ido;
00963     cc_offset = (cc_dim1 << 2) + 1;
00964     cc -= cc_offset;
00965     ch_dim1 = *ido;
00966     ch_dim2 = *l1;
00967     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
00968     ch -= ch_offset;
00969     --wa1;
00970     --wa2;
00971 
00972     /* Function Body */
00973     if (*ido != 2) {
00974         goto L102;
00975     }
00976     i_1 = *l1;
00977     for (k = 1; k <= i_1; ++k) {
00978         tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1];
00979         cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2;
00980         ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2;
00981 
00982         ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2];
00983         ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2;
00984         ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2;
00985 
00986         cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] - cc[(k * 3 + 3) * 
00987                 cc_dim1 + 1]);
00988         ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] - cc[(k * 3 + 3) * 
00989                 cc_dim1 + 2]);
00990         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3;
00991         ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3;
00992         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3;
00993         ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3;
00994 /* L101: */
00995     }
00996     return 0;
00997 L102:
00998     i_1 = *l1;
00999     for (k = 1; k <= i_1; ++k) {
01000         i_2 = *ido;
01001         for (i = 2; i <= i_2; i += 2) {
01002             tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1] + cc[i - 1 + (k * 3 + 3) *
01003                      cc_dim1];
01004             cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2;
01005             ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1) * 
01006                     cc_dim1] + tr2;
01007             ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * 
01008                     cc_dim1];
01009             ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2;
01010             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + 
01011                     ti2;
01012             cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1] - cc[i - 1 + (k * 
01013                     3 + 3) * cc_dim1]);
01014             ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1] - cc[i + (k * 3 + 3) *
01015                      cc_dim1]);
01016             dr2 = cr2 - ci3;
01017             dr3 = cr2 + ci3;
01018             di2 = ci2 + cr3;
01019             di3 = ci2 - cr3;
01020             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 - wa1[i]
01021                      * dr2;
01022             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 + 
01023                     wa1[i] * di2;
01024             ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 - wa2[i] * 
01025                     dr3;
01026             ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 + wa2[
01027                     i] * di3;
01028 /* L103: */
01029         }
01030 /* L104: */
01031     }
01032     return 0;
01033 } /* passf3_ */
01034 
01035 /* Subroutine */ int passf4(int *ido, int *l1, double *cc, 
01036         double *ch, double *wa1, double *wa2, double *wa3)
01037 {
01038     /* System generated locals */
01039     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
01040 
01041     /* Local variables */
01042     static int i, k;
01043     static double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, 
01044             tr2, tr3, tr4;
01045 
01046     /* Parameter adjustments */
01047     cc_dim1 = *ido;
01048     cc_offset = cc_dim1 * 5 + 1;
01049     cc -= cc_offset;
01050     ch_dim1 = *ido;
01051     ch_dim2 = *l1;
01052     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
01053     ch -= ch_offset;
01054     --wa1;
01055     --wa2;
01056     --wa3;
01057 
01058     /* Function Body */
01059     if (*ido != 2) {
01060         goto L102;
01061     }
01062     i_1 = *l1;
01063     for (k = 1; k <= i_1; ++k) {
01064         ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] - cc[((k << 2) + 3) * cc_dim1 
01065                 + 2];
01066         ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] + cc[((k << 2) + 3) * cc_dim1 
01067                 + 2];
01068         tr4 = cc[((k << 2) + 2) * cc_dim1 + 2] - cc[((k << 2) + 4) * cc_dim1 
01069                 + 2];
01070         ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] + cc[((k << 2) + 4) * cc_dim1 
01071                 + 2];
01072         tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - cc[((k << 2) + 3) * cc_dim1 
01073                 + 1];
01074         tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] + cc[((k << 2) + 3) * cc_dim1 
01075                 + 1];
01076         ti4 = cc[((k << 2) + 4) * cc_dim1 + 1] - cc[((k << 2) + 2) * cc_dim1 
01077                 + 1];
01078         tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] + cc[((k << 2) + 4) * cc_dim1 
01079                 + 1];
01080         ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3;
01081         ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3;
01082         ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3;
01083         ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3;
01084         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4;
01085         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4;
01086         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4;
01087         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4;
01088 /* L101: */
01089     }
01090     return 0;
01091 L102:
01092     i_1 = *l1;
01093     for (k = 1; k <= i_1; ++k) {
01094         i_2 = *ido;
01095         for (i = 2; i <= i_2; i += 2) {
01096             ti1 = cc[i + ((k << 2) + 1) * cc_dim1] - cc[i + ((k << 2) + 3) * 
01097                     cc_dim1];
01098             ti2 = cc[i + ((k << 2) + 1) * cc_dim1] + cc[i + ((k << 2) + 3) * 
01099                     cc_dim1];
01100             ti3 = cc[i + ((k << 2) + 2) * cc_dim1] + cc[i + ((k << 2) + 4) * 
01101                     cc_dim1];
01102             tr4 = cc[i + ((k << 2) + 2) * cc_dim1] - cc[i + ((k << 2) + 4) * 
01103                     cc_dim1];
01104             tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] - cc[i - 1 + ((k << 2) 
01105                     + 3) * cc_dim1];
01106             tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] + cc[i - 1 + ((k << 2) 
01107                     + 3) * cc_dim1];
01108             ti4 = cc[i - 1 + ((k << 2) + 4) * cc_dim1] - cc[i - 1 + ((k << 2) 
01109                     + 2) * cc_dim1];
01110             tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] + cc[i - 1 + ((k << 2) 
01111                     + 4) * cc_dim1];
01112             ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3;
01113             cr3 = tr2 - tr3;
01114             ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3;
01115             ci3 = ti2 - ti3;
01116             cr2 = tr1 + tr4;
01117             cr4 = tr1 - tr4;
01118             ci2 = ti1 + ti4;
01119             ci4 = ti1 - ti4;
01120             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2 + 
01121                     wa1[i] * ci2;
01122             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2 - wa1[i]
01123                      * cr2;
01124             ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3 + wa2[
01125                     i] * ci3;
01126             ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3 - wa2[i] * 
01127                     cr3;
01128             ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4 + 
01129                     wa3[i] * ci4;
01130             ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4 - wa3[i]
01131                      * cr4;
01132 /* L103: */
01133         }
01134 /* L104: */
01135     }
01136     return 0;
01137 } /* passf4_ */
01138 
01139 /* Subroutine */ int passf5(int *ido, int *l1, double *cc, 
01140         double *ch, double *wa1, double *wa2, double *wa3, 
01141         double *wa4)
01142 {
01143     /* Initialized data */
01144 
01145     static double tr11 = .309016994374947;
01146     static double ti11 = -.951056516295154;
01147     static double tr12 = -.809016994374947;
01148     static double ti12 = -.587785252292473;
01149 
01150     /* System generated locals */
01151     int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2;
01152 
01153     /* Local variables */
01154     static int i, k;
01155     static double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, 
01156             cr4, ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
01157 
01158     /* Parameter adjustments */
01159     cc_dim1 = *ido;
01160     cc_offset = cc_dim1 * 6 + 1;
01161     cc -= cc_offset;
01162     ch_dim1 = *ido;
01163     ch_dim2 = *l1;
01164     ch_offset = ch_dim1 * (ch_dim2 + 1) + 1;
01165     ch -= ch_offset;
01166     --wa1;
01167     --wa2;
01168     --wa3;
01169     --wa4;
01170 
01171     /* Function Body */
01172     if (*ido != 2) {
01173         goto L102;
01174     }
01175     i_1 = *l1;
01176     for (k = 1; k <= i_1; ++k) {
01177         ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2];
01178         ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2];
01179         ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2];
01180         ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2];
01181         tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1];
01182         tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1];
01183         tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1];
01184         tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1];
01185         ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] + tr2 
01186                 + tr3;
01187         ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] + ti2 
01188                 + ti3;
01189         cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3;
01190         ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3;
01191         cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3;
01192         ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3;
01193         cr5 = ti11 * tr5 + ti12 * tr4;
01194         ci5 = ti11 * ti5 + ti12 * ti4;
01195         cr4 = ti12 * tr5 - ti11 * tr4;
01196         ci4 = ti12 * ti5 - ti11 * ti4;
01197         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5;
01198         ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5;
01199         ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5;
01200         ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4;
01201         ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4;
01202         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4;
01203         ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4;
01204         ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5;
01205 /* L101: */
01206     }
01207     return 0;
01208 L102:
01209     i_1 = *l1;
01210     for (k = 1; k <= i_1; ++k) {
01211         i_2 = *ido;
01212         for (i = 2; i <= i_2; i += 2) {
01213             ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * 
01214                     cc_dim1];
01215             ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * 
01216                     cc_dim1];
01217             ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * 
01218                     cc_dim1];
01219             ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * 
01220                     cc_dim1];
01221             tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1] - cc[i - 1 + (k * 5 + 5) *
01222                      cc_dim1];
01223             tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1] + cc[i - 1 + (k * 5 + 5) *
01224                      cc_dim1];
01225             tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1] - cc[i - 1 + (k * 5 + 4) *
01226                      cc_dim1];
01227             tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1] + cc[i - 1 + (k * 5 + 4) *
01228                      cc_dim1];
01229             ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1) * 
01230                     cc_dim1] + tr2 + tr3;
01231             ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1) * cc_dim1] + 
01232                     ti2 + ti3;
01233             cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3;
01234 
01235             ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3;
01236             cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3;
01237 
01238             ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3;
01239             cr5 = ti11 * tr5 + ti12 * tr4;
01240             ci5 = ti11 * ti5 + ti12 * ti4;
01241             cr4 = ti12 * tr5 - ti11 * tr4;
01242             ci4 = ti12 * ti5 - ti11 * ti4;
01243             dr3 = cr3 - ci4;
01244             dr4 = cr3 + ci4;
01245             di3 = ci3 + cr4;
01246             di4 = ci3 - cr4;
01247             dr5 = cr2 + ci5;
01248             dr2 = cr2 - ci5;
01249             di5 = ci2 - cr5;
01250             di2 = ci2 + cr5;
01251             ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 + 
01252                     wa1[i] * di2;
01253             ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 - wa1[i]
01254                      * dr2;
01255             ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 + wa2[
01256                     i] * di3;
01257             ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 - wa2[i] * 
01258                     dr3;
01259             ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4 + 
01260                     wa3[i] * di4;
01261             ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4 - wa3[i]
01262                      * dr4;
01263             ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5 + wa4[
01264                     i] * di5;
01265             ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5 - wa4[i] * 
01266                     dr5;
01267 /* L103: */
01268         }
01269 /* L104: */
01270     }
01271     return 0;
01272 } /* passf5_ */
01273 
01274 /* Subroutine */ int passb4(int *ido, int *l1, double *cc,
01275         double *ch, double *wa1, double *wa2, double *wa3);
01276 
01277 /* Subroutine */ int cfftb1(int *n, double *c, double *ch, 
01278         double *wa, int *ifac)
01279 {
01280     /* System generated locals */
01281     int i_1;
01282 
01283     /* Local variables */
01284     static int idot, i;
01285     static int k1, l1, l2, n2;
01286     static int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
01287 
01288     /* Parameter adjustments */
01289     --c;
01290     --ch;
01291     --wa;
01292     --ifac;
01293 
01294     /* Function Body */
01295     nf = ifac[2];
01296     na = 0;
01297     l1 = 1;
01298     iw = 1;
01299     i_1 = nf;
01300     for (k1 = 1; k1 <= i_1; ++k1) {
01301         ip = ifac[k1 + 2];
01302         l2 = ip * l1;
01303         ido = *n / l2;
01304         idot = ido + ido;
01305         idl1 = idot * l1;
01306         if (ip != 4) {
01307             goto L103;
01308         }
01309         ix2 = iw + idot;
01310         ix3 = ix2 + idot;
01311         if (na != 0) {
01312             goto L101;
01313         }
01314         passb4(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
01315         goto L102;
01316 L101:
01317         passb4(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3]);
01318 L102:
01319         na = 1 - na;
01320         goto L115;
01321 L103:
01322         if (ip != 2) {
01323             goto L106;
01324         }
01325         if (na != 0) {
01326             goto L104;
01327         }
01328         passb2(&idot, &l1, &c[1], &ch[1], &wa[iw]);
01329         goto L105;
01330 L104:
01331         passb2(&idot, &l1, &ch[1], &c[1], &wa[iw]);
01332 L105:
01333         na = 1 - na;
01334         goto L115;
01335 L106:
01336         if (ip != 3) {
01337             goto L109;
01338         }
01339         ix2 = iw + idot;
01340         if (na != 0) {
01341             goto L107;
01342         }
01343         passb3(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2]);
01344         goto L108;
01345 L107:
01346         passb3(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2]);
01347 L108:
01348         na = 1 - na;
01349         goto L115;
01350 L109:
01351         if (ip != 5) {
01352             goto L112;
01353         }
01354         ix2 = iw + idot;
01355         ix3 = ix2 + idot;
01356         ix4 = ix3 + idot;
01357         if (na != 0) {
01358             goto L110;
01359         }
01360         passb5(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
01361                 ix4]);
01362         goto L111;
01363 L110:
01364         passb5(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
01365                 ix4]);
01366 L111:
01367         na = 1 - na;
01368         goto L115;
01369 L112:
01370         if (na != 0) {
01371             goto L113;
01372         }
01373         passb(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], &ch[
01374                 1], &wa[iw]);
01375         goto L114;
01376 L113:
01377         passb(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], &
01378                 c[1], &wa[iw]);
01379 L114:
01380         if (nac != 0) {
01381             na = 1 - na;
01382         }
01383 L115:
01384         l1 = l2;
01385         iw += (ip - 1) * idot;
01386 /* L116: */
01387     }
01388     if (na == 0) {
01389         return 0;
01390     }
01391     n2 = *n + *n;
01392     i_1 = n2;
01393     for (i = 1; i <= i_1; ++i) {
01394         c[i] = ch[i];
01395 /* L117: */
01396     }
01397     return 0;
01398 } /* cfftb1_ */
01399 
01400 
01401 /* Subroutine */ int cfftb(int *n, double *c, double *wsave)
01402 {
01403     static int iw1, iw2;
01404 
01405     /* Parameter adjustments */
01406     --c;
01407     --wsave;
01408 
01409     /* Function Body */
01410     if (*n == 1) {
01411         return 0;
01412     }
01413     iw1 = *n + *n + 1;
01414     iw2 = iw1 + *n + *n;
01415     cfftb1(n, &c[1], &wsave[1], &wsave[iw1], (int*) &wsave[iw2]);
01416     return 0;
01417 } /* cfftb_ */
01418 
01419 
01420 /* Subroutine */ int cfftf1(int *n, double *c, double *ch, 
01421         double *wa, int *ifac)
01422 {
01423     /* System generated locals */
01424     int i_1;
01425 
01426     /* Local variables */
01427     static int idot, i;
01428     static int k1, l1, l2, n2;
01429     static int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
01430 
01431     /* Parameter adjustments */
01432     --c;
01433     --ch;
01434     --wa;
01435     --ifac;
01436 
01437     /* Function Body */
01438     nf = ifac[2];
01439     na = 0;
01440     l1 = 1;
01441     iw = 1;
01442     i_1 = nf;
01443     for (k1 = 1; k1 <= i_1; ++k1) {
01444         ip = ifac[k1 + 2];
01445         l2 = ip * l1;
01446         ido = *n / l2;
01447         idot = ido + ido;
01448         idl1 = idot * l1;
01449         if (ip != 4) {
01450             goto L103;
01451         }
01452         ix2 = iw + idot;
01453         ix3 = ix2 + idot;
01454         if (na != 0) {
01455             goto L101;
01456         }
01457         passf4(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
01458         goto L102;
01459 L101:
01460         passf4(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3]);
01461 L102:
01462         na = 1 - na;
01463         goto L115;
01464 L103:
01465         if (ip != 2) {
01466             goto L106;
01467         }
01468         if (na != 0) {
01469             goto L104;
01470         }
01471         passf2(&idot, &l1, &c[1], &ch[1], &wa[iw]);
01472         goto L105;
01473 L104:
01474         passf2(&idot, &l1, &ch[1], &c[1], &wa[iw]);
01475 L105:
01476         na = 1 - na;
01477         goto L115;
01478 L106:
01479         if (ip != 3) {
01480             goto L109;
01481         }
01482         ix2 = iw + idot;
01483         if (na != 0) {
01484             goto L107;
01485         }
01486         passf3(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2]);
01487         goto L108;
01488 L107:
01489         passf3(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2]);
01490 L108:
01491         na = 1 - na;
01492         goto L115;
01493 L109:
01494         if (ip != 5) {
01495             goto L112;
01496         }
01497         ix2 = iw + idot;
01498         ix3 = ix2 + idot;
01499         ix4 = ix3 + idot;
01500         if (na != 0) {
01501             goto L110;
01502         }
01503         passf5(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
01504                 ix4]);
01505         goto L111;
01506 L110:
01507         passf5(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
01508                 ix4]);
01509 L111:
01510         na = 1 - na;
01511         goto L115;
01512 L112:
01513         if (na != 0) {
01514             goto L113;
01515         }
01516         passf(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], &ch[
01517                 1], &wa[iw]);
01518         goto L114;
01519 L113:
01520         passf(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], &
01521                 c[1], &wa[iw]);
01522 L114:
01523         if (nac != 0) {
01524             na = 1 - na;
01525         }
01526 L115:
01527         l1 = l2;
01528         iw += (ip - 1) * idot;
01529 /* L116: */
01530     }
01531     if (na == 0) {
01532         return 0;
01533     }
01534     n2 = *n + *n;
01535     i_1 = n2;
01536     for (i = 1; i <= i_1; ++i) {
01537         c[i] = ch[i];
01538 /* L117: */
01539     }
01540     return 0;
01541 } /* cfftf1_ */
01542 
01543 
01544 /* Subroutine */ int cfftf(int *n, double *c, double *wsave)
01545 {
01546     static int iw1, iw2;
01547 
01548     /* Parameter adjustments */
01549     --c;
01550     --wsave;
01551 
01552     /* Function Body */
01553     if (*n == 1) {
01554         return 0;
01555     }
01556     iw1 = *n + *n + 1;
01557     iw2 = iw1 + *n + *n;
01558     cfftf1(n, &c[1], &wsave[1], &wsave[iw1], (int*) &wsave[iw2]);
01559     return 0;
01560 } /* cfftf_ */
01561 
01562 
01563 /* Subroutine */ int cffti1(int *n, double *wa, int *ifac)
01564 {
01565     /* Initialized data */
01566 
01567     static int ntryh[4] = { 3,4,2,5 };
01568 
01569     /* System generated locals */
01570     int i_1, i_2, i_3;
01571 
01572     /* Local variables */
01573     static double argh;
01574     static int idot, ntry, i, j;
01575     static double argld;
01576     static int i1, k1, l1, l2, ib;
01577     static double fi;
01578     static int ld, ii, nf, ip, nl, nq, nr;
01579     static double arg;
01580     static int ido, ipm;
01581     static double tpi;
01582 
01583     /* Parameter adjustments */
01584     --wa;
01585     --ifac;
01586 
01587     /* Function Body */
01588     nl = *n;
01589     nf = 0;
01590     j = 0;
01591 L101:
01592     ++j;
01593     if (j - 4 <= 0) {
01594         goto L102;
01595     } else {
01596         goto L103;
01597     }
01598 L102:
01599     ntry = ntryh[j - 1];
01600     goto L104;
01601 L103:
01602     ntry += 2;
01603 L104:
01604     nq = nl / ntry;
01605     nr = nl - ntry * nq;
01606     if (nr != 0) {
01607         goto L101;
01608     } else {
01609         goto L105;
01610     }
01611 L105:
01612     ++nf;
01613     ifac[nf + 2] = ntry;
01614     nl = nq;
01615     if (ntry != 2) {
01616         goto L107;
01617     }
01618     if (nf == 1) {
01619         goto L107;
01620     }
01621     i_1 = nf;
01622     for (i = 2; i <= i_1; ++i) {
01623         ib = nf - i + 2;
01624         ifac[ib + 2] = ifac[ib + 1];
01625 /* L106: */
01626     }
01627     ifac[3] = 2;
01628 L107:
01629     if (nl != 1) {
01630         goto L104;
01631     }
01632     ifac[1] = *n;
01633     ifac[2] = nf;
01634     tpi = 6.28318530717959;
01635     argh = tpi / (double) (*n);
01636     i = 2;
01637     l1 = 1;
01638     i_1 = nf;
01639     for (k1 = 1; k1 <= i_1; ++k1) {
01640         ip = ifac[k1 + 2];
01641         ld = 0;
01642         l2 = l1 * ip;
01643         ido = *n / l2;
01644         idot = ido + ido + 2;
01645         ipm = ip - 1;
01646         i_2 = ipm;
01647         for (j = 1; j <= i_2; ++j) {
01648             i1 = i;
01649             wa[i - 1] = 1.;
01650             wa[i] = 0.;
01651             ld += l1;
01652             fi = 0.;
01653             argld = (double) ld * argh;
01654             i_3 = idot;
01655             for (ii = 4; ii <= i_3; ii += 2) {
01656                 i += 2;
01657                 fi += 1.;
01658                 arg = fi * argld;
01659                 wa[i - 1] = cos(arg);
01660                 wa[i] = sin(arg);
01661 /* L108: */
01662             }
01663             if (ip <= 5) {
01664                 goto L109;
01665             }
01666             wa[i1 - 1] = wa[i - 1];
01667             wa[i1] = wa[i];
01668 L109:
01669         ;}
01670         l1 = l2;
01671 /* L110: */
01672     }
01673     return 0;
01674 } /* cffti1_ */
01675 
01676 /* Subroutine */ int cffti(int *n, double *wsave)
01677 {
01678     static int iw1, iw2;
01679 
01680     /* Parameter adjustments */
01681     --wsave;
01682 
01683     /* Function Body */
01684     if (*n == 1) {
01685         return 0;
01686     }
01687     iw1 = *n + *n + 1;
01688     iw2 = iw1 + *n + *n;
01689     cffti1(n, &wsave[iw1], (int*) &wsave[iw2]);
01690     return 0;
01691 } /* cffti_ */
01692 
01693 
01694 /* typedef struct { double r, i; } doublecomplex; */
01695 
01696 /****************************************************************************
01697 **/
01698 
01699 /*      3D (slow) Fourier Transform */
01700 /*   this 1d->3d code is brute force approach */
01701 /*   the 1d code is a double precision version of fftpack from netlib */
01702 /*   due to Paul N Swartztrauber at NCAR Boulder Coloraso */
01703 
01704 /****************************************************************************
01705 **/
01706 /* Subroutine */ int pubz3di(int *n1, int *n2, int *n3, 
01707         double *table, int *ntable)
01708 {
01709     /* System generated locals */
01710     int table_dim1, table_offset;
01711 
01712     /* Local variables */
01713 
01714     /* Parameter adjustments */
01715     table_dim1 = *ntable;
01716     table_offset = table_dim1 + 1;
01717     table -= table_offset;
01718 
01719     /* Function Body */
01720 /* ntable should be 4*max(n1,n2,n3) +15 */
01721     cffti(n1, &table[table_dim1 + 1]);
01722     cffti(n2, &table[(table_dim1 << 1) + 1]);
01723     cffti(n3, &table[table_dim1 * 3 + 1]);
01724     return 0;
01725 } /* pubz3di_ */
01726 
01727 /****************************************************************************
01728 **/
01729 /* Subroutine */ int pubz3d(int *isign, int *n1, int *n2, 
01730         int *n3, doublecomplex *w, int *ld1, int *ld2, double 
01731         *table, int *ntable, doublecomplex *work /*, int * DUMMY1 nwork */)
01732 {
01733     /* System generated locals */
01734     int w_dim1, w_dim2, w_offset, table_dim1, table_offset, i_1, i_2, i_3,
01735              i_4, i_5;
01736 
01737     /* Local variables */
01738     static int i, j, k;
01739 
01740     /* Parameter adjustments */
01741     w_dim1 = *ld1;
01742     w_dim2 = *ld2;
01743     w_offset = w_dim1 * (w_dim2 + 1) + 1;
01744     w -= w_offset;
01745     table_dim1 = *ntable;
01746     table_offset = table_dim1 + 1;
01747     table -= table_offset;
01748     --work;
01749 
01750     /* Function Body */
01751 /* ntable should be 4*max(n1,n2,n3) +15 */
01752 /* nwork should be max(n1,n2,n3) */
01753 
01754 /*   transform along X  first ... */
01755 
01756     i_1 = *n3;
01757     for (k = 1; k <= i_1; ++k) {
01758         i_2 = *n2;
01759         for (j = 1; j <= i_2; ++j) {
01760             i_3 = *n1;
01761             for (i = 1; i <= i_3; ++i) {
01762                 i_4 = i;
01763                 i_5 = i + (j + k * w_dim2) * w_dim1;
01764                 work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
01765 /* L70: */
01766             }
01767             if (*isign == -1) {
01768                 cfftf(n1, (double*) &work[1], &table[table_dim1 + 1]);
01769             }
01770             if (*isign == 1) {
01771                 cfftb(n1, (double*) &work[1], &table[table_dim1 + 1]);
01772             }
01773             i_3 = *n1;
01774             for (i = 1; i <= i_3; ++i) {
01775                 i_4 = i + (j + k * w_dim2) * w_dim1;
01776                 i_5 = i;
01777                 w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
01778 /* L80: */
01779             }
01780 /* L90: */
01781         }
01782 /* L100: */
01783     }
01784 
01785 /*   transform along Y then ... */
01786 
01787     i_1 = *n3;
01788     for (k = 1; k <= i_1; ++k) {
01789         i_2 = *n1;
01790         for (i = 1; i <= i_2; ++i) {
01791             i_3 = *n2;
01792             for (j = 1; j <= i_3; ++j) {
01793                 i_4 = j;
01794                 i_5 = i + (j + k * w_dim2) * w_dim1;
01795                 work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
01796 /* L170: */
01797             }
01798             if (*isign == -1) {
01799                 cfftf(n2, (double*) &work[1], &table[(table_dim1 << 1) + 1]);
01800             }
01801             if (*isign == 1) {
01802                 cfftb(n2, (double*) &work[1], &table[(table_dim1 << 1) + 1]);
01803             }
01804             i_3 = *n2;
01805             for (j = 1; j <= i_3; ++j) {
01806                 i_4 = i + (j + k * w_dim2) * w_dim1;
01807                 i_5 = j;
01808                 w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
01809 /* L180: */
01810             }
01811 /* L190: */
01812         }
01813 /* L200: */
01814     }
01815 
01816 /*   transform along Z finally ... */
01817 
01818     i_1 = *n1;
01819     for (i = 1; i <= i_1; ++i) {
01820         i_2 = *n2;
01821         for (j = 1; j <= i_2; ++j) {
01822             i_3 = *n3;
01823             for (k = 1; k <= i_3; ++k) {
01824                 i_4 = k;
01825                 i_5 = i + (j + k * w_dim2) * w_dim1;
01826                 work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
01827 /* L270: */
01828             }
01829             if (*isign == -1) {
01830                 cfftf(n3, (double*) &work[1], &table[table_dim1 * 3 + 1]);
01831             }
01832             if (*isign == 1) {
01833                 cfftb(n3, (double*) &work[1], &table[table_dim1 * 3 + 1]);
01834             }
01835             i_3 = *n3;
01836             for (k = 1; k <= i_3; ++k) {
01837                 i_4 = i + (j + k * w_dim2) * w_dim1;
01838                 i_5 = k;
01839                 w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
01840 /* L280: */
01841             }
01842 /* L290: */
01843         }
01844 /* L300: */
01845     }
01846     return 0;
01847 } /* pubz3d_ */
01848 
01849 /****************************************************************************/
01850 
01851 int pubd3di(int n1, int n2, int n3, double *table, int ntable) {
01852 
01853   int n1over2;
01854 
01855   n1over2 = n1 / 2;
01856   return pubz3di(&n1over2,&n2,&n3,table,&ntable);
01857 
01858 } /* pubd3di */
01859 
01860 /****************************************************************************/
01861 /* real to complex fft */
01862 
01863 int pubdz3d(int isign, int n1, int n2,
01864    int n3, double *w, int ld1, int ld2, double
01865    *table, int ntable, double *work) {
01866 
01867   int n1over2, ld1over2, rval;
01868   int i, j, j2, k, k2, i1, i2, imax;
01869   double *data, *data2;
01870   double TwoPiOverN, tmp1r, tmp1i, tmp2r, tmp2i;
01871 
01872   /* complex transform */
01873   n1over2 = n1 / 2;
01874   ld1over2 = ld1 / 2;
01875   rval = pubz3d(&isign, &n1over2, &n2, &n3, (doublecomplex*) w,
01876          &ld1over2, &ld2, table, &ntable, (doublecomplex*) work);
01877 
01878   /* rearrange data */
01879   TwoPiOverN = isign * 2.0 * M_PI / n1;
01880   imax = n1/4+1;
01881   for ( i=0; i<imax; ++i ) {
01882     work[2*i] = cos(i * TwoPiOverN);
01883     work[2*i+1] = sin(i * TwoPiOverN);
01884   }
01885   for ( k=0; k<n3; ++k ) {
01886     for ( j=0; j<n2; ++j ) {
01887       data = w + ld1*(ld2*k + j);
01888       data[n1] = data[0];
01889       data[n1+1] = data[1];
01890     }
01891   }
01892   for ( k=0; k<n3; ++k ) {
01893     k2 = k?(n3-k):0;
01894     for ( j=0; j<n2; ++j ) {
01895       j2 = j?(n2-j):0;
01896       data = w + ld1*(ld2*k + j);
01897       data2 = w + ld1*(ld2*k2 + j2);
01898       imax = n1/4;
01899       if ( (n1/2) & 1 ) imax += 1;
01900       else {
01901         if ( (2*j<n2) || (2*j==n2 && 2*k<=n3) ) imax +=1;
01902         if ( j==0 && 2*k>n3 ) imax -=1;
01903       }
01904       for ( i=0; i<imax; ++i ) {
01905         i1 = 2*i;  i2 = n1-i1;
01906         tmp1r = data[i1] - data2[i2];
01907         tmp1i = data[i1+1] + data2[i2+1];
01908         tmp2r = tmp1r * work[i1+1] + tmp1i * work[i1];
01909         tmp2i = tmp1i * work[i1+1] - tmp1r * work[i1];
01910         tmp1r = data[i1] + data2[i2];
01911         tmp1i = data[i1+1] - data2[i2+1];
01912         data[i1] = 0.5 * ( tmp1r + tmp2r );
01913         data[i1+1] = 0.5 * ( tmp1i + tmp2i );
01914         data2[i2] = 0.5 * ( tmp1r - tmp2r );
01915         data2[i2+1] = 0.5 * ( tmp2i - tmp1i );
01916       }
01917     }
01918   }
01919 
01920   return rval;
01921 
01922 } /* pubdz3d */
01923 
01924 /****************************************************************************/
01925 /* complex to real fft */
01926 
01927 int pubzd3d(int isign, int n1, int n2,
01928    int n3, double *w, int ld1, int ld2, double
01929    *table, int ntable, double *work) {
01930 
01931   int n1over2, ld1over2;
01932   int i, j, j2, k, k2, i1, i2, imax;
01933   double *data, *data2;
01934   double TwoPiOverN, tmp1r, tmp1i, tmp2r, tmp2i;
01935 
01936   /* rearrange data */
01937   TwoPiOverN = isign * 2.0 * M_PI / n1;
01938   imax = n1/4+1;
01939   for ( i=0; i<imax; ++i ) {
01940     work[2*i] = -cos(i * TwoPiOverN);
01941     work[2*i+1] = -sin(i * TwoPiOverN);
01942   }
01943   for ( k=0; k<n3; ++k ) {
01944     k2 = k?(n3-k):0;
01945     for ( j=0; j<n2; ++j ) {
01946       j2 = j?(n2-j):0;
01947       data = w + ld1*(ld2*k + j);
01948       data2 = w + ld1*(ld2*k2 + j2);
01949       imax = n1/4;
01950       if ( (n1/2) & 1 ) imax += 1;
01951       else {
01952         if ( (2*j<n2) || (2*j==n2 && 2*k<=n3) ) imax +=1;
01953         if ( j==0 && 2*k>n3 ) imax -=1;
01954       }
01955       for ( i=0; i<imax; ++i ) {
01956         i1 = 2*i;  i2 = n1-i1;
01957         tmp1r = data[i1] - data2[i2];
01958         tmp1i = data[i1+1] + data2[i2+1];
01959         tmp2r = tmp1r * work[i1+1] + tmp1i * work[i1];
01960         tmp2i = tmp1i * work[i1+1] - tmp1r * work[i1];
01961         tmp1r = data[i1] + data2[i2];
01962         tmp1i = data[i1+1] - data2[i2+1];
01963         data[i1] = tmp1r + tmp2r;
01964         data[i1+1] = tmp1i + tmp2i;
01965         data2[i2] = tmp1r - tmp2r;
01966         data2[i2+1] = tmp2i - tmp1i;
01967       }
01968     }
01969   }
01970 
01971   /* complex transform */
01972   n1over2 = n1 / 2;
01973   ld1over2 = ld1 / 2;
01974   return pubz3d(&isign, &n1over2, &n2, &n3, (doublecomplex*) w,
01975          &ld1over2, &ld2, table, &ntable, (doublecomplex*) work);
01976 
01977 } /* pubzd3d */
01978 
01979 /****************************************************************************/
01980 

Generated on Mon Nov 20 01:17:14 2017 for NAMD by  doxygen 1.4.7