00001 #ifdef SPARKY
00002 #undef _POSIX_SOURCE
00003 #endif
00004 #include <sys/types.h>
00005 #include <sys/time.h>
00006 #include <sys/errno.h>
00007 #include <sys/times.h>
00008 #include <limits.h>
00009 
00010 #include "mrilib.h"
00011 
00012 #define FFT_LEN  128
00013 #define FFT_POW  7
00014 #define FFT_NUM  100000
00015 #define MAX_NVEC 256
00016 
00017 double COX_cpu_time(void) ;
00018 
00019 int main( int argc , char * argv[] )
00020 {
00021    int nvec , kk , ii , nvec_now , itop , iarg , fft_len,fft_pow,fft_ops,fft_num,max_nvec ;
00022    complex * cx ;
00023    float fac ;
00024    double cptime , mflops ;
00025 
00026    if( argc > 1 && strncmp(argv[1],"-help",5) == 0 ){
00027       printf("Usage: csfft_cache [-len #] [-num #] [-maxvec #]\n"
00028              "* Does lots of FFTs and determines how many at once are efficient.\n"
00029              "Options:\n"
00030              "  -len #     sets FFT length to '#' (power of 2); default = 128\n"
00031              "  -num #     number of FFTs to do for timing;  default = 100000\n"
00032              "  -maxvec #  maximum number of FFTs to do at once; default = 256\n"
00033             ) ;
00034       exit(0) ;
00035    }
00036 
00037    fft_len  = FFT_LEN ;
00038    fft_pow  = FFT_POW ;
00039    fft_num  = FFT_NUM ;
00040    max_nvec = MAX_NVEC ;
00041    iarg = 1 ;
00042    while( iarg < argc && argv[iarg][0] == '-' ){
00043       if( strncmp(argv[iarg],"-len",4) == 0 ){
00044          fft_len = strtol( argv[++iarg] , NULL , 10 ) ;
00045          if( fft_len < 4 ){fprintf(stderr,"illegal -len\a\n");exit(1);}
00046          for( ii=2,kk=4 ; kk < fft_len ; kk *= 2,ii++ ) ; /* nada */
00047          if( kk != fft_len ){fprintf(stderr,"illegal -len\a\n");exit(1);}
00048          fft_pow = ii ;
00049          iarg++ ; continue ;
00050       }
00051 
00052       if( strncmp(argv[iarg],"-num",4) == 0 ){
00053          fft_num = strtol( argv[++iarg] , NULL , 10 ) ;
00054          if( fft_num < 1 ){fprintf(stderr,"illegal -num\a\n");exit(1);}
00055          iarg++ ; continue ;
00056       }
00057 
00058       if( strncmp(argv[iarg],"-maxvec",5) == 0 ){
00059          max_nvec = strtol( argv[++iarg] , NULL , 10 ) ;
00060          if( max_nvec < 2 ){fprintf(stderr,"illegal -maxvec\a\n");exit(1);}
00061          iarg++ ; continue ;
00062       }
00063    }
00064    fft_ops = (5*fft_pow+2)*fft_len ;
00065    fac = sqrt(1.0/fft_len) ;
00066   
00067    cx = (complex *) malloc( sizeof(complex) * fft_len * max_nvec ) ;
00068    for( ii=0 ; ii < fft_len ; ii++ ){ cx[ii].r = ii ; cx[ii].i = -0.3*ii ; }
00069 
00070    csfft_cox( -1 , fft_len , cx ) ;  /* to initialize trig constants */
00071 
00072    for( nvec=1 ; nvec <= max_nvec ; nvec = (nvec<4) ? (nvec+1) : (1.414*nvec+0.5) ){
00073       for( ii=0 ; ii < fft_len*nvec ; ii++ ){ cx[ii].r = ii ; cx[ii].i = -0.3*ii ; }
00074 
00075       cptime = COX_cpu_time() ;
00076 #if 0
00077       if( nvec == 1 ){
00078 #else
00079       if( nvec == 0 ){
00080 #endif
00081          for( kk=0 ; kk < fft_num ; kk++ ){
00082             csfft_cox( -1 , fft_len , cx ) ;
00083             for( ii=0 ; ii < fft_len ; ii++ ){ cx[ii].r *= fac ; cx[ii].i *= fac ; }
00084          }
00085       } else {
00086          int nvec_now ;
00087          for( kk=0 ; kk < fft_num ; kk+=nvec ){
00088             nvec_now = fft_num - kk ; if( nvec_now > nvec ) nvec_now = nvec ;
00089             csfft_many( -1 , fft_len , nvec_now , cx ) ;
00090             itop = fft_len*nvec_now ;
00091             for( ii=0 ; ii < itop ; ii++ ){ cx[ii].r *= fac ; cx[ii].i *= fac ; }
00092          }
00093       }
00094       cptime = COX_cpu_time() - cptime ;
00095       mflops = fft_num * fft_ops / ( 1.e6 * cptime ) ;
00096       printf("FFT cache size %6d bytes gives %10.2f Mflops"
00097              " (%d FFTs, size %d, %d at a time, in %g sec)\n",
00098              sizeof(complex)*fft_len*nvec , mflops , fft_num,fft_len,nvec,cptime) ;
00099    }
00100 
00101    exit(0) ;
00102 }
00103 
00104 double COX_cpu_time(void)  /* in seconds */
00105 {
00106    struct tms ttt ;
00107    (void) times( &ttt ) ;
00108    return (  (double) (ttt.tms_utime) / (double) CLK_TCK ) ;
00109 }