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