src/gromacs/fft/fft.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * Gromacs 4.0                         Copyright (c) 1991-2003
   5  * Erik Lindahl, David van der Spoel, University of Groningen.
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation; either version 2
  10  * of the License, or (at your option) any later version.
  11  *
  12  * To help us fund GROMACS development, we humbly ask that you cite
  13  * the research papers on the package. Check out http://www.gromacs.org
  14  *
  15  * And Hey:
  16  * Gnomes, ROck Monsters And Chili Sauce
  17  */
  18 #ifdef HAVE_CONFIG_H
  19 #include <config.h>
  20 #endif
  21
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24 #include <string.h>
  25 #include <errno.h>
  26
  27 #include "types/simple.h"
  28 #include "gmxcomplex.h"
  29 #include "gmx_fatal.h"
  30 #include "gromacs/fft/fft.h"
  31
  32
  33 /* This file contains common fft utility functions, but not
  34  * the actual transform implementations. Check the
  35  * files like fft_fftw3.c or fft_mkl.c for that.
  36  */
  37
  38 #ifndef GMX_FFT_FFTW3
  39
  40 struct gmx_many_fft {
  41     int       howmany;
  42     int       dist;
  43     gmx_fft_t fft;
  44 };
  45
  46 typedef struct gmx_many_fft* gmx_many_fft_t;
  47
  48 int
  49 gmx_fft_init_many_1d(gmx_fft_t *        pfft,
  50                      int                nx,
  51                      int                howmany,
  52                      gmx_fft_flag       flags)
  53 {
  54     gmx_many_fft_t fft;
  55     if (pfft == NULL)
  56     {
  57         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
  58         return EINVAL;
  59     }
  60     *pfft = NULL;
  61
  62     if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
  63     {
  64         return ENOMEM;
  65     }
  66
  67     gmx_fft_init_1d(&fft->fft, nx, flags);
  68     fft->howmany = howmany;
  69     fft->dist    = 2*nx;
  70
  71     *pfft = (gmx_fft_t)fft;
  72     return 0;
  73 }
  74
  75 int
  76 gmx_fft_init_many_1d_real(gmx_fft_t *        pfft,
  77                           int                nx,
  78                           int                howmany,
  79                           gmx_fft_flag       flags)
  80 {
  81     gmx_many_fft_t fft;
  82     if (pfft == NULL)
  83     {
  84         gmx_fatal(FARGS, "Invalid opaque FFT datatype pointer.");
  85         return EINVAL;
  86     }
  87     *pfft = NULL;
  88
  89     if ( (fft = (gmx_many_fft_t)malloc(sizeof(struct gmx_many_fft))) == NULL)
  90     {
  91         return ENOMEM;
  92     }
  93
  94     gmx_fft_init_1d_real(&fft->fft, nx, flags);
  95     fft->howmany = howmany;
  96     fft->dist    = 2*(nx/2+1);
  97
  98     *pfft = (gmx_fft_t)fft;
  99     return 0;
 100 }
 101
 102 int
 103 gmx_fft_many_1d     (gmx_fft_t                  fft,
 104                      enum gmx_fft_direction     dir,
 105                      void *                     in_data,
 106                      void *                     out_data)
 107 {
 108     gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
 109     int            i, ret;
 110     for (i = 0; i < mfft->howmany; i++)
 111     {
 112         ret = gmx_fft_1d(mfft->fft, dir, in_data, out_data);
 113         if (ret != 0)
 114         {
 115             return ret;
 116         }
 117         in_data  = (real*)in_data+mfft->dist;
 118         out_data = (real*)out_data+mfft->dist;
 119     }
 120     return 0;
 121 }
 122
 123 int
 124 gmx_fft_many_1d_real     (gmx_fft_t                  fft,
 125                           enum gmx_fft_direction     dir,
 126                           void *                     in_data,
 127                           void *                     out_data)
 128 {
 129     gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
 130     int            i, ret;
 131     for (i = 0; i < mfft->howmany; i++)
 132     {
 133         ret = gmx_fft_1d_real(mfft->fft, dir, in_data, out_data);
 134         if (ret != 0)
 135         {
 136             return ret;
 137         }
 138         in_data  = (real*)in_data+mfft->dist;
 139         out_data = (real*)out_data+mfft->dist;
 140     }
 141     return 0;
 142 }
 143
 144
 145 void
 146 gmx_many_fft_destroy(gmx_fft_t    fft)
 147 {
 148     gmx_many_fft_t mfft = (gmx_many_fft_t)fft;
 149     if (mfft != NULL)
 150     {
 151         if (mfft->fft != NULL)
 152         {
 153             gmx_fft_destroy(mfft->fft);
 154         }
 155         free(mfft);
 156     }
 157 }
 158
 159 #endif //not GMX_FFT_FFTW3
 160
 161 int gmx_fft_transpose_2d(t_complex *          in_data,
 162                          t_complex *          out_data,
 163                          int                  nx,
 164                          int                  ny)
 165 {
 166     int        i, j, k, im, n, ncount, done1, done2;
 167     int        i1, i1c, i2, i2c, kmi, max;
 168
 169     t_complex  tmp1, tmp2, tmp3;
 170     t_complex *data;
 171
 172     /* Use 500 bytes on stack to indicate moves.
 173      * This is just for optimization, it does not limit any dimensions.
 174      */
 175     char          move[500];
 176     int           nmove = 500;
 177
 178     if (nx < 2 || ny < 2)
 179     {
 180         if (in_data != out_data)
 181         {
 182             memcpy(out_data, in_data, sizeof(t_complex)*nx*ny);
 183         }
 184         return 0;
 185     }
 186
 187     /* Out-of-place transposes are easy */
 188     if (in_data != out_data)
 189     {
 190         for (i = 0; i < nx; i++)
 191         {
 192             for (j = 0; j < ny; j++)
 193             {
 194                 out_data[j*nx+i].re = in_data[i*ny+j].re;
 195                 out_data[j*nx+i].im = in_data[i*ny+j].im;
 196             }
 197         }
 198         return 0;
 199     }
 200
 201     /* In-place transform. in_data=out_data=data */
 202     data = in_data;
 203
 204     if (nx == ny)
 205     {
 206         /* trivial case, just swap elements */
 207         for (i = 0; i < nx; i++)
 208         {
 209             for (j = i+1; j < nx; j++)
 210             {
 211                 tmp1.re         = data[i*nx+j].re;
 212                 tmp1.im         = data[i*nx+j].im;
 213                 data[i*nx+j].re = data[j*nx+i].re;
 214                 data[i*nx+j].im = data[j*nx+i].im;
 215                 data[j*nx+i].re = tmp1.re;
 216                 data[j*nx+i].im = tmp1.im;
 217             }
 218         }
 219         return 0;
 220     }
 221
 222     for (i = 0; i < nmove; i++)
 223     {
 224         move[i] = 0;
 225     }
 226
 227     ncount = 2;
 228
 229     if (nx > 2 && ny > 2)
 230     {
 231         i = nx-1;
 232         j = ny-1;
 233         do
 234         {
 235             k = i % j;
 236             i = j;
 237             j = k;
 238         }
 239         while (k);
 240         ncount += i-1;
 241     }
 242
 243     n  = nx*ny;
 244     k  = n - 1;
 245     i  = 1;
 246     im = ny;
 247
 248     done1 = 0;
 249     do
 250     {
 251         i1      = i;
 252         kmi     = k-i;
 253         tmp1.re = data[i1].re;
 254         tmp1.im = data[i1].im;
 255         i1c     = kmi;
 256         tmp2.re = data[i1c].re;
 257         tmp2.im = data[i1c].im;
 258
 259         done2 = 0;
 260         do
 261         {
 262             i2  = ny*i1-k*(i1/nx);
 263             i2c = k-i2;
 264             if (i1 < nmove)
 265             {
 266                 move[i1] = 1;
 267             }
 268             if (i1c < nmove)
 269             {
 270                 move[i1c] = 1;
 271             }
 272             ncount += 2;
 273             if (i2 == i)
 274             {
 275                 done2 = 1;
 276             }
 277             else if (i2 == kmi)
 278             {
 279                 tmp3.re = tmp1.re;
 280                 tmp3.im = tmp1.im;
 281                 tmp1.re = tmp2.re;
 282                 tmp1.im = tmp2.im;
 283                 tmp2.re = tmp3.re;
 284                 tmp2.im = tmp3.im;
 285                 done2   = 1;
 286             }
 287             else
 288             {
 289                 data[i1].re  = data[i2].re;
 290                 data[i1].im  = data[i2].im;
 291                 data[i1c].re = data[i2c].re;
 292                 data[i1c].im = data[i2c].im;
 293                 i1           = i2;
 294                 i1c          = i2c;
 295             }
 296         }
 297         while (!done2);
 298
 299         data[i1].re  = tmp1.re;
 300         data[i1].im  = tmp1.im;
 301         data[i1c].re = tmp2.re;
 302         data[i1c].im = tmp2.im;
 303
 304         if (ncount >= n)
 305         {
 306             done1 = 1;
 307         }
 308         else
 309         {
 310             done2 = 0;
 311             do
 312             {
 313                 max = k-i;
 314                 i++;
 315                 im += ny;
 316                 if (im > k)
 317                 {
 318                     im -= k;
 319                 }
 320                 i2 = im;
 321                 if (i != i2)
 322                 {
 323                     if (i >= nmove)
 324                     {
 325                         while (i2 > i && i2 < max)
 326                         {
 327                             i1 = i2;
 328                             i2 = ny*i1-k*(i1/nx);
 329                         }
 330                         if (i2 == i)
 331                         {
 332                             done2 = 1;
 333                         }
 334                     }
 335                     else if (!move[i])
 336                     {
 337                         done2 = 1;
 338                     }
 339                 }
 340             }
 341             while (!done2);
 342         }
 343     }
 344     while (!done1);
 345
 346     return 0;
 347 }