pub3dfft.C File Reference

#include <math.h>
#include "pub3dfft.h"

Go to the source code of this file.

Defines

#define M_PI   3.14159265358979323846

Functions

int passb (int *nac, int *ido, int *ip, int *l1, int *idl1, double *cc, double *c1, double *c2, double *ch, double *ch2, double *wa)
int passb2 (int *ido, int *l1, double *cc, double *ch, double *wa1)
int passb3 (int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2)
int passb4 (int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3)
int passb5 (int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3, double *wa4)
int passf (int *nac, int *ido, int *ip, int *l1, int *idl1, double *cc, double *c1, double *c2, double *ch, double *ch2, double *wa)
int passf2 (int *ido, int *l1, double *cc, double *ch, double *wa1)
int passf3 (int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2)
int passf4 (int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3)
int passf5 (int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2, double *wa3, double *wa4)
int cfftb1 (int *n, double *c, double *ch, double *wa, int *ifac)
int cfftb (int *n, double *c, double *wsave)
int cfftf1 (int *n, double *c, double *ch, double *wa, int *ifac)
int cfftf (int *n, double *c, double *wsave)
int cffti1 (int *n, double *wa, int *ifac)
int cffti (int *n, double *wsave)
int pubz3di (int *n1, int *n2, int *n3, double *table, int *ntable)
int pubz3d (int *isign, int *n1, int *n2, int *n3, doublecomplex *w, int *ld1, int *ld2, double *table, int *ntable, doublecomplex *work)
int pubd3di (int n1, int n2, int n3, double *table, int ntable)
int pubdz3d (int isign, int n1, int n2, int n3, double *w, int ld1, int ld2, double *table, int ntable, double *work)
int pubzd3d (int isign, int n1, int n2, int n3, double *w, int ld1, int ld2, double *table, int ntable, double *work)


Define Documentation

#define M_PI   3.14159265358979323846

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 15 of file pub3dfft.C.


Function Documentation

int cfftb ( int *  n,
double *  c,
double *  wsave 
)

Definition at line 1401 of file pub3dfft.C.

References cfftb1().

Referenced by pubz3d().

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_ */

int cfftb1 ( int *  n,
double *  c,
double *  ch,
double *  wa,
int *  ifac 
)

Definition at line 1277 of file pub3dfft.C.

References passb(), passb2(), passb3(), passb4(), and passb5().

Referenced by cfftb().

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_ */

int cfftf ( int *  n,
double *  c,
double *  wsave 
)

Definition at line 1544 of file pub3dfft.C.

References cfftf1().

Referenced by pubz3d().

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_ */

int cfftf1 ( int *  n,
double *  c,
double *  ch,
double *  wa,
int *  ifac 
)

Definition at line 1420 of file pub3dfft.C.

References passf(), passf2(), passf3(), passf4(), and passf5().

Referenced by cfftf().

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_ */

int cffti ( int *  n,
double *  wsave 
)

Definition at line 1676 of file pub3dfft.C.

References cffti1().

Referenced by pubz3di().

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_ */

int cffti1 ( int *  n,
double *  wa,
int *  ifac 
)

Definition at line 1563 of file pub3dfft.C.

References j.

Referenced by cffti().

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_ */

int passb ( int *  nac,
int *  ido,
int *  ip,
int *  l1,
int *  idl1,
double *  cc,
double *  c1,
double *  c2,
double *  ch,
double *  ch2,
double *  wa 
)

Definition at line 18 of file pub3dfft.C.

References j.

Referenced by cfftb1().

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_ */

int passb2 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1 
)

Definition at line 257 of file pub3dfft.C.

Referenced by cfftb1().

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_ */

int passb3 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1,
double *  wa2 
)

Definition at line 319 of file pub3dfft.C.

Referenced by cfftb1().

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_ */

int passb4 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1,
double *  wa2,
double *  wa3 
)

Definition at line 408 of file pub3dfft.C.

Referenced by cfftb1().

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_ */

int passb5 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1,
double *  wa2,
double *  wa3,
double *  wa4 
)

Definition at line 512 of file pub3dfft.C.

Referenced by cfftb1().

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_ */

int passf ( int *  nac,
int *  ido,
int *  ip,
int *  l1,
int *  idl1,
double *  cc,
double *  c1,
double *  c2,
double *  ch,
double *  ch2,
double *  wa 
)

Definition at line 647 of file pub3dfft.C.

References j.

Referenced by cfftf1().

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_ */

int passf2 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1 
)

Definition at line 885 of file pub3dfft.C.

Referenced by cfftf1().

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_ */

int passf3 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1,
double *  wa2 
)

Definition at line 946 of file pub3dfft.C.

Referenced by cfftf1().

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_ */

int passf4 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1,
double *  wa2,
double *  wa3 
)

Definition at line 1035 of file pub3dfft.C.

Referenced by cfftf1().

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_ */

int passf5 ( int *  ido,
int *  l1,
double *  cc,
double *  ch,
double *  wa1,
double *  wa2,
double *  wa3,
double *  wa4 
)

Definition at line 1139 of file pub3dfft.C.

Referenced by cfftf1().

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_ */

int pubd3di ( int  n1,
int  n2,
int  n3,
double *  table,
int  ntable 
)

Definition at line 1851 of file pub3dfft.C.

References pubz3di().

01851                                                                {
01852 
01853   int n1over2;
01854 
01855   n1over2 = n1 / 2;
01856   return pubz3di(&n1over2,&n2,&n3,table,&ntable);
01857 
01858 } /* pubd3di */

int pubdz3d ( int  isign,
int  n1,
int  n2,
int  n3,
double *  w,
int  ld1,
int  ld2,
double *  table,
int  ntable,
double *  work 
)

Definition at line 1863 of file pub3dfft.C.

References j, M_PI, and pubz3d().

01865                                      {
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 */

int pubz3d ( int *  isign,
int *  n1,
int *  n2,
int *  n3,
doublecomplex w,
int *  ld1,
int *  ld2,
double *  table,
int *  ntable,
doublecomplex work 
)

Definition at line 1729 of file pub3dfft.C.

References cfftb(), cfftf(), doublecomplex::i, j, and doublecomplex::r.

Referenced by pubdz3d(), and pubzd3d().

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_ */

int pubz3di ( int *  n1,
int *  n2,
int *  n3,
double *  table,
int *  ntable 
)

Definition at line 1706 of file pub3dfft.C.

References cffti().

Referenced by pubd3di().

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_ */

int pubzd3d ( int  isign,
int  n1,
int  n2,
int  n3,
double *  w,
int  ld1,
int  ld2,
double *  table,
int  ntable,
double *  work 
)

Definition at line 1927 of file pub3dfft.C.

References j, M_PI, and pubz3d().

01929                                      {
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 */


Generated on Thu Nov 23 01:17:16 2017 for NAMD by  doxygen 1.4.7