home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #3 / NN_1993_3.iso / spool / comp / benchmar / 1989 < prev    next >
Encoding:
Text File  |  1993-01-22  |  11.6 KB  |  503 lines
  1. Newsgroups: comp.benchmarks
  2. Path: sparky!uunet!spool.mu.edu!agate!dog.ee.lbl.gov!news!marlin!aburto
  3. From: aburto@nosc.mil (Alfred A. Aburto)
  4. Subject: Test FFT DP C (tfftdp.c) Source Code
  5. Message-ID: <1993Jan22.165132.26692@nosc.mil>
  6. Organization: Naval Ocean Systems Center, San Diego
  7. Distribution: comp.benchmarks
  8. Date: Fri, 22 Jan 1993 16:51:32 GMT
  9. Expires: Mon, 15 Feb 1993 08:00:00 GMT
  10. Lines: 491
  11.  
  12. /*
  13.  C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c)
  14.  by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
  15. */
  16.  
  17. #include <stdio.h>
  18. #include <math.h>
  19.  
  20. /***************************************************************/
  21. /* Timer options. You MUST uncomment one of the options below  */
  22. /* or compile, for example, with the '-DUNIX' option.          */
  23. /*                                                             */
  24. /* Example compilation:                                        */
  25. /* UNIX systems:                                               */
  26. /* cc -DUNIX -O tfftdp.c -o fft -lm                            */
  27. /***************************************************************/
  28. /* #define Amiga       */
  29. /* #define UNIX        */
  30. /* #define UNIX_Old    */
  31. /* #define VMS         */
  32. /* #define BORLAND_C   */
  33. /* #define MSC         */
  34. /* #define MAC         */
  35. /* #define IPSC        */
  36. /* #define FORTRAN_SEC */
  37. /* #define GTODay      */
  38. /* #define CTimer      */
  39. /* #define UXPM        */
  40.  
  41. #ifndef PI
  42. #define PI  3.14159265358979323846
  43. #endif
  44.  
  45. double xr[262144],xi[262144];
  46.  
  47. void main()
  48. {
  49. int i,j,np,npm,lp,n2,kr,ki;
  50. double a,enp,t,x,y,zr,zi,zm;
  51. double t1,t2,t3,benchtime,VAX_REF,dtime();
  52. int dfft();
  53.  
  54. np = 8;
  55.       printf("\nFFT benchmark - Double Precision - V1.0 -  05 Jan 1993\n\n");
  56.       printf("FFT size   Time(sec)   max error\n");
  57.  
  58.       t3 = dtime();
  59.       for (lp = 1; lp <= 15; lp++ ){ 
  60.       np = np * 2; 
  61.       printf(" %7d ",np);
  62.       enp = np; 
  63.       npm = np / 2  - 1;  
  64.       t = PI / enp;
  65.       xr[0] = (enp - 1.0) * 0.5;  
  66.       xi[0] = 0;
  67.       n2 = np / 2;  
  68.       xr[n2] = -0.5; 
  69.       xi[n2] = 0.0;
  70.       
  71.       for (i = 1; i <= npm; i++) {
  72.           j = np - i;
  73.           xr[i] = -0.5; 
  74.           xr[j] = -0.5;
  75.           y = t * i;  
  76.           y = -0.5*(cos(y)/sin(y));
  77.           xi[i] = y; xi[j] = -y;  }
  78.       
  79.       t1 = dtime();
  80.       dfft(xr,xi,np);
  81.       t2 = dtime() - t1;
  82.       if (t2 < 0.0) t2 = 0.0;
  83.       printf(" %10.4lf ",t2);
  84. /*
  85.       for (i=0; i<=15; i++) printf("%d  %g  %g\n",i,xr[i],xi[i]);
  86. */
  87.       zr = 0.0; zi = 0.0; kr = 0; ki = 0; npm = np-1;
  88.       for (i = 0; i <= npm; i++ ) { 
  89.           a = fabs( xr[i] - i );
  90.           if (zr < a) { zr = a; kr = i; }
  91.           a = fabs(xi[i]);
  92.           if (zi < a) { zi = a; ki = i; }
  93.           }
  94.       zm = zr;
  95.       if ( fabs(zr) < fabs(zi) ) zm = zi;
  96.       printf(" %10.2g %10.4\n",zm,t2);
  97.  
  98. /*
  99.       printf("re %7d  %10.2g  im %7d  %10.2g\n",kr,zr,ki,zi);
  100. */
  101.       }
  102.       benchtime = dtime() - t3;
  103.       VAX_REF = 140.658;
  104.       printf("\nBenchTime (sec) = %10.4lf\n",benchtime);
  105.       printf("VAX_FFTs = %10.3lf\n\n",VAX_REF/benchtime);
  106.  
  107. }
  108.  
  109.  
  110.  
  111. /*
  112.        A Duhamel-Hollman split-radix dif fft
  113.        Ref: Electronics Letters, Jan. 5, 1984
  114.        Complex input and output data in arrays x and y
  115.        Length is n.
  116. */
  117.  
  118. int dfft( x, y, np )
  119. double x[2]; double y[2]; int np ;
  120. {
  121. double *px,*py;
  122. int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4 ;
  123. double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt ;
  124.   px = x - 1; py = y - 1;
  125.   i = 2; m = 1; while (i < np) { i = i+i; m = m+1; };
  126.   n = i; if (n != np) {
  127.   for (i = np+1; i <= n; i++)  { *(px + i) = 0.0; *(py + i) = 0.0; };
  128.   printf("\nuse %d point fft",n); }
  129.  
  130.   n2 = n+n;
  131.   for (k = 1;  k <= m-1; k++ ) {
  132.     n2 = n2 / 2; n4 = n2 / 4; e = 2.0 * PI / n2; a = 0.0;
  133.     for (j = 1; j<= n4 ; j++) {
  134.       a3 = 3.0*a; cc1 = cos(a); ss1 = sin(a);
  135.       cc3 = cos(a3); ss3 = sin(a3); a = j*e; is = j; id = 2*n2;
  136.       while ( is < n ) {
  137.       for (i0 = is; i0 <= n-1; i0 = i0 + id) {
  138.          i1 = i0 + n4; i2 = i1 + n4; i3 = i2 + n4;
  139.          r1 = *(px+i0) - *(px+i2);
  140.          *(px+i0) = *(px+i0) + *(px+i2);
  141.          r2 = *(px+i1) - *(px+i3);
  142.          *(px+i1) = *(px+i1) + *(px+i3);
  143.          s1 = *(py+i0) - *(py+i2);
  144.          *(py+i0) = *(py+i0) + *(py+i2);
  145.          s2 = *(py+i1) - *(py+i3);
  146.          *(py+i1) = *(py+i1) + *(py+i3);
  147.          s3 = r1 - s2; r1 = r1 + s2; s2 = r2 - s1; r2 = r2 + s1;
  148.          *(px+i2) = r1*cc1 - s2*ss1; 
  149.          *(py+i2) = -s2*cc1 - r1*ss1;
  150.          *(px+i3) = s3*cc3 + r2*ss3;
  151.              *(py+i3) = r2*cc3 - s3*ss3;
  152.      }
  153.        is = 2*id - n2 + j; id = 4*id;
  154.     }
  155.     }
  156.   }
  157.  
  158. /*
  159. ---------------------Last stage, length=2 butterfly---------------------
  160. */
  161.   is = 1; id = 4;
  162.   while ( is < n) {
  163.   for (i0 = is; i0 <= n; i0 = i0 + id) {
  164.       i1 = i0 + 1; r1 = *(px+i0);
  165.       *(px+i0) = r1 + *(px+i1);
  166.       *(px+i1) = r1 - *(px+i1);
  167.       r1 = *(py+i0);
  168.       *(py+i0) = r1 + *(py+i1);
  169.       *(py+i1) = r1 - *(py+i1);
  170.     }
  171.   is = 2*id - 1; id = 4 * id; }
  172.  
  173. /*
  174. c--------------------------Bit reverse counter
  175. */
  176.   j = 1; n1 = n - 1;
  177.   for (i = 1; i <= n1; i++) {
  178.     if (i < j) {
  179.       xt = *(px+j);
  180.       *(px+j) = *(px+i); *(px+i) = xt;
  181.       xt = *(py+j); *(py+j) = *(py+i);
  182.       *(py+i) = xt;
  183.     }
  184.     k = n / 2; while (k < j) { j = j - k; k = k / 2; }
  185.     j = j + k;
  186.   }
  187.  
  188. /*
  189.   for (i = 1; i<=16; i++) printf("%d  %g   %gn",i,*(px+i),(py+i));
  190. */
  191.   
  192.   return(n);
  193.   }
  194.  
  195.  
  196. /*****************************************************/
  197. /* Various timer routines.                           */
  198. /* Al Aburto, aburto@marlin.nosc.mil, 20 Dec 1992    */
  199. /*                                                   */
  200. /* t = dtime() outputs the current time in seconds.  */
  201. /* Use CAUTION as some of these routines will mess   */
  202. /* up when timing across the hour mark!!!            */
  203. /*                                                   */
  204. /* For timing I use the 'user' time whenever         */
  205. /* possible. Using 'user+sys' time is a separate     */
  206. /* issue.                                            */
  207. /*                                                   */
  208. /* Example Usage:                                    */
  209. /* [timer options added here]                        */
  210. /* main()                                            */
  211. /* {                                                 */
  212. /* double starttime,benchtime,dtime();               */
  213. /*                                                   */
  214. /* starttime = dtime();                              */ 
  215. /* [routine to time]                                 */
  216. /* benchtime = dtime() - starttime;                  */
  217. /* }                                                 */
  218. /*                                                   */
  219. /* [timer code below added here]                     */
  220. /*****************************************************/
  221.  
  222. /*********************************/
  223. /* Timer code.                   */
  224. /*********************************/
  225. /*******************/
  226. /*  Amiga dtime()  */
  227. /*******************/
  228. #ifdef Amiga
  229. #include <ctype.h>
  230. #define HZ 50
  231.  
  232. double dtime()
  233. {
  234.    double q;
  235.  
  236.    struct   tt {
  237.       long  days;
  238.       long  minutes;
  239.       long  ticks;
  240.    } tt;
  241.  
  242.    DateStamp(&tt);
  243.  
  244.    q = ((double)(tt.ticks + (tt.minutes * 60L * 50L))) / (double)HZ;
  245.  
  246.    return q;
  247. }
  248. #endif
  249.  
  250. /*****************************************************/
  251. /*  UNIX dtime(). This is the preferred UNIX timer.  */
  252. /*  Provided by: Markku Kolkka, mk59200@cc.tut.fi    */
  253. /*  HP-UX Addition by: Bo Thide', bt@irfu.se         */
  254. /*****************************************************/
  255. #ifdef UNIX
  256. #include <sys/time.h>
  257. #include <sys/resource.h>
  258.  
  259. #ifdef hpux
  260. #include <sys/syscall.h>
  261. #define getrusage(a,b) syscall(SYS_getrusage,a,b)
  262. #endif
  263.  
  264. struct rusage rusage;
  265.  
  266. double dtime()
  267. {
  268.    double q;
  269.  
  270.    getrusage(RUSAGE_SELF,&rusage);
  271.  
  272.    q = (double)(rusage.ru_utime.tv_sec);
  273.    q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
  274.    
  275.    return q;
  276. }
  277. #endif
  278.  
  279. /***************************************************/
  280. /*  UNIX_Old dtime(). This is the old UNIX timer.  */
  281. /*  Use only if absolutely necessary as HZ may be  */
  282. /*  ill defined on your system.                    */
  283. /***************************************************/
  284. #ifdef UNIX_Old
  285. #include <sys/types.h>
  286. #include <sys/times.h>
  287. #include <sys/param.h>
  288.  
  289. #ifndef HZ
  290. #define HZ 60
  291. #endif
  292.  
  293. struct tms tms;
  294.  
  295. double dtime()
  296. {
  297.    double q;
  298.  
  299.    times(&tms);
  300.  
  301.    q = (double)(tms.tms_utime) / (double)HZ;
  302.    
  303.    return q;
  304. }
  305. #endif
  306.  
  307. /*********************************************************/
  308. /*  VMS dtime() for VMS systems.                         */
  309. /*  Provided by: RAMO@uvphys.phys.UVic.CA                */
  310. /*  Some people have run into problems with this timer.  */
  311. /*********************************************************/
  312. #ifdef VMS
  313. #include time
  314.  
  315. #ifndef HZ
  316. #define HZ 100
  317. #endif
  318.  
  319. struct tbuffer_t
  320.        {
  321.     int proc_user_time;
  322.     int proc_system_time;
  323.     int child_user_time;
  324.     int child_system_time;
  325.        };
  326. struct tbuffer_t tms;
  327.  
  328. double dtime()
  329. {
  330.    double q;
  331.  
  332.    times(&tms);
  333.  
  334.    q = (double)(tms.proc_user_time) / (double)HZ;
  335.    
  336.    return q;
  337. }
  338. #endif
  339.  
  340. /******************************/
  341. /*  BORLAND C dtime() for DOS */
  342. /******************************/
  343. #ifdef BORLAND_C
  344. #include <ctype.h>
  345. #include <dos.h>
  346. #include <time.h>
  347.  
  348. #define HZ 100
  349. struct time tnow;
  350.  
  351. double dtime()
  352. {
  353.    double q;
  354.  
  355.    gettime(&tnow);
  356.  
  357.    q = 60.0 * (double)(tnow.ti_min);
  358.    q = q + (double)(tnow.ti_sec);
  359.    q = q + (double)(tnow.ti_hund)/(double)HZ;
  360.    
  361.    return q;
  362. }
  363. #endif
  364.  
  365. /**************************************/
  366. /*  Microsoft C (MSC) dtime() for DOS */
  367. /**************************************/
  368. #ifdef MSC
  369. #include <time.h>
  370. #include <ctype.h>
  371.  
  372. #define HZ CLK_TCK
  373. clock_t tnow;
  374.  
  375. double dtime()
  376. {
  377.    double q;
  378.  
  379.    tnow = clock();
  380.  
  381.    q = (double)tnow / (double)HZ;
  382.    
  383.    return q;
  384. }
  385. #endif
  386.  
  387. /*************************************/
  388. /*  Macintosh (MAC) Think C dtime()  */
  389. /*************************************/
  390. #ifdef MAC
  391. #include <time.h>
  392.  
  393. #define HZ 60
  394.  
  395. double dtime()
  396. {
  397.    double q;
  398.  
  399.    q = (double)clock() / (double)HZ;
  400.    
  401.    return q;
  402. }
  403. #endif
  404.  
  405. /************************************************************/
  406. /*  iPSC/860 (IPSC) dtime() for i860.                       */
  407. /*  Provided by: Dan Yergeau, yergeau@gloworm.Stanford.EDU  */
  408. /************************************************************/
  409. #ifdef IPSC
  410. extern double dclock();
  411.  
  412. double dtime()
  413. {
  414.    double q;
  415.  
  416.    q = dclock();
  417.    
  418.    return q;
  419. }
  420. #endif
  421.  
  422. /**************************************************/
  423. /*  FORTRAN dtime() for Cray type systems.        */
  424. /*  This is the preferred timer for Cray systems. */
  425. /**************************************************/
  426. #ifdef FORTRAN_SEC
  427.  
  428. fortran double second();
  429.  
  430. double dtime()
  431. {
  432.    double q;
  433.  
  434.    second(&q);
  435.    
  436.    return q;
  437. }
  438. #endif
  439.  
  440. /***********************************************************/
  441. /*  UNICOS C dtime() for Cray UNICOS systems.  Don't use   */
  442. /*  unless absolutely necessary as returned time includes  */
  443. /*  'user+system' time.  Provided by: R. Mike Dority,      */
  444. /*  dority@craysea.cray.com                                */
  445. /***********************************************************/
  446. #ifdef CTimer
  447. #include <time.h>
  448.  
  449. double dtime()
  450. {
  451.    double    q;
  452.    clock_t   t;
  453.  
  454.        t = clock();
  455.  
  456.        q = (double)t / (double)CLOCKS_PER_SEC;
  457.  
  458.        return q;
  459. }
  460. #endif
  461.  
  462. /********************************************/
  463. /* Another UNIX timer using gettimeofday(). */
  464. /* However, getrusage() is preferred.       */
  465. /********************************************/
  466. #ifdef GTODay
  467. #include <sys/time.h>
  468.  
  469. struct timeval tnow;
  470.  
  471. double dtime()
  472. {
  473.    double q;
  474.  
  475.    gettimeofday(&tnow,NULL);
  476.    q = (double)tnow.tv_sec + (double)tnow.tv_usec * 1.0e-6;
  477.  
  478.    return q;
  479. }
  480. #endif
  481.  
  482. /*****************************************************/
  483. /*  Fujitsu UXP/M timer.                             */
  484. /*  Provided by: Mathew Lim, ANUSF, M.Lim@anu.edu.au */
  485. /*****************************************************/
  486. #ifdef UXPM
  487. #include <sys/types.h>
  488. #include <sys/timesu.h>
  489. struct tmsu rusage;
  490.  
  491. double dtime()
  492. {
  493.    double q;
  494.  
  495.    timesu(&rusage);
  496.  
  497.    q = (double)(rusage.tms_utime) * 1.0e-06;
  498.    
  499.    return q;
  500. }
  501. #endif
  502.  
  503.