pub3dfft.h File Reference

Go to the source code of this file.

Classes

struct  doublecomplex

Functions

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)


Function Documentation

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 Mon Nov 20 01:17:16 2017 for NAMD by  doxygen 1.4.7