src/mdlib/gmx_fft_acml.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Gromacs 4.0                         Copyright (c) 1991-2003
   5  * David van der Spoel, Erik Lindahl, University of Groningen.
   6  * Copyright (c) 2012, by the GROMACS development team, led by
   7  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   8  * others, as listed in the AUTHORS file in the top-level source
   9  * directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #ifdef GMX_FFT_ACML
  42
  43 #include <errno.h>
  44 #include <stdlib.h>
  45 #include <memory.h>
  46 #include <math.h>
  47 #include <assert.h>
  48
  49 #include "acml.h"
  50
  51
  52 #include "gmx_fft.h"
  53 #include "gmx_fatal.h"
  54
  55 /*      Since c has no support of templates, we get around the type restrictions
  56         with macros.  Our custom macro names obscure the float vs double interfaces
  57         */
  58 #ifdef GMX_DOUBLE
  59 #define ACML_FFT1DX zfft1dx
  60 #define ACML_FFT2DX zfft2dx
  61 #define ACML_FFT3DX zfft3dy
  62 #define ACML_FFT1DMX zfft1mx
  63
  64 #define ACML_RCFFT1D dzfft
  65 #define ACML_CRFFT1D zdfft
  66 #define ACML_RCFFT1DM dzfftm
  67 #define ACML_CRFFT1DM zdfftm
  68
  69 #define acmlComplex doublecomplex
  70 #else
  71
  72 #define ACML_FFT1DX cfft1dx
  73 #define ACML_FFT2DX cfft2dx
  74 #define ACML_FFT3DX cfft3dy
  75 #define ACML_FFT1DMX cfft1mx
  76
  77 #define ACML_RCFFT1D scfft
  78 #define ACML_CRFFT1D csfft
  79 #define ACML_RCFFT1DM scfftm
  80 #define ACML_CRFFT1DM csfftm
  81
  82 #define acmlComplex complex
  83 #endif
  84
  85 void cPrintArray( t_complex* arr, int outer, int mid, int inner )
  86 {
  87         int i, j, k;
  88
  89         printf("\n");
  90         for( i = 0; i < outer; ++i )
  91         {
  92                 for( j = 0; j < mid; ++j )
  93                 {
  94                         for( k = 0; k < inner; ++k )
  95                         {
  96                                 printf("%f,%f ", arr[i*inner*mid+j*inner+k].re, arr[i*inner*mid+j*inner+k].im );
  97                         }
  98                         printf("\n");
  99                 }
 100                 printf("\n");
 101         }
 102
 103         return;
 104 }
 105
 106 void rPrintArray( real* arr, int outer, int mid, int inner )
 107 {
 108         int i, j, k;
 109
 110         printf("\n");
 111         for( i = 0; i < outer; ++i )
 112         {
 113                 for( j = 0; j < mid; ++j )
 114                 {
 115                         for( k = 0; k < inner; ++k )
 116                         {
 117                                 printf("%f ", arr[i*inner*mid+j*inner+k] );
 118                         }
 119                         printf("\n");
 120                 }
 121                 printf("\n");
 122         }
 123
 124         return;
 125 }
 126
 127 /* Contents of the AMD ACML FFT fft datatype.
 128  *
 129  * Note that this is one of several possible implementations of gmx_fft_t.
 130  *
 131  *  The ACML _API_ supports 1D,2D, and 3D transforms, including real-to-complex.
 132  *  Unfortunately the actual library implementation does not support 2D&3D real
 133  *  transforms as of version 4.2  In addition, ACML outputs their hermitian data
 134  *  in such a form that it can't be fed straight back into their complex API's.
 135  *
 136  *  To work around this we roll our own 2D and 3D real-to-complex transforms,
 137  *  using separate X/Y/Z handles defined to perform (ny*nz), (nx*nz), and
 138  *  (nx*ny) transforms at once when necessary. To perform strided multiple
 139  *  transforms out-of-place (i.e., without padding in the last dimension)
 140  *  on the fly we also need to separate the forward and backward
 141  *  handles for real-to-complex/complex-to-real data permutation.
 142  *
 143  *  So, the handles are enumerated as follows:
 144  *
 145  *  1D FFT (real too):    Index 0 is the handle for the entire FFT
 146  *  2D complex FFT:       Index 0 is the handle for the entire FFT
 147  *  3D complex FFT:       Index 0 is the handle for the entire FFT
 148  *  2D, real FFT:                       0=FFTx, 1=FFTrc, 2=FFTcr handle
 149  *  3D, real FFT:                       0=FFTx, 1=FFTy, 2=FFTrc, 3=FFTcr handle
 150  *
 151  *  AMD people reading this: Learn from FFTW what a good interface looks like :-)
 152  *
 153  */
 154 struct gmx_fft
 155 {
 156     /* Work arrays;
 157      * c2c = index0
 158          * r2c = index0, c2r = index1
 159      */
 160     void*               comm[4];
 161     void*               realScratch;
 162
 163     int                ndim;              /**< Number of dimensions in FFT  */
 164     int                nx;                /**< Length of X transform        */
 165     int                ny;                /**< Length of Y transform        */
 166     int                nz;                /**< Length of Z transform        */
 167     int                real_fft;          /**< 1 if real FFT, otherwise 0   */
 168 };
 169
 170 /*      ACML's real FFT support leaves the data in a swizzled 'hermitian' format.  The
 171         problem with this is that you can't just feed this data back into ACML :(  A
 172         pre-processing step is required to transform the hermitian swizzled complex
 173         values into unswizzled complex values
 174         -This function assumes that the complex conjugates are not expanded; the
 175         calling function should not assume that they exist.
 176         -This function pads all real numbers with 0's for imaginary components.
 177         */
 178 void    hermitianUnPacking( float scale, int row, int col, void* src, int srcPitch, void* dst )
 179 {
 180         int mid                 = col/2;
 181         int loopCount   = (col-1)/2;
 182         int hermLength  = mid + 1;
 183
 184         int nFFT;
 185
 186         /* Currently this function is not written to handle aliased input/output pointers */
 187         assert( src != dst );
 188
 189         /* These two pointers could be aliased on top of each other; be careful */
 190         real*   realData                = (real*)src;
 191         t_complex*      hermData        = (t_complex*)dst;
 192
 193         /*      We have to expand the data from real to complex, which will clobber data if you start
 194          *      from the beginning.  For this reason, we start at the last array and work backwards,
 195          *      however, we skip the very first row (nFFT == 0) because we still have data collision
 196          *      on that row.  We copy the first row into a temporary buffer.
 197          */
 198         for( nFFT = row-1; nFFT >= 0; --nFFT )
 199         {
 200                 realData                = (real*)src + (nFFT*srcPitch);
 201                 hermData                = (t_complex*)dst + (nFFT*hermLength);
 202
 203                 /* The first complex number is real valued */
 204                 t_complex       tmp = { realData[0]*scale, 0 };
 205                 hermData[0] = tmp;
 206
 207                 int i;
 208                 for(i = 1; i <= loopCount; ++i)
 209                 {
 210                         t_complex       tmp = { realData[i]*scale, realData[col-i]*scale };
 211                         hermData[i] = tmp;
 212                 }
 213
 214                 /* A little cleanup for even lengths */
 215                 if( loopCount != mid )
 216                 {
 217                         /* The n/2 number is real valued */
 218                         t_complex       tmp = { realData[mid]*scale, 0 };
 219                         hermData[mid] = tmp;
 220                 }
 221         }
 222
 223         return;
 224 }
 225
 226 /*      ACML's real FFT support requires the data in a swizzled 'hermitian' format.
 227         This is a non-standard format, and requires us to manually swizzle the data :(
 228         A pre-processing step is required to transform the complex values into the
 229         swizzled hermitian format
 230         */
 231 void    hermitianPacking( float scale, int row, int col, void* src, void* dst, int dstPitch )
 232 {
 233         int mid                 = col/2;
 234         int loopCount   = (col-1)/2;
 235         int hermLength  = mid + 1;
 236
 237         int nFFT;
 238
 239         real*   realData;
 240         t_complex*      hermData;
 241
 242         /* Currently this function is not written to handle aliased input/output pointers */
 243         assert( src != dst );
 244
 245         for( nFFT = 0; nFFT < row; ++nFFT )
 246         {
 247                 realData                = (real*)dst + (nFFT*dstPitch);
 248                 hermData                = (t_complex*)src + (nFFT*hermLength);
 249
 250                 /* The first complex number is real valued */
 251                 realData[0] = hermData[0].re * scale;
 252
 253                 /* ACML's complex to real function is documented as only handling 'forward'
 254                  * transforms, which mean we have to manually conjugate the data before doing
 255                  * the backwards transform. */
 256                 int i;
 257                 for(i = 1; i <= loopCount; ++i)
 258                 {
 259                         realData[i]     = hermData[i].re * scale;
 260                         realData[col-i] = -hermData[i].im * scale;
 261                 }
 262
 263                 /* A little cleanup for even lengths */
 264                 if( loopCount != mid )
 265                 {
 266                         /* The n/2 number is real valued */
 267                         realData[mid] = hermData[mid].re * scale;
 268                 }
 269
 270         }
 271
 272         return;
 273 }
 274
 275 /*      Gromacs expects results that are returned to be in a packed form; No
 276         complex conjugate's.  This routine takes a 2D array of complex data,
 277         and compresses it to fit into 'real' form by removing all the complex
 278         conjugates and any known imaginary 0's.  Data is expected in row-major
 279         form, with ny describing the length of a row.
 280         */
 281 void    compressConj2D( int nX, int nY, void* src, void* dst )
 282 {
 283         /* These two pointers are aliased on top of each other; be careful */
 284         real*   realData        = (real*)dst;
 285         t_complex*      complexData;
 286
 287         int halfRows    = nX/2 + 1;
 288         int halfCols    = nY/2 + 1;
 289         int rowOdd              = nX & 0x1;
 290         int colOdd              = nY & 0x1;
 291
 292
 293         int nRow;
 294         for( nRow = 0; nRow < nX; ++nRow )
 295         {
 296                 complexData     = (t_complex*)src + (nRow*halfCols);
 297
 298                 int nCol;
 299                 for( nCol = 0; nCol < halfCols; ++nCol )
 300                 {
 301                         if(( nRow == 0) && (nCol == 0 ))
 302                         {
 303                                 /* The complex number is real valued */
 304                                 realData[0] = complexData[nCol].re;
 305                                 realData += 1;
 306                         }
 307                         /* Last column is real if even column length */
 308                         else if( (nRow == 0) && ( !colOdd && ( nCol == (halfCols-1) ) ) )
 309                         {
 310                                 /* The complex number is real valued */
 311                                 realData[0] = complexData[nCol].re;
 312                                 realData += 1;
 313                         }
 314                         /* The middle value of the first column is real if we have an even number of rows and
 315                          * the last column of middle row is real if we have an even number of rows and columns */
 316                         else if( !rowOdd && ( ( nRow == (halfRows-1) ) && ( ( nCol == 0 ) || ( !colOdd && ( nCol == (halfCols-1) ) ) ) ) )
 317                         {
 318                                 /* The complex number is real valued */
 319                                 realData[0] = complexData[nCol].re;
 320                                 realData += 1;
 321                         }
 322                         else if( (nCol == 0) || ( !colOdd && ( nCol == (halfCols-1) ) ) )
 323                         {
 324                                 /* Last half of real columns are conjugates */
 325                                 if( nRow < halfRows )
 326                                 {
 327                                         /* Copy a complex number, which is two reals */
 328                                         realData[0] = complexData[nCol].re;
 329                                         realData[1] = complexData[nCol].im;
 330                                         realData += 2;
 331                                 }
 332                         }
 333                         else
 334                         {
 335                                 /* Copy a complex number, which is two reals */
 336                                 realData[0] = complexData[nCol].re;
 337                                 realData[1] = complexData[nCol].im;
 338                                 realData += 2;
 339                         }
 340                 }
 341         }
 342
 343         return;
 344 }
 345
 346 /*      Gromacs expects results that are returned to be in a packed form; No
 347         complex conjugate's.  This routine takes a 2D array of real data,
 348         and uncompresses it to fit into 'complex' form by creating all the complex
 349         conjugates and any known imaginary 0's.  Data is expected in row-major
 350         form, with ny describing the length of a row.
 351         */
 352 void    unCompressConj2D( int nX, int nY, void* src, void* dst )
 353 {
 354         /* These two pointers are aliased on top of each other; be careful */
 355         real*           realData        = (real*)src;
 356         t_complex*      complexData     = (t_complex*)dst;
 357
 358         int halfRows    = nX/2 + 1;
 359         int halfCols    = nY/2 + 1;
 360         int rowOdd              = nX & 0x1;
 361         int colOdd              = nY & 0x1;
 362
 363         int nRow;
 364         for( nRow = 0; nRow < nX; ++nRow )
 365         {
 366                 int nCol;
 367                 for( nCol = 0; nCol < halfCols; ++nCol )
 368                 {
 369                         int     iIndex  = nRow*halfCols + nCol;
 370
 371                         if(( nRow == 0) && (nCol == 0 ))
 372                         {
 373                                 /* The complex number is real valued */
 374                                 complexData[iIndex].re = realData[0];
 375                                 complexData[iIndex].im =  0;
 376                                 realData += 1;
 377                         }
 378                         /* Last column is real if even column length */
 379                         else if( (nRow == 0) && ( !colOdd && ( nCol == (halfCols-1) ) ) )
 380                         {
 381                                 /* The complex number is real valued */
 382                                 complexData[iIndex].re = realData[0];
 383                                 complexData[iIndex].im =  0;
 384                                 realData += 1;
 385                         }
 386                         /* The middle value of the first column is real if we have an even number of rows and
 387                          * the last column of middle row is real if we have an even number of rows and columns */
 388                         else if( !rowOdd && ( ( nRow == (halfRows-1) ) && ( ( nCol == 0 ) || ( !colOdd && ( nCol == (halfCols-1) ) ) ) ) )
 389                         {
 390                                 /* The complex number is real valued */
 391                                 complexData[iIndex].re = realData[0];
 392                                 complexData[iIndex].im =  0;
 393                                 realData += 1;
 394                         }
 395                         else if( (nCol == 0) || (!colOdd && ( nCol == (halfCols-1) ) ) )
 396                         {
 397                                 /* Last half of real columns are conjugates */
 398                                 if( nRow < halfRows )
 399                                 {
 400                                         /* Copy a complex number, which is two reals */
 401                                         complexData[iIndex].re = realData[0];
 402                                         complexData[iIndex].im = realData[1];
 403                                         realData += 2;
 404                                 }
 405                                 else
 406                                 {
 407                                         int     oIndex  = (nX-nRow)*halfCols + nCol;
 408                                         complexData[iIndex].re = complexData[oIndex].re;
 409                                         complexData[iIndex].im = -complexData[oIndex].im;
 410                                 }
 411                         }
 412                         else
 413                         {
 414                                 /* Copy a complex number, which is two reals */
 415                                 complexData[iIndex].re = realData[0];
 416                                 complexData[iIndex].im = realData[1];
 417                                 realData += 2;
 418                         }
 419                 }
 420         }
 421
 422         return;
 423 }
 424
 425 /* Support routine for the 1D array tranforms.  ACML does not support a MODE
 426  * flag on the real/complex interface.  It assumes a 'forward' transform both
 427  * directions, so that requires a manual conjugation of the imaginary comps. */
 428 void    negateConj( int len, void* src )
 429 {
 430         real* imag      = (real*)src;
 431         int     half    = len/2 + 1;
 432         int i;
 433
 434         for( i = half; i < len; ++i)
 435         {
 436                 imag[i] = -imag[i];
 437         }
 438
 439         return;
 440 }
 441
 442 int
 443 gmx_fft_init_1d(gmx_fft_t *        pfft,
 444                 int                nx,
 445                 enum gmx_fft_flag  flags)
 446 {
 447     gmx_fft_t      fft;
 448     int                    info = 0;
 449     acmlComplex*           comm = NULL;
 450     int                    commSize     = 0;
 451
 452
 453     if(pfft==NULL)
 454     {
 455         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 456         return EINVAL;
 457     }
 458     *pfft = NULL;
 459
 460     if( (fft = malloc(sizeof(struct gmx_fft))) == NULL)
 461     {
 462         return ENOMEM;
 463     }
 464
 465         //      Single precision requires 5*nx+100
 466         //      Double precision requires 3*nx+100
 467         if( sizeof( acmlComplex ) == 16 )
 468                 commSize        = (3*nx+100)*sizeof( acmlComplex );
 469         else
 470                 commSize        = (5*nx+100)*sizeof( acmlComplex );
 471
 472     // Allocate communication work array
 473         if( (comm = (acmlComplex*)malloc( commSize ) ) == NULL )
 474         {
 475                 return ENOMEM;
 476         }
 477
 478     // Initialize communication work array
 479     ACML_FFT1DX( 100, 1.0f, TRUE, nx, NULL, 1, NULL, 1, comm, &info );
 480
 481     if( info != 0 )
 482     {
 483         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 484         gmx_fft_destroy( fft );
 485         return info;
 486     }
 487
 488     fft->ndim     = 1;
 489     fft->nx       = nx;
 490     fft->real_fft = 0;
 491     fft->comm[0]  = comm;
 492
 493     *pfft = fft;
 494     return 0;
 495 }
 496
 497 int
 498 gmx_fft_init_1d_real(gmx_fft_t *        pfft,
 499                      int                nx,
 500                      enum gmx_fft_flag  flags)
 501 {
 502     gmx_fft_t      fft;
 503     int                    info = 0;
 504     real*          commRC       = NULL;
 505     real*          commCR       = NULL;
 506     int                    commSize     = 0;
 507
 508     if(pfft==NULL)
 509     {
 510         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 511         return EINVAL;
 512     }
 513     *pfft = NULL;
 514
 515     if( (fft = malloc(sizeof(struct gmx_fft))) == NULL)
 516     {
 517         return ENOMEM;
 518     }
 519
 520
 521         commSize        = (3*nx+100)*sizeof( float );
 522
 523     // Allocate communication work array, r2c
 524         if( (commRC = (real*)malloc( commSize ) ) == NULL )
 525         {
 526                 return ENOMEM;
 527         }
 528
 529     // Allocate communication work array, c2r
 530         if( (commCR = (real*)malloc( commSize ) ) == NULL )
 531         {
 532                 return ENOMEM;
 533         }
 534
 535     // Initialize communication work array
 536     ACML_RCFFT1D( 100, nx, NULL, (real*)commRC, &info );
 537
 538     if( info != 0 )
 539     {
 540         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 541         gmx_fft_destroy( fft );
 542         return info;
 543     }
 544
 545     // Initialize communication work array
 546     ACML_CRFFT1D( 100, nx, NULL, (real*)commCR, &info );
 547
 548     if( info != 0 )
 549     {
 550         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 551         gmx_fft_destroy( fft );
 552         return info;
 553     }
 554
 555     /* Allocate scratch work array that ACML uses to splat from hermitian complex format to
 556      * full complex format */
 557         if( (fft->realScratch = (acmlComplex*)malloc( nx*sizeof( acmlComplex ) ) ) == NULL )
 558         {
 559                 return ENOMEM;
 560         }
 561
 562     fft->ndim     = 1;
 563     fft->nx       = nx;
 564     fft->real_fft = 1;
 565     fft->comm[0]  = commRC;
 566     fft->comm[1]  = commCR;
 567
 568     *pfft = fft;
 569     return 0;
 570 }
 571
 572 int
 573 gmx_fft_init_2d(gmx_fft_t *        pfft,
 574                 int                nx,
 575                 int                ny,
 576                 enum gmx_fft_flag  flags)
 577 {
 578     gmx_fft_t      fft;
 579     int                    info = 0;
 580     acmlComplex*           comm = NULL;
 581     int                    commSize     = 0;
 582
 583
 584     if(pfft==NULL)
 585     {
 586         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 587         return EINVAL;
 588     }
 589     *pfft = NULL;
 590
 591     if( (fft = malloc(sizeof(struct gmx_fft))) == NULL)
 592     {
 593         return ENOMEM;
 594     }
 595
 596         //      Single precision requires nx*ny+5*(nx+ny)
 597         //      Double precision requires nx*ny+3*(nx+ny)
 598         if( sizeof( acmlComplex ) == 16 )
 599                 commSize        = (nx*ny+3*(nx+ny)+200)*sizeof( acmlComplex );
 600         else
 601                 commSize        = (nx*ny+5*(nx+ny)+200)*sizeof( acmlComplex );
 602
 603     // Allocate communication work array
 604         if( (comm = (acmlComplex*)malloc( commSize ) ) == NULL )
 605         {
 606                 return ENOMEM;
 607         }
 608
 609     // Initialize communication work array
 610     ACML_FFT2DX( 100, 1.0f, TRUE, TRUE, nx, ny, NULL, 1, nx, NULL, 1, nx, (acmlComplex*)comm, &info );
 611
 612     if( info != 0 )
 613     {
 614         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 615         gmx_fft_destroy( fft );
 616         return info;
 617     }
 618
 619     fft->ndim     = 2;
 620     fft->nx       = nx;
 621     fft->ny       = ny;
 622     fft->real_fft = 0;
 623     fft->comm[0]  = comm;
 624
 625     *pfft = fft;
 626     return 0;
 627 }
 628
 629 int
 630 gmx_fft_init_2d_real(gmx_fft_t *        pfft,
 631                      int                nx,
 632                      int                ny,
 633                      enum gmx_fft_flag  flags)
 634 {
 635     gmx_fft_t      fft;
 636     int                    info         = 0;
 637     acmlComplex*   comm         = NULL;
 638     real*                  commRC       = NULL;
 639     real*                  commCR       = NULL;
 640     int                    commSize     = 0;
 641     int                    nyc          = 0;
 642
 643     if(pfft==NULL)
 644     {
 645         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 646         return EINVAL;
 647     }
 648     *pfft = NULL;
 649
 650     if( (fft = malloc(sizeof(struct gmx_fft))) == NULL)
 651     {
 652         return ENOMEM;
 653     }
 654
 655     nyc = (ny/2 + 1);
 656
 657     /* Roll our own 2D real transform using multiple transforms in ACML,
 658      * since the current ACML versions does not support our storage format,
 659      * and all but the most recent don't even have 2D real FFTs.
 660      */
 661
 662         //      Single precision requires 5*nx+100
 663         //      Double precision requires 3*nx+100
 664         if( sizeof( acmlComplex ) == 16 )
 665                 commSize        = (3*nx+100)*sizeof( acmlComplex );
 666         else
 667                 commSize        = (5*nx+100)*sizeof( acmlComplex );
 668
 669     // Allocate communication work array
 670         if( (comm = (acmlComplex*)malloc( commSize ) ) == NULL )
 671         {
 672                 return ENOMEM;
 673         }
 674
 675     // Initialize communication work array
 676     ACML_FFT1DMX( 100, 1.0f, FALSE, nyc, nx, NULL, nyc, 1, NULL, nyc, 1, comm, &info );
 677
 678     if( info != 0 )
 679     {
 680         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 681         gmx_fft_destroy( fft );
 682         return info;
 683     }
 684
 685         commSize        = (3*ny+100)*sizeof( real );
 686     // Allocate communication work array
 687         if( (commRC = (real*)malloc( commSize ) ) == NULL )
 688         {
 689                 return ENOMEM;
 690         }
 691
 692         //      TODO:  Is there no MODE or PLAN for multiple hermetian sequences?
 693     // Initialize communication work array
 694     ACML_RCFFT1D( 100, ny, NULL, commRC, &info );
 695
 696     if( info != 0 )
 697     {
 698         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 699         gmx_fft_destroy( fft );
 700         return info;
 701     }
 702
 703         commSize        = (3*ny+100)*sizeof( real );
 704     // Allocate communication work array
 705         if( (commCR = (real*)malloc( commSize ) ) == NULL )
 706         {
 707                 return ENOMEM;
 708         }
 709
 710         //      TODO:  Is there no MODE or PLAN for multiple hermetian sequences?
 711     // Initialize communication work array
 712     ACML_CRFFT1D( 100, ny, NULL, commCR, &info );
 713
 714     if( info != 0 )
 715     {
 716         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 717         gmx_fft_destroy( fft );
 718         return info;
 719     }
 720
 721     /* Allocate scratch work array that ACML uses to splat from hermitian complex format to
 722      * full complex format */
 723         if( (fft->realScratch = (acmlComplex*)malloc( (nx*ny)*sizeof( acmlComplex ) ) ) == NULL )
 724         {
 725                 return ENOMEM;
 726         }
 727
 728     fft->ndim     = 2;
 729     fft->nx       = nx;
 730     fft->ny       = ny;
 731     fft->real_fft = 1;
 732     fft->comm[0]  = comm;
 733     fft->comm[1]  = commRC;
 734     fft->comm[2]  = commCR;
 735
 736     *pfft = fft;
 737     return 0;
 738 }
 739
 740 int
 741 gmx_fft_init_3d(gmx_fft_t *        pfft,
 742                 int                nx,
 743                 int                ny,
 744                 int                nz,
 745                 enum gmx_fft_flag  flags)
 746 {
 747     gmx_fft_t      fft;
 748     int                    info = 0;
 749     acmlComplex*           comm = NULL;
 750     int                    commSize     = 0;
 751
 752
 753     if(pfft==NULL)
 754     {
 755         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 756         return EINVAL;
 757     }
 758     *pfft = NULL;
 759
 760     if( (fft = malloc(sizeof(struct gmx_fft))) == NULL)
 761     {
 762         return ENOMEM;
 763     }
 764
 765         commSize        = (nx*ny*nz+4*(nx+ny+nz)+300)*sizeof( acmlComplex );
 766
 767     // Allocate communication work array
 768         if( (comm = (acmlComplex*)malloc( commSize ) ) == NULL )
 769         {
 770                 return ENOMEM;
 771         }
 772
 773     ACML_FFT3DX( 100, 1.0f, TRUE, nx, ny, nz, NULL, 1, nx, nx*ny, NULL, 1, nx, nx*ny, comm, commSize, &info );
 774
 775     fft->ndim     = 3;
 776     fft->nx       = nx;
 777     fft->ny       = ny;
 778     fft->nz       = nz;
 779     fft->real_fft = 0;
 780     fft->comm[0]  = comm;
 781
 782     *pfft = fft;
 783     return 0;
 784 }
 785
 786 int
 787 gmx_fft_init_3d_real(gmx_fft_t *        pfft,
 788                      int                nx,
 789                      int                ny,
 790                      int                nz,
 791                      enum gmx_fft_flag  flags)
 792 {
 793     gmx_fft_t      fft;
 794     int                    info         = 0;
 795     acmlComplex*   commX        = NULL;
 796     acmlComplex*   commY        = NULL;
 797     real*                  commRC       = NULL;
 798     real*                  commCR       = NULL;
 799     int                    commSize     = 0;
 800
 801     if(pfft==NULL)
 802     {
 803         gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
 804         return EINVAL;
 805     }
 806     *pfft = NULL;
 807
 808     /* nzc = (nz/2 + 1); */
 809
 810     if( (fft = malloc(sizeof(struct gmx_fft))) == NULL)
 811     {
 812         return ENOMEM;
 813     }
 814
 815     /* Roll our own 3D real transform using multiple transforms in ACML,
 816      * since the current ACML versions does not support 2D
 817      * or 3D real transforms.
 818      */
 819
 820     /* In-place X FFT.
 821      * ny*nz complex-to-complex transforms, length nx
 822      * transform distance: 1
 823      * element strides: ny*nz
 824      */
 825
 826         /*      Single precision requires 5*nx+100
 827                 Double precision requires 3*nx+100
 828         */
 829         if( sizeof( acmlComplex ) == 16 )
 830                 commSize        = (3*nx+100)*sizeof( acmlComplex );
 831         else
 832                 commSize        = (5*nx+100)*sizeof( acmlComplex );
 833
 834     /* Allocate communication work array */
 835         if( (commX = (acmlComplex*)malloc( commSize ) ) == NULL )
 836         {
 837                 return ENOMEM;
 838         }
 839
 840     /* Initialize communication work array */
 841     ACML_FFT1DMX( 100, 1.0f, TRUE, ny*nz, nx, NULL, ny*nz, 1, NULL, ny*nz, 1, commX, &info );
 842
 843     if( info != 0 )
 844     {
 845         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 846         gmx_fft_destroy( fft );
 847         return info;
 848     }
 849
 850     /* In-place Y FFT.
 851      * We cannot do all NX*NZ transforms at once, so define a handle to do
 852      * NZ transforms, and then execute it NX times.
 853      * nz complex-to-complex transforms, length ny
 854      * transform distance: 1
 855      * element strides: nz
 856      */
 857         /*      Single precision requires 5*nx+100
 858                 Double precision requires 3*nx+100
 859         */
 860         if( sizeof( acmlComplex ) == 16 )
 861                 commSize        = (3*ny+100)*sizeof( acmlComplex );
 862         else
 863                 commSize        = (5*ny+100)*sizeof( acmlComplex );
 864
 865     /* Allocate communication work array */
 866         if( (commY = (acmlComplex*)malloc( commSize ) ) == NULL )
 867         {
 868                 return ENOMEM;
 869         }
 870
 871     /* Initialize communication work array */
 872         /* We want to do multiple 1D FFT's in z-y plane, so we have to loop over x
 873          * dimension recalculating z-y plane for each slice.
 874          */
 875     ACML_FFT1DMX( 100, 1.0f, TRUE, nz, ny, NULL, nz, 1, NULL, nz, 1, commY, &info );
 876
 877     if( info != 0 )
 878     {
 879         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 880         gmx_fft_destroy( fft );
 881         return info;
 882     }
 883
 884     /* In-place Z FFT:
 885      * nx*ny real-to-complex transforms, length nz
 886      * transform distance: nzc*2 -> nzc*2
 887      * element strides: 1
 888      */
 889
 890         commSize        = (3*nz+100)*sizeof( real );
 891     /* Allocate communication work array */
 892         if( (commRC = (real*)malloc( commSize ) ) == NULL )
 893         {
 894                 return ENOMEM;
 895         }
 896
 897         /*      TODO:  Is there no MODE or PLAN for multiple hermetian sequences? */
 898         // Initialize communication work array
 899     ACML_RCFFT1D( 100, nz, NULL, commRC, &info );
 900
 901
 902     if( info != 0 )
 903     {
 904         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 905         gmx_fft_destroy( fft );
 906         return info;
 907     }
 908
 909     /* Out-of-place complex-to-real (affects distance) Z FFT:
 910      * nx*ny real-to-complex transforms, length nz
 911      * transform distance: nzc*2 -> nz
 912      * element STRIDES: 1
 913      */
 914         commSize        = (3*nz+100)*sizeof( real );
 915     /* Allocate communication work array */
 916         if( (commCR = (real*)malloc( commSize ) ) == NULL )
 917         {
 918                 return ENOMEM;
 919         }
 920
 921         // Initialize communication work array
 922     ACML_CRFFT1D( 100, nz, NULL, commCR, &info );
 923
 924     if( info != 0 )
 925     {
 926         gmx_fatal(FARGS,"Error initializing ACML FFT; status=%d", info);
 927         gmx_fft_destroy( fft );
 928         return info;
 929     }
 930
 931     /* Allocate scratch work array that ACML uses to splat from hermitian complex format to
 932      * full complex format */
 933         if( (fft->realScratch = (acmlComplex*)malloc( (nx*ny*nz)*sizeof( acmlComplex ) ) ) == NULL )
 934         {
 935                 return ENOMEM;
 936         }
 937
 938     fft->ndim     = 3;
 939     fft->nx       = nx;
 940     fft->ny       = ny;
 941     fft->nz       = nz;
 942     fft->real_fft = 1;
 943     fft->comm[0]  = commX;
 944     fft->comm[1]  = commY;
 945     fft->comm[2]  = commRC;
 946     fft->comm[3]  = commCR;
 947
 948     *pfft = fft;
 949     return 0;
 950 }
 951
 952 int
 953 gmx_fft_1d(gmx_fft_t                  fft,
 954            enum gmx_fft_direction     dir,
 955            void *                     in_data,
 956            void *                     out_data)
 957 {
 958     int inpl = (in_data == out_data);
 959     int info = 0;
 960         int     mode = (dir == GMX_FFT_FORWARD) ? -1 : 1;
 961
 962     if( (fft->real_fft == 1) || (fft->ndim != 1) ||
 963         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
 964     {
 965         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 966         return EINVAL;
 967     }
 968
 969     ACML_FFT1DX( mode, 1.0f, inpl, fft->nx, in_data, 1, out_data, 1, fft->comm[0], &info );
 970
 971     if( info != 0 )
 972     {
 973         gmx_fatal(FARGS,"Error executing AMD ACML FFT.");
 974         info = -1;
 975     }
 976
 977     return info;
 978 }
 979
 980 int
 981 gmx_fft_1d_real(gmx_fft_t                  fft,
 982                 enum gmx_fft_direction     dir,
 983                 void *                     inPtr,
 984                 void *                     outPtr)
 985 {
 986     int info = 0;
 987         int nx = fft->nx;
 988
 989     if( (fft->real_fft != 1) || (fft->ndim != 1) ||
 990         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
 991     {
 992         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
 993         return EINVAL;
 994     }
 995
 996     /* I apply a correction scale to the automatic scaling that was done in the real-complex step */
 997     float       recipCorrection = sqrtf( nx );
 998
 999         /*      ACML needs to do complex2complex intermediate transforms, which will not fit in the amount
1000          *      of memory allocated by the gromacs program, which assumes real.  The real2complex transforms
1001          *      are also only in-place, so we manually do a memcpy() first */
1002     if(dir==GMX_FFT_REAL_TO_COMPLEX)
1003     {
1004         memcpy( fft->realScratch, inPtr, nx*sizeof( real ) );
1005                 ACML_RCFFT1D( 1, nx, fft->realScratch, fft->comm[0], &info );
1006
1007                 hermitianUnPacking( recipCorrection, 1, nx, fft->realScratch, 0, outPtr );
1008     }
1009     else
1010     {
1011         memcpy( fft->realScratch, inPtr, nx*sizeof( t_complex ) );
1012                 hermitianPacking( recipCorrection, 1, nx, fft->realScratch, outPtr, 0 );
1013
1014                 ACML_CRFFT1D( 1, nx, outPtr, fft->comm[1], &info );
1015     }
1016
1017         if( info != 0 )
1018         {
1019         gmx_fatal(FARGS,"Error executing AMD ACML FFT.");
1020         info = -1;
1021         }
1022
1023     return info;
1024 }
1025
1026 int
1027 gmx_fft_2d(gmx_fft_t                  fft,
1028            enum gmx_fft_direction     dir,
1029            void *                     in_data,
1030            void *                     out_data)
1031 {
1032     int inpl = (in_data == out_data);
1033     int info = 0;
1034         int     mode = (dir == GMX_FFT_FORWARD) ? -1 : 1;
1035
1036     if( (fft->real_fft == 1) || (fft->ndim != 2) ||
1037         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
1038     {
1039         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1040         return EINVAL;
1041         }
1042
1043     ACML_FFT2DX( mode, 1.0f, TRUE, inpl, fft->nx, fft->ny,
1044                                 in_data, 1, fft->nx, out_data, 1, fft->nx, fft->comm[0], &info );
1045
1046     if( info != 0 )
1047     {
1048         gmx_fatal(FARGS,"Error executing AMD ACML FFT.");
1049         info = -1;
1050     }
1051
1052     return info;
1053 }
1054
1055
1056 int
1057 gmx_fft_2d_real(gmx_fft_t                  fft,
1058                 enum gmx_fft_direction     dir,
1059                 void *                     inPtr,
1060                 void *                     outPtr )
1061 {
1062     int inpl = (inPtr == outPtr);
1063     int info = 0;
1064     int i       = 0;
1065
1066         int nx  = fft->nx;
1067         int ny  = fft->ny;
1068     int nyc     = (fft->ny/2 + 1);
1069     int nyLen   = 0;
1070
1071     /* Depending on whether we are calculating the FFT in place or not, the rows of the
1072      * 2D FFT will be packed or not.
1073      */
1074     if( inpl )
1075     {
1076         /* Not packed; 1 or 2 extra reals padding at end of each row */
1077         nyLen   = nyc*2;
1078     }
1079     else
1080     {
1081         nyLen   = ny;
1082     }
1083
1084     if( (fft->real_fft != 1) || (fft->ndim != 2) ||
1085         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
1086     {
1087         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1088         return EINVAL;
1089     }
1090
1091     /* I apply a correction scale to the automatic scaling that was done in the real-complex step */
1092     float       recipCorrection = sqrtf( ny );
1093
1094     if(dir==GMX_FFT_REAL_TO_COMPLEX)
1095     {
1096         /*      ACML needs to do complex2complex intermediate transforms, which will not fit in the amount
1097          *      of memory allocated by the gromacs program, which assumes real.  The real2complex transforms
1098          *      are also only in-place, so we manually do a memcpy() first */
1099         if( !inpl )
1100         {
1101                 memcpy( outPtr, inPtr, nx*ny*sizeof( real ) );
1102         }
1103
1104         /* real-to-complex in Y dimension, in-place
1105          * SCFFTM is not valid to call here.  SCFFTM does not take any stride information, and assumes that
1106          * the rows are tightly packed.  GROMACS pads rows with either 1 or 2 extra reals depending
1107          * on even or odd lengths.
1108          */
1109         for( i = 0; i < nx; ++i )
1110         {
1111             if ( info == 0 )
1112                 ACML_RCFFT1D( 1, ny, (real*)outPtr+i*nyLen, fft->comm[1], &info );
1113         }
1114
1115                 hermitianUnPacking( 1.0f, nx, ny, outPtr, nyLen, fft->realScratch );
1116
1117         /* complex-to-complex in X dimension */
1118         if ( info == 0 )
1119                         ACML_FFT1DMX( -1, recipCorrection, FALSE, nyc, nx, fft->realScratch, nyc, 1,
1120                                         outPtr, nyc, 1, fft->comm[0], &info );
1121     }
1122     else
1123     {
1124         /* complex-to-complex in X dimension, in-place */
1125                 ACML_FFT1DMX( 1, recipCorrection, FALSE, nyc, nx, inPtr, nyc, 1,
1126                                         fft->realScratch, nyc, 1, fft->comm[0], &info );
1127
1128                 hermitianPacking( 1.0f, nx, ny, fft->realScratch, outPtr, nyLen );
1129
1130         /* complex-to-real in Y dimension
1131         * CSFFTM is not valid to call here.  SCFFTM does not take any stride information, and assumes that
1132         * the rows are tightly packed.  GROMACS pads rows with either 1 or 2 extra reals depending
1133         * on even or odd lengths.
1134         */
1135         if ( info == 0 )
1136         {
1137                 for( i = 0; i < nx; ++i )
1138                 {
1139                 if ( info == 0 )
1140                         ACML_CRFFT1D( 1, ny, (real*)outPtr+i*nyLen, fft->comm[2], &info );
1141                 }
1142         }
1143     }
1144
1145     if( info != 0 )
1146     {
1147         gmx_fatal(FARGS,"Error executing AMD ACML FFT.");
1148         info = -1;
1149     }
1150
1151     return info;
1152 }
1153
1154
1155 int
1156 gmx_fft_3d(gmx_fft_t                  fft,
1157            enum gmx_fft_direction     dir,
1158            void *                     in_data,
1159            void *                     out_data)
1160 {
1161         int     mode = (dir == GMX_FFT_FORWARD) ? -1 : 1;
1162         int inpl = ( in_data == out_data );
1163
1164     int commSize        = 0;
1165         int nx = fft->nx;
1166         int ny = fft->ny;
1167         int nz = fft->nz;
1168     int info = 0;
1169
1170     if( (fft->real_fft == 1) || (fft->ndim != 3) ||
1171         ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
1172     {
1173         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1174         return EINVAL;
1175     }
1176
1177         commSize        = (nx*ny*nz+4*(nx+ny+nz)+300)*sizeof( acmlComplex );
1178
1179     ACML_FFT3DX( mode, 1.0f, inpl, nx, ny, nz, in_data, 1, nx, nx*ny,
1180                                         out_data, 1, nx, nx*ny, fft->comm[0], commSize, &info );
1181
1182     if( info != 0 )
1183     {
1184         gmx_fatal(FARGS,"Error executing AMD ACML FFT.");
1185         info = -1;
1186     }
1187
1188     return info;
1189 }
1190
1191 int
1192 gmx_fft_3d_real(gmx_fft_t                  fft,
1193                 enum gmx_fft_direction     dir,
1194                 void *                     inPtr,
1195                 void *                     outPtr)
1196 {
1197     int inpl = (inPtr == outPtr);
1198     int info = 0;
1199     int i;
1200     int nx,ny,nz, nzLen = 0;
1201
1202     nx  = fft->nx;
1203     ny  = fft->ny;
1204     nz  = fft->nz;
1205     int nzc     = (nz/2 + 1);
1206
1207     /* Depending on whether we are calculating the FFT in place or not, the rows of the
1208      * 3D FFT will be packed or not.
1209      */
1210     if( inpl )
1211     {
1212         /* Not packed; 1 or 2 extra reals padding at end of each row */
1213         nzLen   = nzc*2;
1214     }
1215     else
1216     {
1217         nzLen   = nz;
1218     }
1219
1220     if( (fft->real_fft != 1) || (fft->ndim != 3) ||
1221         ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
1222     {
1223         gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1224         return EINVAL;
1225     }
1226
1227     /* I apply a correction scale to the automatic scaling that was done in the real-complex step */
1228     float       recipCorrection = sqrtf( nz );
1229
1230     if(dir==GMX_FFT_REAL_TO_COMPLEX)
1231     {
1232         /*      ACML needs to do complex2complex intermediate transforms, which will not fit in the amount
1233          *      of memory allocated by the gromacs program, which assumes real.  The real2complex transforms
1234          *      are also only in-place, so we manually do a memcpy() first */
1235         if( !inpl )
1236         {
1237                 memcpy( outPtr, inPtr, nx*ny*nz*sizeof( real ) );
1238         }
1239
1240                 /* real-to-complex in Z dimension, in-place
1241          * SCFFTM is not valid to call here.  SCFFTM does not take any stride information, and assumes that
1242          * the rows are tightly packed.  GROMACS pads rows with either 1 or 2 extra reals depending
1243          * on even or odd lengths.
1244          */
1245         for( i = 0; i < nx*ny; ++i )
1246         {
1247             if ( info == 0 )
1248                 ACML_RCFFT1D( 1, nz, (real*)outPtr+i*nzLen, fft->comm[2], &info );
1249         }
1250
1251                 hermitianUnPacking( 1.0f, nx*ny, nz, outPtr, nzLen, fft->realScratch );
1252
1253                 /* complex-to-complex in Y dimension, in-place */
1254                 for(i=0;i<nx;i++)
1255                 {
1256                         if ( info == 0 )
1257                                 ACML_FFT1DMX( -1, 1.0f, TRUE, nzc, ny, (acmlComplex*)fft->realScratch+i*ny*nzc, nzc, 1,
1258                                                                 (acmlComplex*)fft->realScratch+i*ny*nzc, nzc, 1, fft->comm[1], &info );
1259                 }
1260
1261                 /* complex-to-complex in X dimension, in-place */
1262                 if ( info == 0 )
1263                         ACML_FFT1DMX( -1, recipCorrection, FALSE, ny*nzc, nx, fft->realScratch, ny*nzc, 1, outPtr, ny*nzc, 1, fft->comm[0], &info );
1264     }
1265     else
1266     {
1267                 /* complex-to-complex in X dimension, from inPtr to work */
1268                 ACML_FFT1DMX( 1, recipCorrection, FALSE, ny*nzc, nx, inPtr, ny*nzc, 1, fft->realScratch, ny*nzc, 1, fft->comm[0], &info );
1269
1270                 /* complex-to-complex in Y dimension, in-place */
1271                 for( i=0; i < nx; i++ )
1272                 {
1273                         if ( info == 0 )
1274                                 ACML_FFT1DMX( 1, 1.0f, TRUE, nzc, ny, (acmlComplex*)fft->realScratch+i*ny*nzc, nzc, 1,
1275                                                         (acmlComplex*)fft->realScratch+i*ny*nzc, nzc, 1, fft->comm[1], &info );
1276                 }
1277
1278                 hermitianPacking( 1.0f, nx*ny, nz, fft->realScratch, outPtr, nzLen );
1279
1280                 /* complex-to-real in Z dimension, in-place
1281          * CSFFTM is not valid to call here.  CSFFTM does not take any stride information, and assumes that
1282          * the rows are tightly packed.  GROMACS pads rows with either 1 or 2 extra reals depending
1283          * on even or odd lengths.
1284          */
1285         for( i = 0; i < nx*ny; ++i )
1286         {
1287             if ( info == 0 )
1288                 ACML_CRFFT1D( 1, nz, (real*)outPtr+i*nzLen, fft->comm[3], &info );
1289         }
1290
1291     }
1292
1293     if( info != 0 )
1294     {
1295         gmx_fatal(FARGS,"Error executing AMD ACML FFT.");
1296         info = -1;
1297     }
1298
1299     return info;
1300 }
1301
1302 void
1303 gmx_fft_destroy(gmx_fft_t    fft)
1304 {
1305     int d;
1306
1307     if(fft != NULL)
1308     {
1309         for(d=0;d<4;d++)
1310         {
1311             if(fft->comm[d] != NULL)
1312             {
1313                                 free(fft->comm[d]);
1314             }
1315         }
1316
1317         if( fft->realScratch != NULL)
1318                 free( fft->realScratch );
1319
1320         free( fft );
1321     }
1322 }
1323
1324 #else
1325 int
1326 gmx_fft_acml_empty;
1327 #endif /* GMX_FFT_ACML */