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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_testfft_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <stdio.h>
  33 #include "typedefs.h"
  34 #include "macros.h"
  35 #include "smalloc.h"
  36 #include "xvgr.h"
  37 #include "complex.h"
  38 #include "fftgrid.h"
  39 #include "mdrun.h"
  40
  41 void testfft(FILE *fp,t_complex ***grid,int nx,int ny,int nz,bool bFirst)
  42 {
  43 #ifdef USE_SGI_FFT
  44 #ifdef DOUBLE
  45   static    zomplex   *coeff;
  46 #else
  47   static    complex   *coeff;
  48 #endif
  49   static    int la1,la2;
  50 #endif
  51   t_complex *cptr;
  52   real      *gptr,*fqqq,fg,fac;
  53   int       ntot,i,j,k,m,n,ndim[4];
  54   int       npppm;
  55
  56   ndim[0] = 0;
  57   ndim[1] = nx;
  58   ndim[2] = ny;
  59   ndim[3] = nz;
  60
  61   ntot    = nx*ny*nz;
  62   cptr    = grid[0][0];
  63   fqqq    = &(grid[0][0][0].re);
  64
  65 #ifdef USE_SGI_FFT
  66   if (bFirst) {
  67     fprintf(fp,"Going to use SGI optimized FFT routines.\n");
  68 #ifdef DOUBLE
  69     coeff  = zfft3di(nx,ny,nz,NULL);
  70 #else
  71     coeff  = cfft3di(nx,ny,nz,NULL);
  72 #endif
  73     bFirst = FALSE;
  74   }
  75   la1 = nx;
  76   la2 = ny;
  77 #ifdef DOUBLE
  78   zfft3d(1,nx,ny,nz,(zomplex *)cptr,la1,la2,coeff);
  79 #else
  80   cfft3d(1,nx,ny,nz,(complex *)cptr,la1,la2,coeff);
  81 #endif
  82 #else
  83   fourn(fqqq-1,ndim,3,1);
  84 #endif
  85
  86 #ifdef USE_SGI_FFT
  87 #ifdef DOUBLE
  88   zfft3d(-1,nx,ny,nz,(zomplex *)cptr,la1,la2,coeff);
  89 #else
  90   cfft3d(-1,nx,ny,nz,(complex *)cptr,la1,la2,coeff);
  91 #endif
  92 #else
  93   fourn(fqqq-1,ndim,3,-1);
  94 #endif
  95 }
  96
  97 void testrft(FILE *fp,real ***grid,int nx,int ny,int nz,bool bFirst)
  98 {
  99 #ifdef USE_SGI_FFT
 100 #ifdef DOUBLE
 101   static    double *coeff;
 102 #else
 103   static    float *coeff;
 104 #endif
 105   static    int la1,la2;
 106 #endif
 107   real      *cptr;
 108   real      *gptr,*fqqq,fg,fac;
 109   int       ntot,i,j,k,m,n,ndim[4];
 110   int       npppm,job;
 111
 112   ndim[0] = 0;
 113   ndim[1] = nx;
 114   ndim[2] = ny;
 115   ndim[3] = nz;
 116
 117   ntot    = nx*ny*nz;
 118   cptr    = grid[0][0];
 119   fqqq    = &(grid[0][0][0]);
 120
 121 #ifdef USE_SGI_FFT
 122   if (bFirst) {
 123     fprintf(fp,"Going to use SGI optimized FFT routines.\n");
 124 #ifdef DOUBLE
 125     coeff  = dfft3di(nx,ny,nz,NULL);
 126 #else
 127     coeff  = sfft3di(nx,ny,nz,NULL);
 128 #endif
 129     bFirst = FALSE;
 130   }
 131   job = 1;
 132   la1 = nx+2;
 133   la2 = ny;
 134 #ifdef DOUBLE
 135   dzfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 136 #else
 137   scfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 138 #endif
 139 #else
 140   fourn(fqqq-1,ndim,3,1);
 141 #endif
 142
 143   job = -1;
 144
 145 #ifdef USE_SGI_FFT
 146 #ifdef DOUBLE
 147   zdfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 148 #else
 149   csfft3d(job,nx,ny,nz,cptr,la1,la2,coeff);
 150 #endif
 151 #else
 152   fourn(fqqq-1,ndim,3,-1);
 153 #endif
 154 }
 155
 156 int main(int argc,char *argv[])
 157 {
 158   FILE      *fp;
 159   int       nnn[] = { 8, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40,
 160                       45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 100 };
 161 #define NNN asize(nnn)
 162   int       *niter;
 163   int       i,j,n,nit,ntot,n3;
 164   double    t,nflop;
 165   double    *rt,*ct;
 166   t_complex ***g;
 167   real      ***h;
 168
 169   snew(rt,NNN);
 170   snew(ct,NNN);
 171   snew(niter,NNN);
 172
 173   for(i=0; (i<NNN); i++) {
 174     n = nnn[i];
 175     fprintf(stderr,"\rReal %d     ",n);
 176     if (n < 16)
 177       niter[i] = 100;
 178     else if (n < 26)
 179       niter[i] = 50;
 180     else if (n < 51)
 181       niter[i] = 10;
 182     else
 183       niter[i] = 5;
 184     nit = niter[i];
 185
 186     h   = mk_rgrid(n+2,n,n);
 187     start_time();
 188     for(j=0; (j<nit); j++) {
 189       testrft(stdout,h,n,n,n,(j==0));
 190     }
 191     update_time();
 192     rt[i] = cpu_time();
 193     free_rgrid(h,n,n);
 194
 195     fprintf(stderr,"\rComplex %d     ",n);
 196     g   = mk_cgrid(n,n,n);
 197     start_time();
 198     for(j=0; (j<nit); j++) {
 199       testfft(stdout,g,n,n,n,(j==0));
 200     }
 201     update_time();
 202     ct[i] = cpu_time();
 203     free_cgrid(g,n,n);
 204   }
 205   fprintf(stderr,"\n");
 206   fp=xvgropen("timing.xvg","FFT timings per grid point","n","t (s)");
 207   for(i=0; (i<NNN); i++) {
 208     n3 = 2*niter[i]*nnn[i]*nnn[i]*nnn[i];
 209     fprintf(fp,"%10d  %10g  %10g\n",nnn[i],rt[i]/n3,ct[i]/n3);
 210   }
 211   fclose(fp);
 212
 213   return 0;
 214 }