embedded_fehlberg_7_8.c 24 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612
  1. // #ifdef __cplusplus
  2. // extern "C" {
  3. // #endif
  4. //========================================================================================================================================================================================================200
  5. // DESCRIPTION
  6. //========================================================================================================================================================================================================200
  7. ////////////////////////////////////////////////////////////////////////////////
  8. // File: embedded_fehlberg_7_8.c //
  9. // Routines: //
  10. // Embedded_Fehlberg_7_8 //
  11. ////////////////////////////////////////////////////////////////////////////////
  12. ////////////////////////////////////////////////////////////////////////////////
  13. // //
  14. // Description: //
  15. // The Runge-Kutta-Fehlberg method is an adaptive procedure for approxi- //
  16. // mating the solution of the differential equation y'(x) = f(x,y) with //
  17. // initial condition y(x0) = c. This implementation evaluates f(x,y) //
  18. // thirteen times per step using embedded seventh order and eight order //
  19. // Runge-Kutta estimates to estimate the not only the solution but also //
  20. // the error. //
  21. // The next step size is then calculated using the preassigned tolerance //
  22. // and error estimate. //
  23. // For step i+1, //
  24. // y[i+1] = y[i] + h * (41/840 * k1 + 34/105 * finavalu_temp[5] + 9/35 * finavalu_temp[6] //
  25. // + 9/35 * finavalu_temp[7] + 9/280 * finavalu_temp[8] + 9/280 finavalu_temp[9] + 41/840 finavalu_temp[10] ) //
  26. // where //
  27. // k1 = f( x[i],y[i] ), //
  28. // finavalu_temp[1] = f( x[i]+2h/27, y[i] + 2h*k1/27), //
  29. // finavalu_temp[2] = f( x[i]+h/9, y[i]+h/36*( k1 + 3 finavalu_temp[1]) ), //
  30. // finavalu_temp[3] = f( x[i]+h/6, y[i]+h/24*( k1 + 3 finavalu_temp[2]) ), //
  31. // finavalu_temp[4] = f( x[i]+5h/12, y[i]+h/48*(20 k1 - 75 finavalu_temp[2] + 75 finavalu_temp[3])), //
  32. // finavalu_temp[5] = f( x[i]+h/2, y[i]+h/20*( k1 + 5 finavalu_temp[3] + 4 finavalu_temp[4] ) ), //
  33. // finavalu_temp[6] = f( x[i]+5h/6, y[i]+h/108*( -25 k1 + 125 finavalu_temp[3] - 260 finavalu_temp[4] + 250 finavalu_temp[5] ) ), //
  34. // finavalu_temp[7] = f( x[i]+h/6, y[i]+h*( 31/300 k1 + 61/225 finavalu_temp[4] - 2/9 finavalu_temp[5] //
  35. // + 13/900 finavalu_temp[6]) ) //
  36. // finavalu_temp[8] = f( x[i]+2h/3, y[i]+h*( 2 k1 - 53/6 finavalu_temp[3] + 704/45 finavalu_temp[4] - 107/9 finavalu_temp[5] //
  37. // + 67/90 finavalu_temp[6] + 3 finavalu_temp[7]) ), //
  38. // finavalu_temp[9] = f( x[i]+h/3, y[i]+h*( -91/108 k1 + 23/108 finavalu_temp[3] - 976/135 finavalu_temp[4] //
  39. // + 311/54 finavalu_temp[5] - 19/60 finavalu_temp[6] + 17/6 finavalu_temp[7] - 1/12 finavalu_temp[8]) ), //
  40. // finavalu_temp[10] = f( x[i]+h, y[i]+h*( 2383/4100 k1 - 341/164 finavalu_temp[3] + 4496/1025 finavalu_temp[4] //
  41. // - 301/82 finavalu_temp[5] + 2133/4100 finavalu_temp[6] + 45/82 finavalu_temp[7] + 45/164 finavalu_temp[8] + 18/41 finavalu_temp[9]) ) //
  42. // finavalu_temp[11] = f( x[i], y[i]+h*( 3/205 k1 - 6/41 finavalu_temp[5] - 3/205 finavalu_temp[6] - 3/41 finavalu_temp[7] //
  43. // + 3/41 finavalu_temp[8] + 6/41 finavalu_temp[9]) ) //
  44. // finavalu_temp[12] = f( x[i]+h, y[i]+h*( -1777/4100 k1 - 341/164 finavalu_temp[3] + 4496/1025 finavalu_temp[4] //
  45. // - 289/82 finavalu_temp[5] + 2193/4100 finavalu_temp[6] + 51/82 finavalu_temp[7] + 33/164 finavalu_temp[8] + //
  46. // 12/41 finavalu_temp[9] + finavalu_temp[11]) ) //
  47. // x[i+1] = x[i] + h. //
  48. // //
  49. // The error is estimated to be //
  50. // err = -41/840 * h * ( k1 + finavalu_temp[10] - finavalu_temp[11] - finavalu_temp[12]) //
  51. // The step size h is then scaled by the scale factor //
  52. // scale = 0.8 * | epsilon * y[i] / [err * (xmax - x[0])] | ^ 1/7 //
  53. // The scale factor is further constrained 0.125 < scale < 4.0. //
  54. // The new step size is h := scale * h. //
  55. ////////////////////////////////////////////////////////////////////////////////
  56. ////////////////////////////////////////////////////////////////////////////////
  57. // static fp Runge_Kutta(fp (*f)(fp,fp), fp *y, fp x0, fp h) //
  58. // //
  59. // Description: //
  60. // This routine uses Fehlberg's embedded 7th and 8th order methods to //
  61. // approximate the solution of the differential equation y'=f(x,y) with //
  62. // the initial condition y = y[0] at x = x0. The value at x + h is //
  63. // returned in y[1]. The function returns err / h ( the absolute error //
  64. // per step size ). //
  65. // //
  66. // Arguments: //
  67. // fp *f Pointer to the function which returns the slope at (x,y) of //
  68. // integral curve of the differential equation y' = f(x,y) //
  69. // which passes through the point (x0,y[0]). //
  70. // fp y[] On input y[0] is the initial value of y at x, on output //
  71. // y[1] is the solution at x + h. //
  72. // fp x Initial value of x. //
  73. // fp h Step size //
  74. // //
  75. // Return Values: //
  76. // This routine returns the err / h. The solution of y(x) at x + h is //
  77. // returned in y[1]. //
  78. // //
  79. ////////////////////////////////////////////////////////////////////////////////
  80. //========================================================================================================================================================================================================200
  81. // DEFINE/INCLUDE
  82. //========================================================================================================================================================================================================200
  83. //======================================================================================================================================================150
  84. // KERNEL
  85. //======================================================================================================================================================150
  86. #include "master.c" // (in directory)
  87. //======================================================================================================================================================150
  88. // LIBRARIES
  89. //======================================================================================================================================================150
  90. #include <math.h> // (in path provided to compiler) needed by pow, fabs
  91. //======================================================================================================================================================150
  92. // END
  93. //======================================================================================================================================================150
  94. //========================================================================================================================================================================================================200
  95. // PARTICULAR SOLVER FUNCTION
  96. //========================================================================================================================================================================================================200
  97. void
  98. embedded_fehlberg_7_8( fp timeinst,
  99. fp h,
  100. fp *initvalu,
  101. fp *finavalu,
  102. fp *error,
  103. fp *parameter,
  104. fp *com,
  105. cl_mem d_initvalu,
  106. cl_mem d_finavalu,
  107. cl_mem d_params,
  108. cl_mem d_com,
  109. cl_command_queue command_queue,
  110. cl_kernel kernel,
  111. long long *timecopyin,
  112. long long *timecopykernel,
  113. long long *timecopyout)
  114. {
  115. //======================================================================================================================================================
  116. // VARIABLES
  117. //======================================================================================================================================================
  118. static const fp c_1_11 = 41.0 / 840.0;
  119. static const fp c6 = 34.0 / 105.0;
  120. static const fp c_7_8= 9.0 / 35.0;
  121. static const fp c_9_10 = 9.0 / 280.0;
  122. static const fp a2 = 2.0 / 27.0;
  123. static const fp a3 = 1.0 / 9.0;
  124. static const fp a4 = 1.0 / 6.0;
  125. static const fp a5 = 5.0 / 12.0;
  126. static const fp a6 = 1.0 / 2.0;
  127. static const fp a7 = 5.0 / 6.0;
  128. static const fp a8 = 1.0 / 6.0;
  129. static const fp a9 = 2.0 / 3.0;
  130. static const fp a10 = 1.0 / 3.0;
  131. static const fp b31 = 1.0 / 36.0;
  132. static const fp b32 = 3.0 / 36.0;
  133. static const fp b41 = 1.0 / 24.0;
  134. static const fp b43 = 3.0 / 24.0;
  135. static const fp b51 = 20.0 / 48.0;
  136. static const fp b53 = -75.0 / 48.0;
  137. static const fp b54 = 75.0 / 48.0;
  138. static const fp b61 = 1.0 / 20.0;
  139. static const fp b64 = 5.0 / 20.0;
  140. static const fp b65 = 4.0 / 20.0;
  141. static const fp b71 = -25.0 / 108.0;
  142. static const fp b74 = 125.0 / 108.0;
  143. static const fp b75 = -260.0 / 108.0;
  144. static const fp b76 = 250.0 / 108.0;
  145. static const fp b81 = 31.0/300.0;
  146. static const fp b85 = 61.0/225.0;
  147. static const fp b86 = -2.0/9.0;
  148. static const fp b87 = 13.0/900.0;
  149. static const fp b91 = 2.0;
  150. static const fp b94 = -53.0/6.0;
  151. static const fp b95 = 704.0 / 45.0;
  152. static const fp b96 = -107.0 / 9.0;
  153. static const fp b97 = 67.0 / 90.0;
  154. static const fp b98 = 3.0;
  155. static const fp b10_1 = -91.0 / 108.0;
  156. static const fp b10_4 = 23.0 / 108.0;
  157. static const fp b10_5 = -976.0 / 135.0;
  158. static const fp b10_6 = 311.0 / 54.0;
  159. static const fp b10_7 = -19.0 / 60.0;
  160. static const fp b10_8 = 17.0 / 6.0;
  161. static const fp b10_9 = -1.0 / 12.0;
  162. static const fp b11_1 = 2383.0 / 4100.0;
  163. static const fp b11_4 = -341.0 / 164.0;
  164. static const fp b11_5 = 4496.0 / 1025.0;
  165. static const fp b11_6 = -301.0 / 82.0;
  166. static const fp b11_7 = 2133.0 / 4100.0;
  167. static const fp b11_8 = 45.0 / 82.0;
  168. static const fp b11_9 = 45.0 / 164.0;
  169. static const fp b11_10 = 18.0 / 41.0;
  170. static const fp b12_1 = 3.0 / 205.0;
  171. static const fp b12_6 = - 6.0 / 41.0;
  172. static const fp b12_7 = - 3.0 / 205.0;
  173. static const fp b12_8 = - 3.0 / 41.0;
  174. static const fp b12_9 = 3.0 / 41.0;
  175. static const fp b12_10 = 6.0 / 41.0;
  176. static const fp b13_1 = -1777.0 / 4100.0;
  177. static const fp b13_4 = -341.0 / 164.0;
  178. static const fp b13_5 = 4496.0 / 1025.0;
  179. static const fp b13_6 = -289.0 / 82.0;
  180. static const fp b13_7 = 2193.0 / 4100.0;
  181. static const fp b13_8 = 51.0 / 82.0;
  182. static const fp b13_9 = 33.0 / 164.0;
  183. static const fp b13_10 = 12.0 / 41.0;
  184. static const fp err_factor = -41.0 / 840.0;
  185. fp h2_7 = a2 * h;
  186. fp timeinst_temp;
  187. fp* initvalu_temp;
  188. fp** finavalu_temp;
  189. int i;
  190. //======================================================================================================================================================
  191. // TEMPORARY STORAGE ALLOCATION
  192. //======================================================================================================================================================
  193. initvalu_temp= (fp *) malloc(EQUATIONS* sizeof(fp));
  194. finavalu_temp= (fp **) malloc(13* sizeof(fp *));
  195. for (i= 0; i<13; i++){
  196. finavalu_temp[i]= (fp *) malloc(EQUATIONS* sizeof(fp));
  197. }
  198. //======================================================================================================================================================
  199. // EVALUATIONS [UNROLLED LOOP] [SEQUENTIAL DEPENDENCY]
  200. //======================================================================================================================================================
  201. //===================================================================================================
  202. // 1
  203. //===================================================================================================
  204. timeinst_temp = timeinst;
  205. for(i=0; i<EQUATIONS; i++){
  206. initvalu_temp[i] = initvalu[i] ;
  207. // printf("initvalu[%d] = %f\n", i, initvalu[i]);
  208. }
  209. master( timeinst_temp,
  210. initvalu_temp,
  211. parameter,
  212. finavalu_temp[0],
  213. com,
  214. d_initvalu,
  215. d_finavalu,
  216. d_params,
  217. d_com,
  218. command_queue,
  219. kernel,
  220. timecopyin,
  221. timecopykernel,
  222. timecopyout);
  223. //===================================================================================================
  224. // 2
  225. //===================================================================================================
  226. timeinst_temp = timeinst+h2_7;
  227. for(i=0; i<EQUATIONS; i++){
  228. initvalu_temp[i] = initvalu[i] + h2_7 * (finavalu_temp[0][i]);
  229. }
  230. master( timeinst_temp,
  231. initvalu_temp,
  232. parameter,
  233. finavalu_temp[1],
  234. com,
  235. d_initvalu,
  236. d_finavalu,
  237. d_params,
  238. d_com,\
  239. command_queue,
  240. kernel,
  241. timecopyin,
  242. timecopykernel,
  243. timecopyout);
  244. //===================================================================================================
  245. // 3
  246. //===================================================================================================
  247. timeinst_temp = timeinst+a3*h;
  248. for(i=0; i<EQUATIONS; i++){
  249. initvalu_temp[i] = initvalu[i] + h * ( b31*finavalu_temp[0][i] + b32*finavalu_temp[1][i]);
  250. }
  251. master( timeinst_temp,
  252. initvalu_temp,
  253. parameter,
  254. finavalu_temp[2],
  255. com,
  256. d_initvalu,
  257. d_finavalu,
  258. d_params,
  259. d_com,
  260. command_queue,
  261. kernel,
  262. timecopyin,
  263. timecopykernel,
  264. timecopyout);
  265. //===================================================================================================
  266. // 4
  267. //===================================================================================================
  268. timeinst_temp = timeinst+a4*h;
  269. for(i=0; i<EQUATIONS; i++){
  270. initvalu_temp[i] = initvalu[i] + h * ( b41*finavalu_temp[0][i] + b43*finavalu_temp[2][i]) ;
  271. }
  272. master( timeinst_temp,
  273. initvalu_temp,
  274. parameter,
  275. finavalu_temp[3],
  276. com,
  277. d_initvalu,
  278. d_finavalu,
  279. d_params,
  280. d_com,
  281. command_queue,
  282. kernel,
  283. timecopyin,
  284. timecopykernel,
  285. timecopyout);
  286. //===================================================================================================
  287. // 5
  288. //===================================================================================================
  289. timeinst_temp = timeinst+a5*h;
  290. for(i=0; i<EQUATIONS; i++){
  291. initvalu_temp[i] = initvalu[i] + h * ( b51*finavalu_temp[0][i] + b53*finavalu_temp[2][i] + b54*finavalu_temp[3][i]) ;
  292. }
  293. master( timeinst_temp,
  294. initvalu_temp,
  295. parameter,
  296. finavalu_temp[4],
  297. com,
  298. d_initvalu,
  299. d_finavalu,
  300. d_params,
  301. d_com,
  302. command_queue,
  303. kernel,
  304. timecopyin,
  305. timecopykernel,
  306. timecopyout);
  307. //===================================================================================================
  308. // 6
  309. //===================================================================================================
  310. timeinst_temp = timeinst+a6*h;
  311. for(i=0; i<EQUATIONS; i++){
  312. initvalu_temp[i] = initvalu[i] + h * ( b61*finavalu_temp[0][i] + b64*finavalu_temp[3][i] + b65*finavalu_temp[4][i]) ;
  313. }
  314. master( timeinst_temp,
  315. initvalu_temp,
  316. parameter,
  317. finavalu_temp[5],
  318. com,
  319. d_initvalu,
  320. d_finavalu,
  321. d_params,
  322. d_com,
  323. command_queue,
  324. kernel,
  325. timecopyin,
  326. timecopykernel,
  327. timecopyout);
  328. //===================================================================================================
  329. // 7
  330. //===================================================================================================
  331. timeinst_temp = timeinst+a7*h;
  332. for(i=0; i<EQUATIONS; i++){
  333. initvalu_temp[i] = initvalu[i] + h * ( b71*finavalu_temp[0][i] + b74*finavalu_temp[3][i] + b75*finavalu_temp[4][i] + b76*finavalu_temp[5][i]);
  334. }
  335. master( timeinst_temp,
  336. initvalu_temp,
  337. parameter,
  338. finavalu_temp[6],
  339. com,
  340. d_initvalu,
  341. d_finavalu,
  342. d_params,
  343. d_com,
  344. command_queue,
  345. kernel,
  346. timecopyin,
  347. timecopykernel,
  348. timecopyout);
  349. //===================================================================================================
  350. // 8
  351. //===================================================================================================
  352. timeinst_temp = timeinst+a8*h;
  353. for(i=0; i<EQUATIONS; i++){
  354. initvalu_temp[i] = initvalu[i] + h * ( b81*finavalu_temp[0][i] + b85*finavalu_temp[4][i] + b86*finavalu_temp[5][i] + b87*finavalu_temp[6][i]);
  355. }
  356. master( timeinst_temp,
  357. initvalu_temp,
  358. parameter,
  359. finavalu_temp[7],
  360. com,
  361. d_initvalu,
  362. d_finavalu,
  363. d_params,
  364. d_com,
  365. command_queue,
  366. kernel,
  367. timecopyin,
  368. timecopykernel,
  369. timecopyout);
  370. //===================================================================================================
  371. // 9
  372. //===================================================================================================
  373. timeinst_temp = timeinst+a9*h;
  374. for(i=0; i<EQUATIONS; i++){
  375. initvalu_temp[i] = initvalu[i] + h * ( b91*finavalu_temp[0][i] + b94*finavalu_temp[3][i] + b95*finavalu_temp[4][i] + b96*finavalu_temp[5][i] + b97*finavalu_temp[6][i]+ b98*finavalu_temp[7][i]) ;
  376. }
  377. master( timeinst_temp,
  378. initvalu_temp,
  379. parameter,
  380. finavalu_temp[8],
  381. com,
  382. d_initvalu,
  383. d_finavalu,
  384. d_params,
  385. d_com,
  386. command_queue,
  387. kernel,
  388. timecopyin,
  389. timecopykernel,
  390. timecopyout);
  391. //===================================================================================================
  392. // 10
  393. //===================================================================================================
  394. timeinst_temp = timeinst+a10*h;
  395. for(i=0; i<EQUATIONS; i++){
  396. initvalu_temp[i] = initvalu[i] + h * ( b10_1*finavalu_temp[0][i] + b10_4*finavalu_temp[3][i] + b10_5*finavalu_temp[4][i] + b10_6*finavalu_temp[5][i] + b10_7*finavalu_temp[6][i] + b10_8*finavalu_temp[7][i] + b10_9*finavalu_temp[8] [i]) ;
  397. }
  398. master( timeinst_temp,
  399. initvalu_temp,
  400. parameter,
  401. finavalu_temp[9],
  402. com,
  403. d_initvalu,
  404. d_finavalu,
  405. d_params,
  406. d_com,
  407. command_queue,
  408. kernel,
  409. timecopyin,
  410. timecopykernel,
  411. timecopyout);
  412. //===================================================================================================
  413. // 11
  414. //===================================================================================================
  415. timeinst_temp = timeinst+h;
  416. for(i=0; i<EQUATIONS; i++){
  417. initvalu_temp[i] = initvalu[i] + h * ( b11_1*finavalu_temp[0][i] + b11_4*finavalu_temp[3][i] + b11_5*finavalu_temp[4][i] + b11_6*finavalu_temp[5][i] + b11_7*finavalu_temp[6][i] + b11_8*finavalu_temp[7][i] + b11_9*finavalu_temp[8][i]+ b11_10 * finavalu_temp[9][i]);
  418. }
  419. master( timeinst_temp,
  420. initvalu_temp,
  421. parameter,
  422. finavalu_temp[10],
  423. com,
  424. d_initvalu,
  425. d_finavalu,
  426. d_params,
  427. d_com,
  428. command_queue,
  429. kernel,
  430. timecopyin,
  431. timecopykernel,
  432. timecopyout);
  433. //===================================================================================================
  434. // 12
  435. //===================================================================================================
  436. timeinst_temp = timeinst;
  437. for(i=0; i<EQUATIONS; i++){
  438. initvalu_temp[i] = initvalu[i] + h * ( b12_1*finavalu_temp[0][i] + b12_6*finavalu_temp[5][i] + b12_7*finavalu_temp[6][i] + b12_8*finavalu_temp[7][i] + b12_9*finavalu_temp[8][i] + b12_10 * finavalu_temp[9][i]) ;
  439. }
  440. master( timeinst_temp,
  441. initvalu_temp,
  442. parameter,
  443. finavalu_temp[11],
  444. com,
  445. d_initvalu,
  446. d_finavalu,
  447. d_params,
  448. d_com,
  449. command_queue,
  450. kernel,
  451. timecopyin,
  452. timecopykernel,
  453. timecopyout);
  454. //===================================================================================================
  455. // 13
  456. //===================================================================================================
  457. timeinst_temp = timeinst+h;
  458. for(i=0; i<EQUATIONS; i++){
  459. initvalu_temp[i] = initvalu[i] + h * ( b13_1*finavalu_temp[0][i] + b13_4*finavalu_temp[3][i] + b13_5*finavalu_temp[4][i] + b13_6*finavalu_temp[5][i] + b13_7*finavalu_temp[6][i] + b13_8*finavalu_temp[7][i] + b13_9*finavalu_temp[8][i] + b13_10*finavalu_temp[9][i] + finavalu_temp[11][i]) ;
  460. }
  461. master( timeinst_temp,
  462. initvalu_temp,
  463. parameter,
  464. finavalu_temp[12],
  465. com,
  466. d_initvalu,
  467. d_finavalu,
  468. d_params,
  469. d_com,
  470. command_queue,
  471. kernel,
  472. timecopyin,
  473. timecopykernel,
  474. timecopyout);
  475. //======================================================================================================================================================
  476. // FINAL VALUE
  477. //======================================================================================================================================================
  478. for(i=0; i<EQUATIONS; i++){
  479. finavalu[i]= initvalu[i] + h * (c_1_11 * (finavalu_temp[0][i] + finavalu_temp[10][i]) + c6 * finavalu_temp[5][i] + c_7_8 * (finavalu_temp[6][i] + finavalu_temp[7][i]) + c_9_10 * (finavalu_temp[8][i] + finavalu_temp[9][i]) );
  480. }
  481. //======================================================================================================================================================
  482. // RETURN
  483. //======================================================================================================================================================
  484. for(i=0; i<EQUATIONS; i++){
  485. error[i] = fabs(err_factor * (finavalu_temp[0][i] + finavalu_temp[10][i] - finavalu_temp[11][i] - finavalu_temp[12][i]));
  486. }
  487. //======================================================================================================================================================
  488. // DEALLOCATION
  489. //======================================================================================================================================================
  490. free(initvalu_temp);
  491. free(finavalu_temp);
  492. }
  493. //========================================================================================================================================================================================================200
  494. // END
  495. //========================================================================================================================================================================================================200
  496. // #ifdef __cplusplus
  497. // }
  498. // #endif