-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
*
*
* Gromacs 4.0 Copyright (c) 1991-2003
*
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org
- *
+ *
* And Hey:
* Gnomes, ROck Monsters And Chili Sauce
*/
#include "gmx_fatal.h"
#include "external/fftpack/fftpack.h"
-/** Contents of the FFTPACK fft datatype.
+/** Contents of the FFTPACK fft datatype.
*
- * FFTPACK only does 1d transforms, so we use a pointers to another fft for
+ * FFTPACK only does 1d transforms, so we use a pointers to another fft for
* the transform in the next dimension.
* Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains
* a pointer to a 1d. The 1d structure has next==NULL.
*/
struct gmx_fft
{
- int ndim; /**< Dimensions, including our subdimensions. */
- int n; /**< Number of points in this dimension. */
- int ifac[15]; /**< 15 bytes needed for cfft and rfft */
- struct gmx_fft *next; /**< Pointer to next dimension, or NULL. */
- real * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */
+ int ndim; /**< Dimensions, including our subdimensions. */
+ int n; /**< Number of points in this dimension. */
+ int ifac[15]; /**< 15 bytes needed for cfft and rfft */
+ struct gmx_fft *next; /**< Pointer to next dimension, or NULL. */
+ real * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */
};
#include <math.h>
int flags)
{
gmx_fft_t fft;
-
- if(pfft==NULL)
+
+ if (pfft == NULL)
{
- gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
+ gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
return EINVAL;
}
*pfft = NULL;
-
- if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
+
+ if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
{
return ENOMEM;
- }
-
+ }
+
fft->next = NULL;
fft->n = nx;
-
+
/* Need 4*n storage for 1D complex FFT */
- if( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
+ if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx))) == NULL)
{
free(fft);
return ENOMEM;
}
- if(fft->n>1)
- fftpack_cffti1(nx,fft->work,fft->ifac);
-
+ if (fft->n > 1)
+ {
+ fftpack_cffti1(nx, fft->work, fft->ifac);
+ }
+
*pfft = fft;
return 0;
};
int flags)
{
gmx_fft_t fft;
-
- if(pfft==NULL)
+
+ if (pfft == NULL)
{
- gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
+ gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
return EINVAL;
}
*pfft = NULL;
- if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
+ if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
{
return ENOMEM;
- }
+ }
fft->next = NULL;
fft->n = nx;
/* Need 2*n storage for 1D real FFT */
- if((fft->work = (real *)malloc(sizeof(real)*(2*nx)))==NULL)
+ if ((fft->work = (real *)malloc(sizeof(real)*(2*nx))) == NULL)
{
free(fft);
return ENOMEM;
- }
-
- if(fft->n>1)
- fftpack_rffti1(nx,fft->work,fft->ifac);
-
+ }
+
+ if (fft->n > 1)
+ {
+ fftpack_rffti1(nx, fft->work, fft->ifac);
+ }
+
*pfft = fft;
return 0;
}
gmx_fft_t fft;
int nyc = (ny/2 + 1);
int rc;
-
- if(pfft==NULL)
+
+ if (pfft == NULL)
{
- gmx_fatal(FARGS,"Invalid FFT opaque type pointer.");
+ gmx_fatal(FARGS, "Invalid FFT opaque type pointer.");
return EINVAL;
}
*pfft = NULL;
-
+
/* Create the X transform */
- if( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
+ if ( (fft = (struct gmx_fft *)malloc(sizeof(struct gmx_fft))) == NULL)
{
return ENOMEM;
- }
-
+ }
+
fft->n = nx;
-
+
/* Need 4*nx storage for 1D complex FFT, and another
* 2*nx*nyc elements for complex-to-real storage in our high-level routine.
*/
- if( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
+ if ( (fft->work = (real *)malloc(sizeof(real)*(4*nx+2*nx*nyc))) == NULL)
{
free(fft);
return ENOMEM;
}
- fftpack_cffti1(nx,fft->work,fft->ifac);
-
+ fftpack_cffti1(nx, fft->work, fft->ifac);
+
/* Create real Y transform as a link from X */
- if( (rc=gmx_fft_init_1d_real(&(fft->next),ny,flags)) != 0)
+ if ( (rc = gmx_fft_init_1d_real(&(fft->next), ny, flags)) != 0)
{
free(fft);
return rc;
void * in_data,
void * out_data)
{
- int i,n;
- real * p1;
- real * p2;
+ int i, n;
+ real * p1;
+ real * p2;
- n=fft->n;
+ n = fft->n;
- if(n==1)
+ if (n == 1)
{
- p1 = (real *)in_data;
- p2 = (real *)out_data;
+ p1 = (real *)in_data;
+ p2 = (real *)out_data;
p2[0] = p1[0];
p2[1] = p1[1];
}
-
+
/* FFTPACK only does in-place transforms, so emulate out-of-place
* by copying data to the output array first.
*/
- if( in_data != out_data )
+ if (in_data != out_data)
{
p1 = (real *)in_data;
p2 = (real *)out_data;
-
+
/* n complex = 2*n real elements */
- for(i=0;i<2*n;i++)
+ for (i = 0; i < 2*n; i++)
{
p2[i] = p1[i];
}
}
-
+
/* Elements 0 .. 2*n-1 in work are used for ffac values,
* Elements 2*n .. 4*n-1 are internal FFTPACK work space.
*/
-
- if(dir == GMX_FFT_FORWARD)
+
+ if (dir == GMX_FFT_FORWARD)
{
- fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, -1);
+ fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, -1);
}
- else if(dir == GMX_FFT_BACKWARD)
+ else if (dir == GMX_FFT_BACKWARD)
{
- fftpack_cfftf1(n,(real *)out_data,fft->work+2*n,fft->work,fft->ifac, 1);
+ fftpack_cfftf1(n, (real *)out_data, fft->work+2*n, fft->work, fft->ifac, 1);
}
else
{
- gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
+ gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
return EINVAL;
- }
+ }
return 0;
}
-int
+int
gmx_fft_1d_real (gmx_fft_t fft,
enum gmx_fft_direction dir,
void * in_data,
void * out_data)
{
- int i,n;
- real * p1;
- real * p2;
+ int i, n;
+ real * p1;
+ real * p2;
n = fft->n;
-
- if(n==1)
+
+ if (n == 1)
{
- p1 = (real *)in_data;
- p2 = (real *)out_data;
+ p1 = (real *)in_data;
+ p2 = (real *)out_data;
p2[0] = p1[0];
- if(dir == GMX_FFT_REAL_TO_COMPLEX)
+ if (dir == GMX_FFT_REAL_TO_COMPLEX)
+ {
p2[1] = 0.0;
+ }
}
-
- if(dir == GMX_FFT_REAL_TO_COMPLEX)
+
+ if (dir == GMX_FFT_REAL_TO_COMPLEX)
{
/* FFTPACK only does in-place transforms, so emulate out-of-place
* by copying data to the output array first. This works fine, since
* the complex array must be larger than the real.
*/
- if( in_data != out_data )
+ if (in_data != out_data)
{
p1 = (real *)in_data;
p2 = (real *)out_data;
-
- for(i=0;i<2*(n/2+1);i++)
+
+ for (i = 0; i < 2*(n/2+1); i++)
{
p2[i] = p1[i];
}
/* Elements 0 .. n-1 in work are used for ffac values,
* Elements n .. 2*n-1 are internal FFTPACK work space.
*/
- fftpack_rfftf1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
+ fftpack_rfftf1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
/*
- * FFTPACK has a slightly more compact storage than we, time to
- * convert it: ove most of the array one step up to make room for
- * zero imaginary parts.
+ * FFTPACK has a slightly more compact storage than we, time to
+ * convert it: ove most of the array one step up to make room for
+ * zero imaginary parts.
*/
p2 = (real *)out_data;
- for(i=n-1;i>0;i--)
+ for (i = n-1; i > 0; i--)
{
p2[i+1] = p2[i];
}
/* imaginary zero freq. */
- p2[1] = 0;
-
+ p2[1] = 0;
+
/* Is n even? */
- if( (n & 0x1) == 0 )
+ if ( (n & 0x1) == 0)
{
p2[n+1] = 0;
}
-
+
}
- else if(dir == GMX_FFT_COMPLEX_TO_REAL)
+ else if (dir == GMX_FFT_COMPLEX_TO_REAL)
{
/* FFTPACK only does in-place transforms, and we cannot just copy
* input to output first here since our real array is smaller than
- * the complex one. However, since the FFTPACK complex storage format
+ * the complex one. However, since the FFTPACK complex storage format
* is more compact than ours (2 reals) it will fit, so compact it
* and copy on-the-fly to the output array.
*/
p2 = (real *)out_data;
p2[0] = p1[0];
- for(i=1;i<n;i++)
+ for (i = 1; i < n; i++)
{
p2[i] = p1[i+1];
}
- fftpack_rfftb1(n,(real *)out_data,fft->work+n,fft->work,fft->ifac);
- }
+ fftpack_rfftb1(n, (real *)out_data, fft->work+n, fft->work, fft->ifac);
+ }
else
{
- gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
+ gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
return EINVAL;
- }
-
+ }
+
return 0;
}
-int
+int
gmx_fft_2d_real (gmx_fft_t fft,
enum gmx_fft_direction dir,
void * in_data,
void * out_data)
{
- int i,j,nx,ny,nyc;
- t_complex * data;
- real * work;
- real * p1;
- real * p2;
-
- nx=fft->n;
- ny=fft->next->n;
+ int i, j, nx, ny, nyc;
+ t_complex * data;
+ real * work;
+ real * p1;
+ real * p2;
+
+ nx = fft->n;
+ ny = fft->next->n;
/* Number of complex elements in y direction */
- nyc=(ny/2+1);
-
+ nyc = (ny/2+1);
+
work = fft->work+4*nx;
-
- if(dir==GMX_FFT_REAL_TO_COMPLEX)
+
+ if (dir == GMX_FFT_REAL_TO_COMPLEX)
{
/* If we are doing an in-place transform the 2D array is already
* properly padded by the user, and we are all set.
* does in-place FFTs internally, so we need to start by copying
* data from the input to the padded (larger) output array.
*/
- if( in_data != out_data )
+ if (in_data != out_data)
{
p1 = (real *)in_data;
p2 = (real *)out_data;
-
- for(i=0;i<nx;i++)
+
+ for (i = 0; i < nx; i++)
{
- for(j=0;j<ny;j++)
+ for (j = 0; j < ny; j++)
{
p2[i*nyc*2+j] = p1[i*ny+j];
}
data = (t_complex *)out_data;
/* y real-to-complex FFTs */
- for(i=0;i<nx;i++)
+ for (i = 0; i < nx; i++)
{
- gmx_fft_1d_real(fft->next,GMX_FFT_REAL_TO_COMPLEX,data+i*nyc,data+i*nyc);
+ gmx_fft_1d_real(fft->next, GMX_FFT_REAL_TO_COMPLEX, data+i*nyc, data+i*nyc);
}
-
+
/* Transform to get X data in place */
- gmx_fft_transpose_2d(data,data,nx,nyc);
-
+ gmx_fft_transpose_2d(data, data, nx, nyc);
+
/* Complex-to-complex X FFTs */
- for(i=0;i<nyc;i++)
+ for (i = 0; i < nyc; i++)
{
- gmx_fft_1d(fft,GMX_FFT_FORWARD,data+i*nx,data+i*nx);
+ gmx_fft_1d(fft, GMX_FFT_FORWARD, data+i*nx, data+i*nx);
}
-
+
/* Transpose back */
- gmx_fft_transpose_2d(data,data,nyc,nx);
-
+ gmx_fft_transpose_2d(data, data, nyc, nx);
+
}
- else if(dir==GMX_FFT_COMPLEX_TO_REAL)
+ else if (dir == GMX_FFT_COMPLEX_TO_REAL)
{
/* An in-place complex-to-real transform is straightforward,
* since the output array must be large enough for the padding to fit.
* In this case there's nothing to do but employing temporary work data,
* starting at work+4*nx and using nx*nyc*2 elements.
*/
- if(in_data != out_data)
+ if (in_data != out_data)
{
- memcpy(work,in_data,sizeof(t_complex)*nx*nyc);
+ memcpy(work, in_data, sizeof(t_complex)*nx*nyc);
data = (t_complex *)work;
}
else
}
/* Transpose to get X arrays */
- gmx_fft_transpose_2d(data,data,nx,nyc);
-
+ gmx_fft_transpose_2d(data, data, nx, nyc);
+
/* Do X iFFTs */
- for(i=0;i<nyc;i++)
+ for (i = 0; i < nyc; i++)
{
- gmx_fft_1d(fft,GMX_FFT_BACKWARD,data+i*nx,data+i*nx);
+ gmx_fft_1d(fft, GMX_FFT_BACKWARD, data+i*nx, data+i*nx);
}
-
+
/* Transpose to get Y arrays */
- gmx_fft_transpose_2d(data,data,nyc,nx);
-
+ gmx_fft_transpose_2d(data, data, nyc, nx);
+
/* Do Y iFFTs */
- for(i=0;i<nx;i++)
+ for (i = 0; i < nx; i++)
{
- gmx_fft_1d_real(fft->next,GMX_FFT_COMPLEX_TO_REAL,data+i*nyc,data+i*nyc);
+ gmx_fft_1d_real(fft->next, GMX_FFT_COMPLEX_TO_REAL, data+i*nyc, data+i*nyc);
}
-
- if( in_data != out_data )
+
+ if (in_data != out_data)
{
/* Output (pointed to by data) is now in padded format.
* Pack it into out_data if we were doing an out-of-place transform.
*/
p1 = (real *)data;
p2 = (real *)out_data;
-
- for(i=0;i<nx;i++)
+
+ for (i = 0; i < nx; i++)
{
- for(j=0;j<ny;j++)
+ for (j = 0; j < ny; j++)
{
p2[i*ny+j] = p1[i*nyc*2+j];
}
}
else
{
- gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
+ gmx_fatal(FARGS, "FFT plan mismatch - bad plan or direction.");
return EINVAL;
- }
-
+ }
+
return 0;
}
void
gmx_fft_destroy(gmx_fft_t fft)
{
- if(fft != NULL)
+ if (fft != NULL)
{
free(fft->work);
- if(fft->next != NULL)
+ if (fft->next != NULL)
+ {
gmx_fft_destroy(fft->next);
+ }
free(fft);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff