diff options
| -rw-r--r-- | fourier.c | 304 | ||||
| -rw-r--r-- | fpemulation.c | 50 |
2 files changed, 142 insertions, 212 deletions
@@ -4,6 +4,7 @@ #include <string.h> #include <math.h> #include <limits.h> +#include <time.h> #include "nmglobal.h" #include "nbench1.h" @@ -12,16 +13,16 @@ ** FOURIER COEFFICIENTS ** *************************/ -static unsigned long DoFPUTransIteration(double *abase, +static clock_t DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize); static double TrapezoidIntegrate(double x0, double x1, int nsteps, - double omegan, + double omega_n, int select); static double thefunction(double x, - double omegan, + double omega_n, int select); /************** @@ -32,95 +33,80 @@ static double thefunction(double x, ** fourier coefficients of the function (x+1)^x defined ** on the interval 0,2. */ -void DoFourier(void) +void +DoFourier(void) { -FourierStruct *locfourierstruct=&global_fourierstruct; /* Local fourier struct */ -double *abase = NULL; /* Base of A[] coefficients array */ -double *bbase = NULL; /* Base of B[] coefficients array */ -unsigned long accumtime; /* Accumulated time in ticks */ -double iterations; /* # of iterations */ -char *context = "FPU:Transcendental"; /* Error context string pointer */ -int systemerror; /* For error code */ - -/* -** See if we need to do self-adjustment code. -*/ -if(locfourierstruct->adjust==0) -{ - locfourierstruct->arraysize=100L; /* Start at 100 elements */ - while(1) - { - abase = realloc(abase, locfourierstruct->arraysize * sizeof(double)); - if (!abase) { + const char* context = "FPU:Transcendental"; + FourierStruct* locfourierstruct = &global_fourierstruct; + clock_t total_time = 0; + int iterations = 0; + double* abase = NULL; + double* bbase = NULL; + + /* + ** See if we need to do self-adjustment code. + */ + if (locfourierstruct->adjust == FALSE) { + locfourierstruct->adjust = TRUE; + locfourierstruct->arraysize = 100L; /* Start at 100 elements */ + + while (1) { + abase = realloc(abase, locfourierstruct->arraysize * sizeof(double)); + if (!abase) { + fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context); + exit(1); + } + + bbase = realloc(bbase, locfourierstruct->arraysize * sizeof(double)); + if (!bbase) { + fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context); + free(abase); + exit(1); + } + + /* + ** Do an iteration of the tests. If the elapsed time is + ** less than or equal to the permitted minimum, re-allocate + ** larger arrays and try again. + */ + + if (DoFPUTransIteration(abase,bbase, locfourierstruct->arraysize) > global_min_ticks) { + break; + } + + /* + ** Make bigger arrays and try again. + */ + locfourierstruct->arraysize += 50L; + } + } else { + /* + ** Don't need self-adjustment. Just allocate the + ** arrays, and go. + */ + abase = malloc(locfourierstruct->arraysize * sizeof(double)); + if (!abase) { fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context); exit(1); - } + } - bbase = realloc(bbase, locfourierstruct->arraysize * sizeof(double)); + bbase = malloc(locfourierstruct->arraysize * sizeof(double)); if (!bbase) { fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context); free(abase); exit(1); } - /* - ** Do an iteration of the tests. If the elapsed time is - ** less than or equal to the permitted minimum, re-allocate - ** larger arrays and try again. - */ - if(DoFPUTransIteration(abase,bbase, - locfourierstruct->arraysize)>global_min_ticks) - break; /* We're ok...exit */ - - /* - ** Make bigger arrays and try again. - */ - locfourierstruct->arraysize+=50L; - } -} -else -{ /* - ** Don't need self-adjustment. Just allocate the - ** arrays, and go. - */ - abase = malloc(locfourierstruct->arraysize * sizeof(double)); - if (!abase) { - fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context); - exit(1); - } - - bbase = malloc(locfourierstruct->arraysize * sizeof(double)); - if (!bbase) { - fprintf(stderr, "Error in %s, could not allocate memory. Exitting...\n", context); - free(abase); - exit(1); - } -} -/* -** All's well if we get here. Repeatedly perform integration -** tests until the accumulated time is greater than the -** # of seconds requested. -*/ -accumtime=0L; -iterations=(double)0.0; -do { - accumtime+=DoFPUTransIteration(abase,bbase,locfourierstruct->arraysize); - iterations+=(double)locfourierstruct->arraysize*(double)2.0-(double)1.0; -} while(TicksToSecs(accumtime)<locfourierstruct->request_secs); - - -/* -** Clean up, calculate results, and go home. -** Also set adjustment flag to indicate no adjust code needed. -*/ -free(abase); -free(bbase); + } -locfourierstruct->fflops=iterations/(double)TicksToFracSecs(accumtime); + do { + total_time += DoFPUTransIteration(abase,bbase,locfourierstruct->arraysize); + iterations += locfourierstruct->arraysize * 2 - 1; + } while (total_time < locfourierstruct->request_secs * CLOCKS_PER_SEC); -if(locfourierstruct->adjust==0) - locfourierstruct->adjust=1; + free(abase); + free(bbase); -return; + locfourierstruct->fflops = (double)(iterations * CLOCKS_PER_SEC) / (double)total_time; } /************************ @@ -133,75 +119,47 @@ return; ** NOTE: The # of integration steps is fixed at ** 200. */ -static unsigned long DoFPUTransIteration(double *abase, /* A coeffs. */ - double *bbase, /* B coeffs. */ - unsigned long arraysize) /* # of coeffs */ +static clock_t +DoFPUTransIteration(double *abase, double *bbase, unsigned long arraysize) { -double omega; /* Fundamental frequency */ -unsigned long i; /* Index */ -unsigned long elapsed; /* Elapsed time */ - -/* -** Start the stopwatch -*/ -elapsed=StartStopwatch(); - -/* -** Calculate the fourier series. Begin by -** calculating A[0]. -*/ - -*abase=TrapezoidIntegrate((double)0.0, - (double)2.0, - 200, - (double)0.0, /* No omega * n needed */ - 0 )/(double)2.0; - -/* -** Calculate the fundamental frequency. -** ( 2 * pi ) / period...and since the period -** is 2, omega is simply pi. -*/ -omega=(double)3.1415926535897932; - -for(i=1;i<arraysize;i++) -{ - - /* - ** Calculate A[i] terms. Note, once again, that we - ** can ignore the 2/period term outside the integral - ** since the period is 2 and the term cancels itself - ** out. - */ - *(abase+i)=TrapezoidIntegrate((double)0.0, - (double)2.0, - 200, - omega * (double)i, - 1); - - /* - ** Calculate the B[i] terms. - */ - *(bbase+i)=TrapezoidIntegrate((double)0.0, - (double)2.0, - 200, - omega * (double)i, - 2); - -} -#ifdef DEBUG -{ - int i; - printf("\nA[i]=\n"); - for (i=0;i<arraysize;i++) printf("%7.3g ",abase[i]); - printf("\nB[i]=\n(undefined) "); - for (i=1;i<arraysize;i++) printf("%7.3g ",bbase[i]); -} -#endif -/* -** All done, stop the stopwatch -*/ -return(StopStopwatch(elapsed)); + clock_t start, stop; + double omega; /* Fundamental frequency */ + unsigned long i; /* Index */ + + start = clock(); + + /* + ** Calculate the fourier series. Begin by + ** calculating A[0]. + */ + *abase = TrapezoidIntegrate(0.0, 2.0, 200, 0.0, 0) / 2.0; + + /* + ** Calculate the fundamental frequency. + ** ( 2 * pi ) / period...and since the period + ** is 2, omega is simply pi. + */ + omega = M_PI; + + for(i = 1; i < arraysize; i++) { + + /* + ** Calculate A[i] terms. Note, once again, that we + ** can ignore the 2/period term outside the integral + ** since the period is 2 and the term cancels itself + ** out. + */ + *(abase + i) = TrapezoidIntegrate(0.0, 2.0, 200, omega * (double)i, 1); + + /* + ** Calculate the B[i] terms. + */ + *(bbase + i) = TrapezoidIntegrate(0.0, 2.0, 200, omega * (double)i, 2); + } + + stop = clock(); + + return stop - start; } /*********************** @@ -212,37 +170,23 @@ return(StopStopwatch(elapsed)); ** x0,x1 set the lower and upper bounds of the ** integration. ** nsteps indicates # of trapezoidal sections -** omegan is the fundamental frequency times +** omega_n is the fundamental frequency times ** the series member # ** select = 0 for the A[0] term, 1 for cosine terms, and ** 2 for sine terms. ** Returns the value. -*/ -static double TrapezoidIntegrate( double x0, /* Lower bound */ - double x1, /* Upper bound */ - int nsteps, /* # of steps */ - double omegan, /* omega * n */ - int select) +** double x0 - lower bound +** double x1 - upper bound +** int nsteps - number of steps +** double omega_n - omega * n +** int select - select functions FIXME: this is dumb +*/ +static double +TrapezoidIntegrate( double x0, double x1, int nsteps, double omega_n, int select) { -double x; /* Independent variable */ -double dx; /* Stepsize */ -double rvalue; /* Return value */ - - -/* -** Initialize independent variable -*/ -x=x0; - -/* -** Calculate stepsize -*/ -dx=(x1 - x0) / (double)nsteps; - -/* -** Initialize the return value. -*/ -rvalue=thefunction(x0,omegan,select)/(double)2.0; + double x = x0; /* Independent variable */ + double dx = (x1 - x0) / (double)nsteps; /* Stepsize */ + double rvalue = thefunction(x0, omega_n, select) / 2.0; /* ** Compute the other terms of the integral. @@ -252,13 +196,13 @@ if(nsteps!=1) while(--nsteps ) { x+=dx; - rvalue+=thefunction(x,omegan,select); + rvalue+=thefunction(x,omega_n,select); } } /* ** Finish computation */ -rvalue=(rvalue+thefunction(x1,omegan,select)/(double)2.0)*dx; +rvalue=(rvalue+thefunction(x1,omega_n,select)/(double)2.0)*dx; return(rvalue); } @@ -269,12 +213,12 @@ return(rvalue); ** This routine selects the function to be used ** in the Trapezoid integration. ** x is the independent variable -** omegan is omega * n +** omega_n is omega * n ** select chooses which of the sine/cosine functions ** are used. note the special case for select=0. */ static double thefunction(double x, /* Independent variable */ - double omegan, /* Omega * term */ + double omega_n, /* Omega * term */ int select) /* Choose term */ { @@ -285,9 +229,9 @@ switch(select) { case 0: return(pow(x+(double)1.0,x)); - case 1: return(pow(x+(double)1.0,x) * cos(omegan * x)); + case 1: return(pow(x+(double)1.0,x) * cos(omega_n * x)); - case 2: return(pow(x+(double)1.0,x) * sin(omegan * x)); + case 2: return(pow(x+(double)1.0,x) * sin(omega_n * x)); } /* diff --git a/fpemulation.c b/fpemulation.c index e851cd9..61b016b 100644 --- a/fpemulation.c +++ b/fpemulation.c @@ -4,6 +4,7 @@ #include <string.h> #include <math.h> #include <limits.h> +#include <time.h> #include "nmglobal.h" #include "nbench1.h" @@ -23,18 +24,14 @@ void DoEmFloat(void) { - /* Error context string pointer */ - const char *errorcontext = "CPU:Floating Emulation"; - /* Local structure */ - EmFloatStruct *locemfloatstruct = &global_emfloatstruct; - - InternalFPF *abase = NULL; /* Base of A array */ - InternalFPF *bbase = NULL; /* Base of B array */ - InternalFPF *cbase = NULL; /* Base of C array */ - unsigned long accumtime; /* Accumulated time in ticks */ - double iterations; /* # of iterations */ - unsigned long tickcount; /* # of ticks */ - unsigned long loops; /* # of loops */ + const char* errorcontext = "CPU:Floating Emulation"; + EmFloatStruct* locemfloatstruct = &global_emfloatstruct; + InternalFPF* abase = NULL; + InternalFPF* bbase = NULL; + InternalFPF* cbase = NULL; + clock_t total_time = 0; + int iterations = 0; + unsigned long loops = 1; abase = malloc(locemfloatstruct->arraysize * sizeof(InternalFPF)); if (!abase) { @@ -61,7 +58,9 @@ DoEmFloat(void) SetupCPUEmFloatArrays(abase, bbase, cbase, locemfloatstruct->arraysize); /* FIXME: ugly */ /* See if we need to do self-adjusting code.*/ - if (locemfloatstruct->adjust == 0) { + if (locemfloatstruct->adjust == FALSE) { + locemfloatstruct->adjust = TRUE; + locemfloatstruct->loops = 0; /* @@ -69,10 +68,8 @@ DoEmFloat(void) ** less than minimum, increase the loop count and try ** again. */ - for (loops = 1; loops < CPUEMFLOATLOOPMAX; loops += loops) { - tickcount = DoEmFloatIteration(abase, bbase, cbase, /* FIXME: ugly */ - locemfloatstruct->arraysize, loops); - if (tickcount > global_min_ticks) { + for (; loops < CPUEMFLOATLOOPMAX; loops += loops) { + if (DoEmFloatIteration(abase, bbase, cbase, locemfloatstruct->arraysize, loops) > global_min_ticks) { locemfloatstruct->loops = loops; break; } @@ -96,25 +93,14 @@ DoEmFloat(void) ** # of seconds requested. ** Each iteration performs arraysize * 3 operations. */ - accumtime = 0L; - iterations = 0.0; do { - accumtime += DoEmFloatIteration(abase, bbase, cbase, /* FIXME: ugly */ - locemfloatstruct->arraysize, locemfloatstruct->loops); - iterations += 1.0; - } while (TicksToSecs(accumtime) < locemfloatstruct->request_secs); + total_time += DoEmFloatIteration(abase, bbase, cbase, locemfloatstruct->arraysize, locemfloatstruct->loops); + ++iterations; + } while (total_time < locemfloatstruct->request_secs * CLOCKS_PER_SEC); - /* - ** Clean up, calculate results, and go home. - ** Also, indicate that adjustment is done. - */ free(abase); free(bbase); free(cbase); - locemfloatstruct->emflops = (iterations * (double)locemfloatstruct->loops) - / (double)TicksToFracSecs(accumtime); - if (locemfloatstruct->adjust == 0) { - locemfloatstruct->adjust = 1; - } + locemfloatstruct->emflops = (double)(iterations * locemfloatstruct->loops * CLOCKS_PER_SEC) / (double)total_time; } |