diff options
Diffstat (limited to 'fourier.c')
| -rw-r--r-- | fourier.c | 304 |
1 files changed, 124 insertions, 180 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)); } /* |