src/gmxlib/testfft.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #include <math.h>
  37 #include <stdio.h>
  38 #include "typedefs.h"
  39 #include "macros.h"
  40 #include "smalloc.h"
  41 #include "xvgr.h"
  42 #include "complex.h"
  43 #include "fftgrid.h"
  44 #include "mdrun.h"
  45
  46 void testfft(FILE *fp,t_complex ***grid,int nx,int ny,int nz,bool bFirst)
  47 {
  48 #ifdef USE_SGI_FFT
  49 #ifdef DOUBLE
  50   static    zomplex   *coeff;
  51 #else
  52   static    complex   *coeff;
  53 #endif
  54   static    int la1,la2;
  55 #endif
  56   t_complex *cptr;
  57   real      *gptr,*fqqq,fg,fac;
  58   int       ntot,i,j,k,m,n,ndim[4];
  59   int       npppm;
  60
  61   ndim[0] = 0;
  62   ndim[1] = nx;
  63   ndim[2] = ny;
  64   ndim[3] = nz;
  65
  66   ntot    = nx*ny*nz;
  67   cptr    = grid[0][0];
  68   fqqq    = &(grid[0][0][0].re);
  69
  70 #ifdef USE_SGI_FFT
  71   if (bFirst) {
  72     fprintf(fp,"Going to use SGI optimized FFT routines.\n");
  73 #ifdef DOUBLE
  74     coeff  = zfft3di(nx,ny,nz,NULL);
  75 #else
  76     coeff  = cfft3di(nx,ny,nz,NULL);
  77 #endif
  78     bFirst = FALSE;
  79   }
  80   la1 = nx;
  81   la2 = ny;
  82 #ifdef DOUBLE
  83   zfft3d(1,nx,ny,nz,(zomplex *)cptr,la1,la2,coeff);
  84 #else
  85   cfft3d(1,nx,ny,nz,(complex *)cptr,la1,la2,coeff);
  86 #endif
  87 #else
  88   fourn(fqqq-1,ndim,3,1);
  89 #endif
  90
  91 #ifdef USE_SGI_FFT
  92 #ifdef DOUBLE
  93   zfft3d(-1,nx,ny,nz,(zomplex *)cptr,la1,la2,coeff);
  94 #else
  95   cfft3d(-1,nx,ny,nz,(complex *)cptr,la1,la2,coeff);
  96 #endif
  97 #else
  98   fourn(fqqq-1,ndim,3,-1);
  99 #endif
 100 }
 101
 102 void testrft(FILE *fp,real ***grid,int nx,int ny,int nz,bool bFirst)
 103 {
 104 #ifdef USE_SGI_FFT
 105 #ifdef DOUBLE
 106   static    double *coeff;
 107 #else
 108   static    float *coeff;
 109 #endif
 110   static    int la1,la2;
 111 #endif
 112   real      *cptr;
 113   real      *gptr,*fqqq,fg,fac;
 114   int       ntot,i,j,k,m,n,ndim[4];
 115   int       npppm,job;
 116
 117   ndim[0] = 0;
 118   ndim[1] = nx;
 119   ndim[2] = ny;
 120   ndim[3] = nz;
 121
 122   ntot    = nx*ny*nz;
 123   cptr    = grid[0][0];
 124   fqqq    = &(grid[0][0][0]);
 125
 126 #ifdef USE_SGI_FFT
 127   if (bFirst) {
 128     fprintf(fp,"Going to use SGI optimized FFT routines.\n");
 129 #ifdef DOUBLE
 130     coeff  = dfft3di(nx,ny,nz,NULL);
 131 #else
 132     coeff  = sfft3di(nx,ny,nz,NULL);
 133 #endif
 134     bFirst = FALSE;
 135   }
 136   job = 1;
 137   la1 = nx+2;
 138   la2 = ny;
 139 #ifdef DOUBLE
 140   dzfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 141 #else
 142   scfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 143 #endif
 144 #else
 145   fourn(fqqq-1,ndim,3,1);
 146 #endif
 147
 148   job = -1;
 149
 150 #ifdef USE_SGI_FFT
 151 #ifdef DOUBLE
 152   zdfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 153 #else
 154   csfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 155 #endif
 156 #else
 157   fourn(fqqq-1,ndim,3,-1);
 158 #endif
 159 }
 160
 161 int main(int argc,char *argv[])
 162 {
 163   FILE      *fp;
 164   int       nnn[] = { 8, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40,
 165                       45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 100 };
 166 #define NNN asize(nnn)
 167   int       *niter;
 168   int       i,j,n,nit,ntot,n3;
 169   double    t,nflop;
 170   double    *rt,*ct;
 171   t_complex ***g;
 172   real      ***h;
 173
 174   snew(rt,NNN);
 175   snew(ct,NNN);
 176   snew(niter,NNN);
 177
 178   for(i=0; (i<NNN); i++) {
 179     n = nnn[i];
 180     fprintf(stderr,"\rReal %d     ",n);
 181     if (n < 16)
 182       niter[i] = 100;
 183     else if (n < 26)
 184       niter[i] = 50;
 185     else if (n < 51)
 186       niter[i] = 10;
 187     else
 188       niter[i] = 5;
 189     nit = niter[i];
 190
 191     h   = mk_rgrid(n+2,n,n);
 192     start_time();
 193     for(j=0; (j<nit); j++) {
 194       testrft(stdout,h,n,n,n,(j==0));
 195     }
 196     update_time();
 197     rt[i] = node_time();
 198     free_rgrid(h,n,n);
 199
 200     fprintf(stderr,"\rComplex %d     ",n);
 201     g   = mk_cgrid(n,n,n);
 202     start_time();
 203     for(j=0; (j<nit); j++) {
 204       testfft(stdout,g,n,n,n,(j==0));
 205     }
 206     update_time();
 207     ct[i] = node_time();
 208     free_cgrid(g,n,n);
 209   }
 210   fprintf(stderr,"\n");
 211   fp=xvgropen("timing.xvg","FFT timings per grid point","n","t (s)");
 212   for(i=0; (i<NNN); i++) {
 213     n3 = 2*niter[i]*nnn[i]*nnn[i]*nnn[i];
 214     fprintf(fp,"%10d  %10g  %10g\n",nnn[i],rt[i]/n3,ct[i]/n3);
 215   }
 216   fclose(fp);
 217
 218   return 0;
 219 }