// #ifdef __cplusplus // extern "C" { // #endif //========================================================================================================================================================================================================200 // DESCRIPTION //========================================================================================================================================================================================================200 //======================================================================================================================================================150 // UPDATE //======================================================================================================================================================150 // Summary of changes by Lukasz G. Szafaryn: // 1) The original code was obtained from: Mathematics Source Library (http://mymathlib.webtrellis.net/index.html) // 2) This solver and particular solving algorithm used with it (embedded_fehlberg_7_8) were adapted to work with a set of equations, not just one like in original version. // 3) In order for solver to provide deterministic number of steps (needed for particular amount of memore previousely allocated for results), every next step is incremented by 1 time unit (h_init). // 4) Function assumes that time interval starts at 0 (xmin) and ends at integer value (xmax) specified by the uses as a parameter on command line. // 5) The appropriate amount of memory is previousely allocated for that range (y). // 5) This setup in 3) - 5) allows solver to adjust the step ony from current time instance to current time instance + 0.9. The next time instance is current time instance + 1; // 6) Solver also takes parameters (params) that it then passes to the equations. // 7) The original solver cannot handle cases when equations return NAN and INF values due to discontinuities and /0. That is why equations provided by user need to make sure that no NAN and INF are returned. // Last update: 15 DEC 09 //======================================================================================================================================================150 // DESCRIPTION //======================================================================================================================================================150 // int solver( fp (*f)(fp, fp), fp y[], // // fp x, fp h, fp xmax, fp *h_next, fp tolerance ) // // // // Description: // // This function solves the differential equation y'=f(x,y) with the // // initial condition y(x) = y[0]. The value at xmax is returned in y[1]. // // The function returns 0 if successful or -1 if it fails. // // // // Arguments: // // fp *f Pointer to the function which returns the slope at (x,y) of // // integral curve of the differential equation y' = f(x,y) // // which passes through the point (x0,y0) corresponding to the // // initial condition y(x0) = y0. // // fp y[] On input y[0] is the initial value of y at x, on output // // y[1] is the solution at xmax. // // fp x The initial value of x. // // fp h Initial step size. // // fp xmax The endpoint of x. // // fp *h_next A pointer to the estimated step size for successive // // calls to solver. // // fp tolerance The tolerance of y(xmax), i.e. a solution is sought // // so that the relative error < tolerance. // // // // Return Values: // // 0 The solution of y' = f(x,y) from x to xmax is stored y[1] and // // h_next has the value to the next size to try. // // -1 The solution of y' = f(x,y) from x to xmax failed. // // -2 Failed because either xmax < x or the step size h <= 0. // // -3 Memory limit allocated for results was reached // //========================================================================================================================================================================================================200 // DEFINE/INCLUDE //========================================================================================================================================================================================================200 //======================================================================================================================================================150 // COMMON //======================================================================================================================================================150 #include "../common.h" // (in path provided here) //======================================================================================================================================================150 // DEFINE //======================================================================================================================================================150 #define max(x,y) ( (x) < (y) ? (y) : (x) ) #define min(x,y) ( (x) < (y) ? (x) : (y) ) #define ATTEMPTS 12 #define MIN_SCALE_FACTOR 0.125 #define MAX_SCALE_FACTOR 4.0 //======================================================================================================================================================150 // KERNEL //======================================================================================================================================================150 #include "./embedded_fehlberg_7_8.c" // (in path provided here) //======================================================================================================================================================150 // LIBRARIES //======================================================================================================================================================150 #include // (in path known to compiler) needed by malloc, free #include // (in path known to compiler) needed by pow, fabs //======================================================================================================================================================150 // END //======================================================================================================================================================150 //========================================================================================================================================================================================================200 // FUNCTION //========================================================================================================================================================================================================200 int solver( fp **y, fp *x, int xmax, fp *params, fp *com, cl_mem d_initvalu, cl_mem d_finavalu, cl_mem d_params, cl_mem d_com, cl_command_queue command_queue, cl_kernel kernel, long long *timecopyin, long long *timecopykernel, long long *timecopyout) { //======================================================================================================================== // VARIABLES //======================================================================================================================== // solver parameters fp err_exponent; int error; int outside; fp h; fp h_init; fp tolerance; int xmin; // memory fp scale_min; fp scale_fina; fp* err= (fp *) malloc(EQUATIONS* sizeof(fp)); fp* scale= (fp *) malloc(EQUATIONS* sizeof(fp)); fp* yy= (fp *) malloc(EQUATIONS* sizeof(fp)); // counters int i, j, k; //======================================================================================================================== // INITIAL SETUP //======================================================================================================================== // solver parameters err_exponent = 1.0 / 7.0; h_init = 1; h = h_init; xmin = 0; tolerance = 10 / (fp)(xmax-xmin); // save value for initial time instance x[0] = 0; //======================================================================================================================== // CHECKING //======================================================================================================================== // Verify that the step size is positive and that the upper endpoint of integration is greater than the initial enpoint. // if (xmax < xmin || h <= 0.0){ return -2; } // If the upper endpoint of the independent variable agrees with the initial value of the independent variable. Set the value of the dependent variable and return success. // if (xmax == xmin){ return 0; } // Insure that the step size h is not larger than the length of the integration interval. // if (h > (xmax - xmin) ) { h = (fp)xmax - (fp)xmin; } //======================================================================================================================== // SOLVING //======================================================================================================================== printf("Time Steps: "); fflush(0); for(k=1; k<=xmax; k++) { // start after initial value x[k] = k-1; h = h_init; //========================================================================================== // REINITIALIZE VARIABLES //========================================================================================== scale_fina = 1.0; //========================================================================================== // MAKE ATTEMPTS TO MINIMIZE ERROR //========================================================================================== // make attempts to minimize error for (j = 0; j < ATTEMPTS; j++) { //============================================================ // REINITIALIZE VARIABLES //============================================================ error = 0; outside = 0; scale_min = MAX_SCALE_FACTOR; //============================================================ // EVALUATE ALL EQUATIONS //============================================================ embedded_fehlberg_7_8( x[k], h, y[k-1], y[k], err, params, com, d_initvalu, d_finavalu, d_params, d_com, command_queue, kernel, timecopyin, timecopykernel, timecopyout); //============================================================ // IF THERE WAS NO ERROR FOR ANY OF EQUATIONS, SET SCALE AND LEAVE THE LOOP //============================================================ for(i=0; i 0){ error = 1; } } if (error != 1) { scale_fina = MAX_SCALE_FACTOR; break; } //============================================================ // FIGURE OUT SCALE AS THE MINIMUM OF COMPONENT SCALES //============================================================ for(i=0; i ( tolerance * yy[i] ) ){ outside = 1; } } if (outside == 0){ break; } //============================================================ // ...OTHERWISE, ADJUST STEP FOR NEXT ATTEMPT //============================================================ // scale next step in a default way h = h * scale_fina; // limit step to 0.9, because when it gets close to 1, it no longer makes sense, as 1 is already the next time instance (added to original algorithm) if (h >= 0.9) { h = 0.9; } // if instance+step exceeds range limit, limit to that range if ( x[k] + h > (fp)xmax ){ h = (fp)xmax - x[k]; } // if getting closer to range limit, decrease step else if ( x[k] + h + 0.5 * h > (fp)xmax ){ h = 0.5 * h; } } //========================================================================================== // SAVE TIME INSTANCE THAT SOLVER ENDED UP USING //========================================================================================== x[k] = x[k] + h; //========================================================================================== // IF MAXIMUM NUMBER OF ATTEMPTS REACHED AND CANNOT GIVE SOLUTION, EXIT PROGRAM WITH ERROR //========================================================================================== if ( j >= ATTEMPTS ) { return -1; } printf("%d ", k); fflush(0); } printf("\n"); fflush(0); //======================================================================================================================== // FREE MEMORY //======================================================================================================================== free(err); free(scale); free(yy); //======================================================================================================================== // FINAL RETURN //======================================================================================================================== return 0; //====================================================================================================================================================== //====================================================================================================================================================== // END OF SOLVER FUNCTION //====================================================================================================================================================== //====================================================================================================================================================== } //========================================================================================================================================================================================================200 // END //========================================================================================================================================================================================================200 // #ifdef __cplusplus // } // #endif