NAMD
Classes | Functions
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().

1851  {
1852 
1853  int n1over2;
1854 
1855  n1over2 = n1 / 2;
1856  return pubz3di(&n1over2,&n2,&n3,table,&ntable);
1857 
1858 } /* pubd3di */
int pubz3di(int *n1, int *n2, int *n3, double *table, int *ntable)
Definition: pub3dfft.C:1706
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 M_PI, and pubz3d().

1865  {
1866 
1867  int n1over2, ld1over2, rval;
1868  int i, j, j2, k, k2, i1, i2, imax;
1869  double *data, *data2;
1870  double TwoPiOverN, tmp1r, tmp1i, tmp2r, tmp2i;
1871 
1872  /* complex transform */
1873  n1over2 = n1 / 2;
1874  ld1over2 = ld1 / 2;
1875  rval = pubz3d(&isign, &n1over2, &n2, &n3, (doublecomplex*) w,
1876  &ld1over2, &ld2, table, &ntable, (doublecomplex*) work);
1877 
1878  /* rearrange data */
1879  TwoPiOverN = isign * 2.0 * M_PI / n1;
1880  imax = n1/4+1;
1881  for ( i=0; i<imax; ++i ) {
1882  work[2*i] = cos(i * TwoPiOverN);
1883  work[2*i+1] = sin(i * TwoPiOverN);
1884  }
1885  for ( k=0; k<n3; ++k ) {
1886  for ( j=0; j<n2; ++j ) {
1887  data = w + ld1*(ld2*k + j);
1888  data[n1] = data[0];
1889  data[n1+1] = data[1];
1890  }
1891  }
1892  for ( k=0; k<n3; ++k ) {
1893  k2 = k?(n3-k):0;
1894  for ( j=0; j<n2; ++j ) {
1895  j2 = j?(n2-j):0;
1896  data = w + ld1*(ld2*k + j);
1897  data2 = w + ld1*(ld2*k2 + j2);
1898  imax = n1/4;
1899  if ( (n1/2) & 1 ) imax += 1;
1900  else {
1901  if ( (2*j<n2) || (2*j==n2 && 2*k<=n3) ) imax +=1;
1902  if ( j==0 && 2*k>n3 ) imax -=1;
1903  }
1904  for ( i=0; i<imax; ++i ) {
1905  i1 = 2*i; i2 = n1-i1;
1906  tmp1r = data[i1] - data2[i2];
1907  tmp1i = data[i1+1] + data2[i2+1];
1908  tmp2r = tmp1r * work[i1+1] + tmp1i * work[i1];
1909  tmp2i = tmp1i * work[i1+1] - tmp1r * work[i1];
1910  tmp1r = data[i1] + data2[i2];
1911  tmp1i = data[i1+1] - data2[i2+1];
1912  data[i1] = 0.5 * ( tmp1r + tmp2r );
1913  data[i1+1] = 0.5 * ( tmp1i + tmp2i );
1914  data2[i2] = 0.5 * ( tmp1r - tmp2r );
1915  data2[i2+1] = 0.5 * ( tmp2i - tmp1i );
1916  }
1917  }
1918  }
1919 
1920  return rval;
1921 
1922 } /* pubdz3d */
#define M_PI
Definition: pub3dfft.C:15
int pubz3d(int *isign, int *n1, int *n2, int *n3, doublecomplex *w, int *ld1, int *ld2, double *table, int *ntable, doublecomplex *work)
Definition: pub3dfft.C:1729
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, and doublecomplex::r.

Referenced by pubdz3d(), and pubzd3d().

