00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include "mrilib.h"
00008 #include <stdlib.h>
00009 
00010 /**** test FFT speeds ****/
00011 
00012 int main( int argc , char * argv[] )
00013 {
00014    int num , len , ii , type_ux=0 , kk , iv , nvec , nvec_in ;
00015    complex * cx ;
00016    float fac ;
00017    double tcpu , tclock ;
00018    int narg=1 , quiet=0 ;
00019    float kbytes ;
00020 
00021    if( argc < 4 ){printf("Usage: fftest [-q] len num nvec\n");exit(0);}
00022 
00023    if( strcmp(argv[narg],"-q") == 0 ){ quiet++ ; narg++ ; }
00024    if( strcmp(argv[narg],"-q") == 0 ){ quiet++ ; narg++ ; }
00025    if( strcmp(argv[narg],"-q") == 0 ){ quiet++ ; narg++ ; }
00026    if( strcmp(argv[narg],"-q") == 0 ){ quiet++ ; narg++ ; }
00027 
00028    (void) my_getenv("TMPDIR") ;
00029 
00030    len = strtol( argv[narg++] , NULL , 10 ) ;
00031    num = strtol( argv[narg++] , NULL , 10 ) ;
00032    nvec_in = strtol( argv[narg++] , NULL , 10 ) ;
00033 
00034    if( len > 0 && len != csfft_nextup(len) ){
00035       fprintf(stderr,"Can't do FFT of length %d; try %d\n",len,csfft_nextup(len));
00036       exit(1) ;
00037    } else if ( len < 0 ){
00038       len = -len ;
00039    }
00040 
00041 #ifdef USE_FFTW
00042    { int use_fftw = (nvec_in < 1) ;
00043      if( use_fftw ) nvec_in = 1 ;
00044      else           csfft_use_fftw(0) ;
00045    }
00046 #else
00047    if( nvec_in < 1 ){fprintf(stderr,"Illegal nvec value!\n"); exit(1) ;}
00048 #endif
00049 
00050    cx = (complex *) malloc( sizeof(complex) * len * nvec_in) ;
00051 
00052    kbytes = sizeof(complex) * len * nvec_in / 1024.0 ;
00053 
00054    for( ii=0 ; ii < len*nvec_in ; ii++ ){ cx[ii].r = ii ; cx[ii].i = -0.3*ii ; }
00055 
00056    fac = sqrt(1.0/len) ;
00057    if( ! quiet ) printf("starting %d FFTs of length %d -- %d at a time\n",
00058                         num,len,nvec_in ) ;
00059    tcpu = COX_cpu_time() ;
00060    tclock = COX_clock_time() ;
00061    if( nvec_in == 1 ){
00062       if( ! quiet ) fprintf(stderr,"Using csfft\n") ;
00063       for( kk=0 ; kk < num ; kk++ ){
00064          csfft_cox( -1 , len , cx ) ;
00065          for( ii=0 ; ii < len ; ii++ ){ cx[ii].r *= fac ; cx[ii].i *= fac ; }
00066       }
00067    } else {
00068       if( ! quiet ) fprintf(stderr,"Using csfft_many\n") ;
00069       for( kk=0 ; kk < num ; kk+=nvec_in ){
00070          nvec = num - kk ; if( nvec > nvec_in ) nvec = nvec_in ;
00071          csfft_many( -1 , len , nvec , cx ) ;
00072          for( ii=0 ; ii < len*nvec ; ii++ ){ cx[ii].r *= fac ; cx[ii].i *= fac ; }
00073       }
00074    }
00075    tcpu = COX_cpu_time() - tcpu ;
00076    tclock = COX_clock_time() - tclock ;
00077 
00078    { double flops , m ;
00079      m = log( (double)len ) / log( (double)2.0 ) ;
00080      flops = (5.0 * m - 1.0 ) * len * num / tcpu / 1.e6;
00081 
00082      if(quiet == 1)
00083         printf("%8f %8f\n",kbytes,flops) ;
00084      else if(quiet == 2)
00085         printf("%8f\n",flops) ;
00086      else if(quiet == 3)
00087         printf("%8d %8f\n",len,flops) ;
00088      else
00089         printf(" CPU = %f  ELAPSED = %f  Kbytes = %f  Megaflops = %f\n",
00090                tcpu,tclock,kbytes,flops);
00091      fflush(stdout) ;
00092    }
00093    exit(0) ;
00094 }