1732 {
1733  /* System generated locals */
1734  int w_dim1, w_dim2, w_offset, table_dim1, table_offset, i_1, i_2, i_3,
1735  i_4, i_5;
1736 
1737  /* Local variables */
1738  static int i, j, k;
1739 
1740  /* Parameter adjustments */
1741  w_dim1 = *ld1;
1742  w_dim2 = *ld2;
1743  w_offset = w_dim1 * (w_dim2 + 1) + 1;
1744  w -= w_offset;
1745  table_dim1 = *ntable;
1746  table_offset = table_dim1 + 1;
1747  table -= table_offset;
1748  --work;
1749 
1750  /* Function Body */
1751 /* ntable should be 4*max(n1,n2,n3) +15 */
1752 /* nwork should be max(n1,n2,n3) */
1753 
1754 /* transform along X first ... */
1755 
1756  i_1 = *n3;
1757  for (k = 1; k <= i_1; ++k) {
1758  i_2 = *n2;
1759  for (j = 1; j <= i_2; ++j) {
1760  i_3 = *n1;
1761  for (i = 1; i <= i_3; ++i) {
1762  i_4 = i;
1763  i_5 = i + (j + k * w_dim2) * w_dim1;
1764  work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
1765 /* L70: */
1766  }
1767  if (*isign == -1) {
1768  cfftf(n1, (double*) &work[1], &table[table_dim1 + 1]);
1769  }
1770  if (*isign == 1) {
1771  cfftb(n1, (double*) &work[1], &table[table_dim1 + 1]);
1772  }
1773  i_3 = *n1;
1774  for (i = 1; i <= i_3; ++i) {
1775  i_4 = i + (j + k * w_dim2) * w_dim1;
1776  i_5 = i;
1777  w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
1778 /* L80: */
1779  }
1780 /* L90: */
1781  }
1782 /* L100: */
1783  }
1784 
1785 /* transform along Y then ... */
1786 
1787  i_1 = *n3;
1788  for (k = 1; k <= i_1; ++k) {
1789  i_2 = *n1;
1790  for (i = 1; i <= i_2; ++i) {
1791  i_3 = *n2;
1792  for (j = 1; j <= i_3; ++j) {
1793  i_4 = j;
1794  i_5 = i + (j + k * w_dim2) * w_dim1;
1795  work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
1796 /* L170: */
1797  }
1798  if (*isign == -1) {
1799  cfftf(n2, (double*) &work[1], &table[(table_dim1 << 1) + 1]);
1800  }
1801  if (*isign == 1) {
1802  cfftb(n2, (double*) &work[1], &table[(table_dim1 << 1) + 1]);
1803  }
1804  i_3 = *n2;
1805  for (j = 1; j <= i_3; ++j) {
1806  i_4 = i + (j + k * w_dim2) * w_dim1;
1807  i_5 = j;
1808  w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
1809 /* L180: */
1810  }
1811 /* L190: */
1812  }
1813 /* L200: */
1814  }
1815 
1816 /* transform along Z finally ... */
1817 
1818  i_1 = *n1;
1819  for (i = 1; i <= i_1; ++i) {
1820  i_2 = *n2;
1821  for (j = 1; j <= i_2; ++j) {
1822  i_3 = *n3;
1823  for (k = 1; k <= i_3; ++k) {
1824  i_4 = k;
1825  i_5 = i + (j + k * w_dim2) * w_dim1;
1826  work[i_4].r = w[i_5].r, work[i_4].i = w[i_5].i;
1827 /* L270: */
1828  }
1829  if (*isign == -1) {
1830  cfftf(n3, (double*) &work[1], &table[table_dim1 * 3 + 1]);
1831  }
1832  if (*isign == 1) {
1833  cfftb(n3, (double*) &work[1], &table[table_dim1 * 3 + 1]);
1834  }
1835  i_3 = *n3;
1836  for (k = 1; k <= i_3; ++k) {
1837  i_4 = i + (j + k * w_dim2) * w_dim1;
1838  i_5 = k;
1839  w[i_4].r = work[i_5].r, w[i_4].i = work[i_5].i;
1840 /* L280: */
1841  }
1842 /* L290: */
1843  }
1844 /* L300: */
1845  }
1846  return 0;
1847 } /* pubz3d_ */
int cfftb(int *n, double *c, double *wsave)
Definition: pub3dfft.C:1401
double r
Definition: pub3dfft.h:10
int cfftf(int *n, double *c, double *wsave)
Definition: pub3dfft.C:1544
double i
Definition: pub3dfft.h:10
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().

1708 {
1709  /* System generated locals */
1710  int table_dim1, table_offset;
1711 
1712  /* Local variables */
1713 
1714  /* Parameter adjustments */
1715  table_dim1 = *ntable;
1716  table_offset = table_dim1 + 1;
1717  table -= table_offset;
1718 
1719  /* Function Body */
1720 /* ntable should be 4*max(n1,n2,n3) +15 */
1721  cffti(n1, &table[table_dim1 + 1]);
1722  cffti(n2, &table[(table_dim1 << 1) + 1]);
1723  cffti(n3, &table[table_dim1 * 3 + 1]);
1724  return 0;
1725 } /* pubz3di_ */
int cffti(int *n, double *wsave)
Definition: pub3dfft.C:1676
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 M_PI, and pubz3d().

1929  {
1930 
1931  int n1over2, ld1over2;
1932  int i, j, j2, k, k2, i1, i2, imax;
1933  double *data, *data2;
1934  double TwoPiOverN, tmp1r, tmp1i, tmp2r, tmp2i;
1935 
1936  /* rearrange data */
1937  TwoPiOverN = isign * 2.0 * M_PI / n1;
1938  imax = n1/4+1;
1939  for ( i=0; i<imax; ++i ) {
1940  work[2*i] = -cos(i * TwoPiOverN);
1941  work[2*i+1] = -sin(i * TwoPiOverN);
1942  }
1943  for ( k=0; k<n3; ++k ) {
1944  k2 = k?(n3-k):0;
1945  for ( j=0; j<n2; ++j ) {
1946  j2 = j?(n2-j):0;
1947  data = w + ld1*(ld2*k + j);
1948  data2 = w + ld1*(ld2*k2 + j2);
1949  imax = n1/4;
1950  if ( (n1/2) & 1 ) imax += 1;
1951  else {
1952  if ( (2*j<n2) || (2*j==n2 && 2*k<=n3) ) imax +=1;
1953  if ( j==0 && 2*k>n3 ) imax -=1;
1954  }
1955  for ( i=0; i<imax; ++i ) {
1956  i1 = 2*i; i2 = n1-i1;
1957  tmp1r = data[i1] - data2[i2];
1958  tmp1i = data[i1+1] + data2[i2+1];
1959  tmp2r = tmp1r * work[i1+1] + tmp1i * work[i1];
1960  tmp2i = tmp1i * work[i1+1] - tmp1r * work[i1];
1961  tmp1r = data[i1] + data2[i2];
1962  tmp1i = data[i1+1] - data2[i2+1];
1963  data[i1] = tmp1r + tmp2r;
1964  data[i1+1] = tmp1i + tmp2i;
1965  data2[i2] = tmp1r - tmp2r;
1966  data2[i2+1] = tmp2i - tmp1i;
1967  }
1968  }
1969  }
1970 
1971  /* complex transform */
1972  n1over2 = n1 / 2;
1973  ld1over2 = ld1 / 2;
1974  return pubz3d(&isign, &n1over2, &n2, &n3, (doublecomplex*) w,
1975  &ld1over2, &ld2, table, &ntable, (doublecomplex*) work);
1976 
1977 } /* pubzd3d */
#define M_PI
Definition: pub3dfft.C:15
int pubz3d(int *isign, int *n1, int *n2, int *n3, doublecomplex *w, int *ld1, int *ld2, double *table, int *ntable, doublecomplex *work)
Definition: pub3dfft.C:1729