1 /*************************************************************************
2 ALGLIB 3.15.0 (source code generated 2019-02-20)
3 Copyright (c) Sergey Bochkanov (ALGLIB project).
4
5 >>> SOURCE LICENSE >>>
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation (www.fsf.org); either version 2 of the
9 License, or (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 A copy of the GNU General Public License is available at
17 http://www.fsf.org/licensing/licenses
18 >>> END OF LICENSE >>>
19 *************************************************************************/
20 #ifdef _MSC_VER
21 #define _CRT_SECURE_NO_WARNINGS
22 #endif
23 #include "stdafx.h"
24 #include "diffequations.h"
25
26 // disable some irrelevant warnings
27 #if (AE_COMPILER==AE_MSVC) && !defined(AE_ALL_WARNINGS)
28 #pragma warning(disable:4100)
29 #pragma warning(disable:4127)
30 #pragma warning(disable:4611)
31 #pragma warning(disable:4702)
32 #pragma warning(disable:4996)
33 #endif
34
35 /////////////////////////////////////////////////////////////////////////
36 //
37 // THIS SECTION CONTAINS IMPLEMENTATION OF C++ INTERFACE
38 //
39 /////////////////////////////////////////////////////////////////////////
40 namespace alglib
41 {
42
43 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD)
44
45 #endif
46
47 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD)
48 /*************************************************************************
49
50 *************************************************************************/
_odesolverstate_owner()51 _odesolverstate_owner::_odesolverstate_owner()
52 {
53 jmp_buf _break_jump;
54 alglib_impl::ae_state _state;
55
56 alglib_impl::ae_state_init(&_state);
57 if( setjmp(_break_jump) )
58 {
59 if( p_struct!=NULL )
60 {
61 alglib_impl::_odesolverstate_destroy(p_struct);
62 alglib_impl::ae_free(p_struct);
63 }
64 p_struct = NULL;
65 #if !defined(AE_NO_EXCEPTIONS)
66 _ALGLIB_CPP_EXCEPTION(_state.error_msg);
67 #else
68 _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
69 return;
70 #endif
71 }
72 alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
73 p_struct = NULL;
74 p_struct = (alglib_impl::odesolverstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::odesolverstate), &_state);
75 memset(p_struct, 0, sizeof(alglib_impl::odesolverstate));
76 alglib_impl::_odesolverstate_init(p_struct, &_state, ae_false);
77 ae_state_clear(&_state);
78 }
79
_odesolverstate_owner(const _odesolverstate_owner & rhs)80 _odesolverstate_owner::_odesolverstate_owner(const _odesolverstate_owner &rhs)
81 {
82 jmp_buf _break_jump;
83 alglib_impl::ae_state _state;
84
85 alglib_impl::ae_state_init(&_state);
86 if( setjmp(_break_jump) )
87 {
88 if( p_struct!=NULL )
89 {
90 alglib_impl::_odesolverstate_destroy(p_struct);
91 alglib_impl::ae_free(p_struct);
92 }
93 p_struct = NULL;
94 #if !defined(AE_NO_EXCEPTIONS)
95 _ALGLIB_CPP_EXCEPTION(_state.error_msg);
96 #else
97 _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
98 return;
99 #endif
100 }
101 alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
102 p_struct = NULL;
103 alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: odesolverstate copy constructor failure (source is not initialized)", &_state);
104 p_struct = (alglib_impl::odesolverstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::odesolverstate), &_state);
105 memset(p_struct, 0, sizeof(alglib_impl::odesolverstate));
106 alglib_impl::_odesolverstate_init_copy(p_struct, const_cast<alglib_impl::odesolverstate*>(rhs.p_struct), &_state, ae_false);
107 ae_state_clear(&_state);
108 }
109
operator =(const _odesolverstate_owner & rhs)110 _odesolverstate_owner& _odesolverstate_owner::operator=(const _odesolverstate_owner &rhs)
111 {
112 if( this==&rhs )
113 return *this;
114 jmp_buf _break_jump;
115 alglib_impl::ae_state _state;
116
117 alglib_impl::ae_state_init(&_state);
118 if( setjmp(_break_jump) )
119 {
120 #if !defined(AE_NO_EXCEPTIONS)
121 _ALGLIB_CPP_EXCEPTION(_state.error_msg);
122 #else
123 _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
124 return *this;
125 #endif
126 }
127 alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
128 alglib_impl::ae_assert(p_struct!=NULL, "ALGLIB: odesolverstate assignment constructor failure (destination is not initialized)", &_state);
129 alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: odesolverstate assignment constructor failure (source is not initialized)", &_state);
130 alglib_impl::_odesolverstate_destroy(p_struct);
131 memset(p_struct, 0, sizeof(alglib_impl::odesolverstate));
132 alglib_impl::_odesolverstate_init_copy(p_struct, const_cast<alglib_impl::odesolverstate*>(rhs.p_struct), &_state, ae_false);
133 ae_state_clear(&_state);
134 return *this;
135 }
136
~_odesolverstate_owner()137 _odesolverstate_owner::~_odesolverstate_owner()
138 {
139 if( p_struct!=NULL )
140 {
141 alglib_impl::_odesolverstate_destroy(p_struct);
142 ae_free(p_struct);
143 }
144 }
145
c_ptr()146 alglib_impl::odesolverstate* _odesolverstate_owner::c_ptr()
147 {
148 return p_struct;
149 }
150
c_ptr() const151 alglib_impl::odesolverstate* _odesolverstate_owner::c_ptr() const
152 {
153 return const_cast<alglib_impl::odesolverstate*>(p_struct);
154 }
odesolverstate()155 odesolverstate::odesolverstate() : _odesolverstate_owner() ,needdy(p_struct->needdy),y(&p_struct->y),dy(&p_struct->dy),x(p_struct->x)
156 {
157 }
158
odesolverstate(const odesolverstate & rhs)159 odesolverstate::odesolverstate(const odesolverstate &rhs):_odesolverstate_owner(rhs) ,needdy(p_struct->needdy),y(&p_struct->y),dy(&p_struct->dy),x(p_struct->x)
160 {
161 }
162
operator =(const odesolverstate & rhs)163 odesolverstate& odesolverstate::operator=(const odesolverstate &rhs)
164 {
165 if( this==&rhs )
166 return *this;
167 _odesolverstate_owner::operator=(rhs);
168 return *this;
169 }
170
~odesolverstate()171 odesolverstate::~odesolverstate()
172 {
173 }
174
175
176 /*************************************************************************
177
178 *************************************************************************/
_odesolverreport_owner()179 _odesolverreport_owner::_odesolverreport_owner()
180 {
181 jmp_buf _break_jump;
182 alglib_impl::ae_state _state;
183
184 alglib_impl::ae_state_init(&_state);
185 if( setjmp(_break_jump) )
186 {
187 if( p_struct!=NULL )
188 {
189 alglib_impl::_odesolverreport_destroy(p_struct);
190 alglib_impl::ae_free(p_struct);
191 }
192 p_struct = NULL;
193 #if !defined(AE_NO_EXCEPTIONS)
194 _ALGLIB_CPP_EXCEPTION(_state.error_msg);
195 #else
196 _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
197 return;
198 #endif
199 }
200 alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
201 p_struct = NULL;
202 p_struct = (alglib_impl::odesolverreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::odesolverreport), &_state);
203 memset(p_struct, 0, sizeof(alglib_impl::odesolverreport));
204 alglib_impl::_odesolverreport_init(p_struct, &_state, ae_false);
205 ae_state_clear(&_state);
206 }
207
_odesolverreport_owner(const _odesolverreport_owner & rhs)208 _odesolverreport_owner::_odesolverreport_owner(const _odesolverreport_owner &rhs)
209 {
210 jmp_buf _break_jump;
211 alglib_impl::ae_state _state;
212
213 alglib_impl::ae_state_init(&_state);
214 if( setjmp(_break_jump) )
215 {
216 if( p_struct!=NULL )
217 {
218 alglib_impl::_odesolverreport_destroy(p_struct);
219 alglib_impl::ae_free(p_struct);
220 }
221 p_struct = NULL;
222 #if !defined(AE_NO_EXCEPTIONS)
223 _ALGLIB_CPP_EXCEPTION(_state.error_msg);
224 #else
225 _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
226 return;
227 #endif
228 }
229 alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
230 p_struct = NULL;
231 alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: odesolverreport copy constructor failure (source is not initialized)", &_state);
232 p_struct = (alglib_impl::odesolverreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::odesolverreport), &_state);
233 memset(p_struct, 0, sizeof(alglib_impl::odesolverreport));
234 alglib_impl::_odesolverreport_init_copy(p_struct, const_cast<alglib_impl::odesolverreport*>(rhs.p_struct), &_state, ae_false);
235 ae_state_clear(&_state);
236 }
237
operator =(const _odesolverreport_owner & rhs)238 _odesolverreport_owner& _odesolverreport_owner::operator=(const _odesolverreport_owner &rhs)
239 {
240 if( this==&rhs )
241 return *this;
242 jmp_buf _break_jump;
243 alglib_impl::ae_state _state;
244
245 alglib_impl::ae_state_init(&_state);
246 if( setjmp(_break_jump) )
247 {
248 #if !defined(AE_NO_EXCEPTIONS)
249 _ALGLIB_CPP_EXCEPTION(_state.error_msg);
250 #else
251 _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
252 return *this;
253 #endif
254 }
255 alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
256 alglib_impl::ae_assert(p_struct!=NULL, "ALGLIB: odesolverreport assignment constructor failure (destination is not initialized)", &_state);
257 alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: odesolverreport assignment constructor failure (source is not initialized)", &_state);
258 alglib_impl::_odesolverreport_destroy(p_struct);
259 memset(p_struct, 0, sizeof(alglib_impl::odesolverreport));
260 alglib_impl::_odesolverreport_init_copy(p_struct, const_cast<alglib_impl::odesolverreport*>(rhs.p_struct), &_state, ae_false);
261 ae_state_clear(&_state);
262 return *this;
263 }
264
~_odesolverreport_owner()265 _odesolverreport_owner::~_odesolverreport_owner()
266 {
267 if( p_struct!=NULL )
268 {
269 alglib_impl::_odesolverreport_destroy(p_struct);
270 ae_free(p_struct);
271 }
272 }
273
c_ptr()274 alglib_impl::odesolverreport* _odesolverreport_owner::c_ptr()
275 {
276 return p_struct;
277 }
278
c_ptr() const279 alglib_impl::odesolverreport* _odesolverreport_owner::c_ptr() const
280 {
281 return const_cast<alglib_impl::odesolverreport*>(p_struct);
282 }
odesolverreport()283 odesolverreport::odesolverreport() : _odesolverreport_owner() ,nfev(p_struct->nfev),terminationtype(p_struct->terminationtype)
284 {
285 }
286
odesolverreport(const odesolverreport & rhs)287 odesolverreport::odesolverreport(const odesolverreport &rhs):_odesolverreport_owner(rhs) ,nfev(p_struct->nfev),terminationtype(p_struct->terminationtype)
288 {
289 }
290
operator =(const odesolverreport & rhs)291 odesolverreport& odesolverreport::operator=(const odesolverreport &rhs)
292 {
293 if( this==&rhs )
294 return *this;
295 _odesolverreport_owner::operator=(rhs);
296 return *this;
297 }
298
~odesolverreport()299 odesolverreport::~odesolverreport()
300 {
301 }
302
303 /*************************************************************************
304 Cash-Karp adaptive ODE solver.
305
306 This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys
307 (here Y may be single variable or vector of N variables).
308
309 INPUT PARAMETERS:
310 Y - initial conditions, array[0..N-1].
311 contains values of Y[] at X[0]
312 N - system size
313 X - points at which Y should be tabulated, array[0..M-1]
314 integrations starts at X[0], ends at X[M-1], intermediate
315 values at X[i] are returned too.
316 SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!
317 M - number of intermediate points + first point + last point:
318 * M>2 means that you need both Y(X[M-1]) and M-2 values at
319 intermediate points
320 * M=2 means that you want just to integrate from X[0] to
321 X[1] and don't interested in intermediate values.
322 * M=1 means that you don't want to integrate :)
323 it is degenerate case, but it will be handled correctly.
324 * M<1 means error
325 Eps - tolerance (absolute/relative error on each step will be
326 less than Eps). When passing:
327 * Eps>0, it means desired ABSOLUTE error
328 * Eps<0, it means desired RELATIVE error. Relative errors
329 are calculated with respect to maximum values of Y seen
330 so far. Be careful to use this criterion when starting
331 from Y[] that are close to zero.
332 H - initial step lenth, it will be adjusted automatically
333 after the first step. If H=0, step will be selected
334 automatically (usualy it will be equal to 0.001 of
335 min(x[i]-x[j])).
336
337 OUTPUT PARAMETERS
338 State - structure which stores algorithm state between subsequent
339 calls of OdeSolverIteration. Used for reverse communication.
340 This structure should be passed to the OdeSolverIteration
341 subroutine.
342
343 SEE ALSO
344 AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
345
346
347 -- ALGLIB --
348 Copyright 01.09.2009 by Bochkanov Sergey
349 *************************************************************************/
odesolverrkck(const real_1d_array & y,const ae_int_t n,const real_1d_array & x,const ae_int_t m,const double eps,const double h,odesolverstate & state,const xparams _xparams)350 void odesolverrkck(const real_1d_array &y, const ae_int_t n, const real_1d_array &x, const ae_int_t m, const double eps, const double h, odesolverstate &state, const xparams _xparams)
351 {
352 jmp_buf _break_jump;
353 alglib_impl::ae_state _alglib_env_state;
354 alglib_impl::ae_state_init(&_alglib_env_state);
355 if( setjmp(_break_jump) )
356 {
357 #if !defined(AE_NO_EXCEPTIONS)
358 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
359 #else
360 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
361 return;
362 #endif
363 }
364 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
365 if( _xparams.flags!=0x0 )
366 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
367 alglib_impl::odesolverrkck(const_cast<alglib_impl::ae_vector*>(y.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), m, eps, h, const_cast<alglib_impl::odesolverstate*>(state.c_ptr()), &_alglib_env_state);
368 alglib_impl::ae_state_clear(&_alglib_env_state);
369 return;
370 }
371
372 /*************************************************************************
373 Cash-Karp adaptive ODE solver.
374
375 This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys
376 (here Y may be single variable or vector of N variables).
377
378 INPUT PARAMETERS:
379 Y - initial conditions, array[0..N-1].
380 contains values of Y[] at X[0]
381 N - system size
382 X - points at which Y should be tabulated, array[0..M-1]
383 integrations starts at X[0], ends at X[M-1], intermediate
384 values at X[i] are returned too.
385 SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!
386 M - number of intermediate points + first point + last point:
387 * M>2 means that you need both Y(X[M-1]) and M-2 values at
388 intermediate points
389 * M=2 means that you want just to integrate from X[0] to
390 X[1] and don't interested in intermediate values.
391 * M=1 means that you don't want to integrate :)
392 it is degenerate case, but it will be handled correctly.
393 * M<1 means error
394 Eps - tolerance (absolute/relative error on each step will be
395 less than Eps). When passing:
396 * Eps>0, it means desired ABSOLUTE error
397 * Eps<0, it means desired RELATIVE error. Relative errors
398 are calculated with respect to maximum values of Y seen
399 so far. Be careful to use this criterion when starting
400 from Y[] that are close to zero.
401 H - initial step lenth, it will be adjusted automatically
402 after the first step. If H=0, step will be selected
403 automatically (usualy it will be equal to 0.001 of
404 min(x[i]-x[j])).
405
406 OUTPUT PARAMETERS
407 State - structure which stores algorithm state between subsequent
408 calls of OdeSolverIteration. Used for reverse communication.
409 This structure should be passed to the OdeSolverIteration
410 subroutine.
411
412 SEE ALSO
413 AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
414
415
416 -- ALGLIB --
417 Copyright 01.09.2009 by Bochkanov Sergey
418 *************************************************************************/
419 #if !defined(AE_NO_EXCEPTIONS)
odesolverrkck(const real_1d_array & y,const real_1d_array & x,const double eps,const double h,odesolverstate & state,const xparams _xparams)420 void odesolverrkck(const real_1d_array &y, const real_1d_array &x, const double eps, const double h, odesolverstate &state, const xparams _xparams)
421 {
422 jmp_buf _break_jump;
423 alglib_impl::ae_state _alglib_env_state;
424 ae_int_t n;
425 ae_int_t m;
426
427 n = y.length();
428 m = x.length();
429 alglib_impl::ae_state_init(&_alglib_env_state);
430 if( setjmp(_break_jump) )
431 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
432 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
433 if( _xparams.flags!=0x0 )
434 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
435 alglib_impl::odesolverrkck(const_cast<alglib_impl::ae_vector*>(y.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), m, eps, h, const_cast<alglib_impl::odesolverstate*>(state.c_ptr()), &_alglib_env_state);
436
437 alglib_impl::ae_state_clear(&_alglib_env_state);
438 return;
439 }
440 #endif
441
442 /*************************************************************************
443 This function provides reverse communication interface
444 Reverse communication interface is not documented or recommended to use.
445 See below for functions which provide better documented API
446 *************************************************************************/
odesolveriteration(const odesolverstate & state,const xparams _xparams)447 bool odesolveriteration(const odesolverstate &state, const xparams _xparams)
448 {
449 jmp_buf _break_jump;
450 alglib_impl::ae_state _alglib_env_state;
451 alglib_impl::ae_state_init(&_alglib_env_state);
452 if( setjmp(_break_jump) )
453 {
454 #if !defined(AE_NO_EXCEPTIONS)
455 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
456 #else
457 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
458 return 0;
459 #endif
460 }
461 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
462 if( _xparams.flags!=0x0 )
463 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
464 ae_bool result = alglib_impl::odesolveriteration(const_cast<alglib_impl::odesolverstate*>(state.c_ptr()), &_alglib_env_state);
465 alglib_impl::ae_state_clear(&_alglib_env_state);
466 return *(reinterpret_cast<bool*>(&result));
467 }
468
469
odesolversolve(odesolverstate & state,void (* diff)(const real_1d_array & y,double x,real_1d_array & dy,void * ptr),void * ptr,const xparams _xparams)470 void odesolversolve(odesolverstate &state,
471 void (*diff)(const real_1d_array &y, double x, real_1d_array &dy, void *ptr),
472 void *ptr, const xparams _xparams){
473 jmp_buf _break_jump;
474 alglib_impl::ae_state _alglib_env_state;
475 alglib_impl::ae_state_init(&_alglib_env_state);
476 if( setjmp(_break_jump) )
477 {
478 #if !defined(AE_NO_EXCEPTIONS)
479 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
480 #else
481 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
482 return;
483 #endif
484 }
485 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
486 if( _xparams.flags!=0x0 )
487 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
488 alglib_impl::ae_assert(diff!=NULL, "ALGLIB: error in 'odesolversolve()' (diff is NULL)", &_alglib_env_state);
489 while( alglib_impl::odesolveriteration(state.c_ptr(), &_alglib_env_state) )
490 {
491 _ALGLIB_CALLBACK_EXCEPTION_GUARD_BEGIN
492 if( state.needdy )
493 {
494 diff(state.y, state.x, state.dy, ptr);
495 continue;
496 }
497 goto lbl_no_callback;
498 _ALGLIB_CALLBACK_EXCEPTION_GUARD_END
499 lbl_no_callback:
500 alglib_impl::ae_assert(ae_false, "ALGLIB: unexpected error in 'odesolversolve'", &_alglib_env_state);
501 lbl_user_exception:
502 alglib_impl::ae_assert(ae_false, "ALGLIB: exception generated in user callback", &_alglib_env_state);
503 }
504 alglib_impl::ae_state_clear(&_alglib_env_state);
505 }
506
507
508
509 /*************************************************************************
510 ODE solver results
511
512 Called after OdeSolverIteration returned False.
513
514 INPUT PARAMETERS:
515 State - algorithm state (used by OdeSolverIteration).
516
517 OUTPUT PARAMETERS:
518 M - number of tabulated values, M>=1
519 XTbl - array[0..M-1], values of X
520 YTbl - array[0..M-1,0..N-1], values of Y in X[i]
521 Rep - solver report:
522 * Rep.TerminationType completetion code:
523 * -2 X is not ordered by ascending/descending or
524 there are non-distinct X[], i.e. X[i]=X[i+1]
525 * -1 incorrect parameters were specified
526 * 1 task has been solved
527 * Rep.NFEV contains number of function calculations
528
529 -- ALGLIB --
530 Copyright 01.09.2009 by Bochkanov Sergey
531 *************************************************************************/
odesolverresults(const odesolverstate & state,ae_int_t & m,real_1d_array & xtbl,real_2d_array & ytbl,odesolverreport & rep,const xparams _xparams)532 void odesolverresults(const odesolverstate &state, ae_int_t &m, real_1d_array &xtbl, real_2d_array &ytbl, odesolverreport &rep, const xparams _xparams)
533 {
534 jmp_buf _break_jump;
535 alglib_impl::ae_state _alglib_env_state;
536 alglib_impl::ae_state_init(&_alglib_env_state);
537 if( setjmp(_break_jump) )
538 {
539 #if !defined(AE_NO_EXCEPTIONS)
540 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
541 #else
542 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
543 return;
544 #endif
545 }
546 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
547 if( _xparams.flags!=0x0 )
548 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
549 alglib_impl::odesolverresults(const_cast<alglib_impl::odesolverstate*>(state.c_ptr()), &m, const_cast<alglib_impl::ae_vector*>(xtbl.c_ptr()), const_cast<alglib_impl::ae_matrix*>(ytbl.c_ptr()), const_cast<alglib_impl::odesolverreport*>(rep.c_ptr()), &_alglib_env_state);
550 alglib_impl::ae_state_clear(&_alglib_env_state);
551 return;
552 }
553 #endif
554 }
555
556 /////////////////////////////////////////////////////////////////////////
557 //
558 // THIS SECTION CONTAINS IMPLEMENTATION OF COMPUTATIONAL CORE
559 //
560 /////////////////////////////////////////////////////////////////////////
561 namespace alglib_impl
562 {
563 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD)
564 static double odesolver_odesolvermaxgrow = 3.0;
565 static double odesolver_odesolvermaxshrink = 10.0;
566 static void odesolver_odesolverinit(ae_int_t solvertype,
567 /* Real */ ae_vector* y,
568 ae_int_t n,
569 /* Real */ ae_vector* x,
570 ae_int_t m,
571 double eps,
572 double h,
573 odesolverstate* state,
574 ae_state *_state);
575
576
577 #endif
578
579 #if defined(AE_COMPILE_ODESOLVER) || !defined(AE_PARTIAL_BUILD)
580
581
582 /*************************************************************************
583 Cash-Karp adaptive ODE solver.
584
585 This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys
586 (here Y may be single variable or vector of N variables).
587
588 INPUT PARAMETERS:
589 Y - initial conditions, array[0..N-1].
590 contains values of Y[] at X[0]
591 N - system size
592 X - points at which Y should be tabulated, array[0..M-1]
593 integrations starts at X[0], ends at X[M-1], intermediate
594 values at X[i] are returned too.
595 SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!
596 M - number of intermediate points + first point + last point:
597 * M>2 means that you need both Y(X[M-1]) and M-2 values at
598 intermediate points
599 * M=2 means that you want just to integrate from X[0] to
600 X[1] and don't interested in intermediate values.
601 * M=1 means that you don't want to integrate :)
602 it is degenerate case, but it will be handled correctly.
603 * M<1 means error
604 Eps - tolerance (absolute/relative error on each step will be
605 less than Eps). When passing:
606 * Eps>0, it means desired ABSOLUTE error
607 * Eps<0, it means desired RELATIVE error. Relative errors
608 are calculated with respect to maximum values of Y seen
609 so far. Be careful to use this criterion when starting
610 from Y[] that are close to zero.
611 H - initial step lenth, it will be adjusted automatically
612 after the first step. If H=0, step will be selected
613 automatically (usualy it will be equal to 0.001 of
614 min(x[i]-x[j])).
615
616 OUTPUT PARAMETERS
617 State - structure which stores algorithm state between subsequent
618 calls of OdeSolverIteration. Used for reverse communication.
619 This structure should be passed to the OdeSolverIteration
620 subroutine.
621
622 SEE ALSO
623 AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.
624
625
626 -- ALGLIB --
627 Copyright 01.09.2009 by Bochkanov Sergey
628 *************************************************************************/
odesolverrkck(ae_vector * y,ae_int_t n,ae_vector * x,ae_int_t m,double eps,double h,odesolverstate * state,ae_state * _state)629 void odesolverrkck(/* Real */ ae_vector* y,
630 ae_int_t n,
631 /* Real */ ae_vector* x,
632 ae_int_t m,
633 double eps,
634 double h,
635 odesolverstate* state,
636 ae_state *_state)
637 {
638
639 _odesolverstate_clear(state);
640
641 ae_assert(n>=1, "ODESolverRKCK: N<1!", _state);
642 ae_assert(m>=1, "ODESolverRKCK: M<1!", _state);
643 ae_assert(y->cnt>=n, "ODESolverRKCK: Length(Y)<N!", _state);
644 ae_assert(x->cnt>=m, "ODESolverRKCK: Length(X)<M!", _state);
645 ae_assert(isfinitevector(y, n, _state), "ODESolverRKCK: Y contains infinite or NaN values!", _state);
646 ae_assert(isfinitevector(x, m, _state), "ODESolverRKCK: Y contains infinite or NaN values!", _state);
647 ae_assert(ae_isfinite(eps, _state), "ODESolverRKCK: Eps is not finite!", _state);
648 ae_assert(ae_fp_neq(eps,(double)(0)), "ODESolverRKCK: Eps is zero!", _state);
649 ae_assert(ae_isfinite(h, _state), "ODESolverRKCK: H is not finite!", _state);
650 odesolver_odesolverinit(0, y, n, x, m, eps, h, state, _state);
651 }
652
653
654 /*************************************************************************
655
656 -- ALGLIB --
657 Copyright 01.09.2009 by Bochkanov Sergey
658 *************************************************************************/
odesolveriteration(odesolverstate * state,ae_state * _state)659 ae_bool odesolveriteration(odesolverstate* state, ae_state *_state)
660 {
661 ae_int_t n;
662 ae_int_t m;
663 ae_int_t i;
664 ae_int_t j;
665 ae_int_t k;
666 double xc;
667 double v;
668 double h;
669 double h2;
670 ae_bool gridpoint;
671 double err;
672 double maxgrowpow;
673 ae_int_t klimit;
674 ae_bool result;
675
676
677
678 /*
679 * Reverse communication preparations
680 * I know it looks ugly, but it works the same way
681 * anywhere from C++ to Python.
682 *
683 * This code initializes locals by:
684 * * random values determined during code
685 * generation - on first subroutine call
686 * * values from previous call - on subsequent calls
687 */
688 if( state->rstate.stage>=0 )
689 {
690 n = state->rstate.ia.ptr.p_int[0];
691 m = state->rstate.ia.ptr.p_int[1];
692 i = state->rstate.ia.ptr.p_int[2];
693 j = state->rstate.ia.ptr.p_int[3];
694 k = state->rstate.ia.ptr.p_int[4];
695 klimit = state->rstate.ia.ptr.p_int[5];
696 gridpoint = state->rstate.ba.ptr.p_bool[0];
697 xc = state->rstate.ra.ptr.p_double[0];
698 v = state->rstate.ra.ptr.p_double[1];
699 h = state->rstate.ra.ptr.p_double[2];
700 h2 = state->rstate.ra.ptr.p_double[3];
701 err = state->rstate.ra.ptr.p_double[4];
702 maxgrowpow = state->rstate.ra.ptr.p_double[5];
703 }
704 else
705 {
706 n = 359;
707 m = -58;
708 i = -919;
709 j = -909;
710 k = 81;
711 klimit = 255;
712 gridpoint = ae_false;
713 xc = -788;
714 v = 809;
715 h = 205;
716 h2 = -838;
717 err = 939;
718 maxgrowpow = -526;
719 }
720 if( state->rstate.stage==0 )
721 {
722 goto lbl_0;
723 }
724
725 /*
726 * Routine body
727 */
728
729 /*
730 * prepare
731 */
732 if( state->repterminationtype!=0 )
733 {
734 result = ae_false;
735 return result;
736 }
737 n = state->n;
738 m = state->m;
739 h = state->h;
740 maxgrowpow = ae_pow(odesolver_odesolvermaxgrow, (double)(5), _state);
741 state->repnfev = 0;
742
743 /*
744 * some preliminary checks for internal errors
745 * after this we assume that H>0 and M>1
746 */
747 ae_assert(ae_fp_greater(state->h,(double)(0)), "ODESolver: internal error", _state);
748 ae_assert(m>1, "ODESolverIteration: internal error", _state);
749
750 /*
751 * choose solver
752 */
753 if( state->solvertype!=0 )
754 {
755 goto lbl_1;
756 }
757
758 /*
759 * Cask-Karp solver
760 * Prepare coefficients table.
761 * Check it for errors
762 */
763 ae_vector_set_length(&state->rka, 6, _state);
764 state->rka.ptr.p_double[0] = (double)(0);
765 state->rka.ptr.p_double[1] = (double)1/(double)5;
766 state->rka.ptr.p_double[2] = (double)3/(double)10;
767 state->rka.ptr.p_double[3] = (double)3/(double)5;
768 state->rka.ptr.p_double[4] = (double)(1);
769 state->rka.ptr.p_double[5] = (double)7/(double)8;
770 ae_matrix_set_length(&state->rkb, 6, 5, _state);
771 state->rkb.ptr.pp_double[1][0] = (double)1/(double)5;
772 state->rkb.ptr.pp_double[2][0] = (double)3/(double)40;
773 state->rkb.ptr.pp_double[2][1] = (double)9/(double)40;
774 state->rkb.ptr.pp_double[3][0] = (double)3/(double)10;
775 state->rkb.ptr.pp_double[3][1] = -(double)9/(double)10;
776 state->rkb.ptr.pp_double[3][2] = (double)6/(double)5;
777 state->rkb.ptr.pp_double[4][0] = -(double)11/(double)54;
778 state->rkb.ptr.pp_double[4][1] = (double)5/(double)2;
779 state->rkb.ptr.pp_double[4][2] = -(double)70/(double)27;
780 state->rkb.ptr.pp_double[4][3] = (double)35/(double)27;
781 state->rkb.ptr.pp_double[5][0] = (double)1631/(double)55296;
782 state->rkb.ptr.pp_double[5][1] = (double)175/(double)512;
783 state->rkb.ptr.pp_double[5][2] = (double)575/(double)13824;
784 state->rkb.ptr.pp_double[5][3] = (double)44275/(double)110592;
785 state->rkb.ptr.pp_double[5][4] = (double)253/(double)4096;
786 ae_vector_set_length(&state->rkc, 6, _state);
787 state->rkc.ptr.p_double[0] = (double)37/(double)378;
788 state->rkc.ptr.p_double[1] = (double)(0);
789 state->rkc.ptr.p_double[2] = (double)250/(double)621;
790 state->rkc.ptr.p_double[3] = (double)125/(double)594;
791 state->rkc.ptr.p_double[4] = (double)(0);
792 state->rkc.ptr.p_double[5] = (double)512/(double)1771;
793 ae_vector_set_length(&state->rkcs, 6, _state);
794 state->rkcs.ptr.p_double[0] = (double)2825/(double)27648;
795 state->rkcs.ptr.p_double[1] = (double)(0);
796 state->rkcs.ptr.p_double[2] = (double)18575/(double)48384;
797 state->rkcs.ptr.p_double[3] = (double)13525/(double)55296;
798 state->rkcs.ptr.p_double[4] = (double)277/(double)14336;
799 state->rkcs.ptr.p_double[5] = (double)1/(double)4;
800 ae_matrix_set_length(&state->rkk, 6, n, _state);
801
802 /*
803 * Main cycle consists of two iterations:
804 * * outer where we travel from X[i-1] to X[i]
805 * * inner where we travel inside [X[i-1],X[i]]
806 */
807 ae_matrix_set_length(&state->ytbl, m, n, _state);
808 ae_vector_set_length(&state->escale, n, _state);
809 ae_vector_set_length(&state->yn, n, _state);
810 ae_vector_set_length(&state->yns, n, _state);
811 xc = state->xg.ptr.p_double[0];
812 ae_v_move(&state->ytbl.ptr.pp_double[0][0], 1, &state->yc.ptr.p_double[0], 1, ae_v_len(0,n-1));
813 for(j=0; j<=n-1; j++)
814 {
815 state->escale.ptr.p_double[j] = (double)(0);
816 }
817 i = 1;
818 lbl_3:
819 if( i>m-1 )
820 {
821 goto lbl_5;
822 }
823
824 /*
825 * begin inner iteration
826 */
827 lbl_6:
828 if( ae_false )
829 {
830 goto lbl_7;
831 }
832
833 /*
834 * truncate step if needed (beyond right boundary).
835 * determine should we store X or not
836 */
837 if( ae_fp_greater_eq(xc+h,state->xg.ptr.p_double[i]) )
838 {
839 h = state->xg.ptr.p_double[i]-xc;
840 gridpoint = ae_true;
841 }
842 else
843 {
844 gridpoint = ae_false;
845 }
846
847 /*
848 * Update error scale maximums
849 *
850 * These maximums are initialized by zeros,
851 * then updated every iterations.
852 */
853 for(j=0; j<=n-1; j++)
854 {
855 state->escale.ptr.p_double[j] = ae_maxreal(state->escale.ptr.p_double[j], ae_fabs(state->yc.ptr.p_double[j], _state), _state);
856 }
857
858 /*
859 * make one step:
860 * 1. calculate all info needed to do step
861 * 2. update errors scale maximums using values/derivatives
862 * obtained during (1)
863 *
864 * Take into account that we use scaling of X to reduce task
865 * to the form where x[0] < x[1] < ... < x[n-1]. So X is
866 * replaced by x=xscale*t, and dy/dx=f(y,x) is replaced
867 * by dy/dt=xscale*f(y,xscale*t).
868 */
869 ae_v_move(&state->yn.ptr.p_double[0], 1, &state->yc.ptr.p_double[0], 1, ae_v_len(0,n-1));
870 ae_v_move(&state->yns.ptr.p_double[0], 1, &state->yc.ptr.p_double[0], 1, ae_v_len(0,n-1));
871 k = 0;
872 lbl_8:
873 if( k>5 )
874 {
875 goto lbl_10;
876 }
877
878 /*
879 * prepare data for the next update of YN/YNS
880 */
881 state->x = state->xscale*(xc+state->rka.ptr.p_double[k]*h);
882 ae_v_move(&state->y.ptr.p_double[0], 1, &state->yc.ptr.p_double[0], 1, ae_v_len(0,n-1));
883 for(j=0; j<=k-1; j++)
884 {
885 v = state->rkb.ptr.pp_double[k][j];
886 ae_v_addd(&state->y.ptr.p_double[0], 1, &state->rkk.ptr.pp_double[j][0], 1, ae_v_len(0,n-1), v);
887 }
888 state->needdy = ae_true;
889 state->rstate.stage = 0;
890 goto lbl_rcomm;
891 lbl_0:
892 state->needdy = ae_false;
893 state->repnfev = state->repnfev+1;
894 v = h*state->xscale;
895 ae_v_moved(&state->rkk.ptr.pp_double[k][0], 1, &state->dy.ptr.p_double[0], 1, ae_v_len(0,n-1), v);
896
897 /*
898 * update YN/YNS
899 */
900 v = state->rkc.ptr.p_double[k];
901 ae_v_addd(&state->yn.ptr.p_double[0], 1, &state->rkk.ptr.pp_double[k][0], 1, ae_v_len(0,n-1), v);
902 v = state->rkcs.ptr.p_double[k];
903 ae_v_addd(&state->yns.ptr.p_double[0], 1, &state->rkk.ptr.pp_double[k][0], 1, ae_v_len(0,n-1), v);
904 k = k+1;
905 goto lbl_8;
906 lbl_10:
907
908 /*
909 * estimate error
910 */
911 err = (double)(0);
912 for(j=0; j<=n-1; j++)
913 {
914 if( !state->fraceps )
915 {
916
917 /*
918 * absolute error is estimated
919 */
920 err = ae_maxreal(err, ae_fabs(state->yn.ptr.p_double[j]-state->yns.ptr.p_double[j], _state), _state);
921 }
922 else
923 {
924
925 /*
926 * Relative error is estimated
927 */
928 v = state->escale.ptr.p_double[j];
929 if( ae_fp_eq(v,(double)(0)) )
930 {
931 v = (double)(1);
932 }
933 err = ae_maxreal(err, ae_fabs(state->yn.ptr.p_double[j]-state->yns.ptr.p_double[j], _state)/v, _state);
934 }
935 }
936
937 /*
938 * calculate new step, restart if necessary
939 */
940 if( ae_fp_less_eq(maxgrowpow*err,state->eps) )
941 {
942 h2 = odesolver_odesolvermaxgrow*h;
943 }
944 else
945 {
946 h2 = h*ae_pow(state->eps/err, 0.2, _state);
947 }
948 if( ae_fp_less(h2,h/odesolver_odesolvermaxshrink) )
949 {
950 h2 = h/odesolver_odesolvermaxshrink;
951 }
952 if( ae_fp_greater(err,state->eps) )
953 {
954 h = h2;
955 goto lbl_6;
956 }
957
958 /*
959 * advance position
960 */
961 xc = xc+h;
962 ae_v_move(&state->yc.ptr.p_double[0], 1, &state->yn.ptr.p_double[0], 1, ae_v_len(0,n-1));
963
964 /*
965 * update H
966 */
967 h = h2;
968
969 /*
970 * break on grid point
971 */
972 if( gridpoint )
973 {
974 goto lbl_7;
975 }
976 goto lbl_6;
977 lbl_7:
978
979 /*
980 * save result
981 */
982 ae_v_move(&state->ytbl.ptr.pp_double[i][0], 1, &state->yc.ptr.p_double[0], 1, ae_v_len(0,n-1));
983 i = i+1;
984 goto lbl_3;
985 lbl_5:
986 state->repterminationtype = 1;
987 result = ae_false;
988 return result;
989 lbl_1:
990 result = ae_false;
991 return result;
992
993 /*
994 * Saving state
995 */
996 lbl_rcomm:
997 result = ae_true;
998 state->rstate.ia.ptr.p_int[0] = n;
999 state->rstate.ia.ptr.p_int[1] = m;
1000 state->rstate.ia.ptr.p_int[2] = i;
1001 state->rstate.ia.ptr.p_int[3] = j;
1002 state->rstate.ia.ptr.p_int[4] = k;
1003 state->rstate.ia.ptr.p_int[5] = klimit;
1004 state->rstate.ba.ptr.p_bool[0] = gridpoint;
1005 state->rstate.ra.ptr.p_double[0] = xc;
1006 state->rstate.ra.ptr.p_double[1] = v;
1007 state->rstate.ra.ptr.p_double[2] = h;
1008 state->rstate.ra.ptr.p_double[3] = h2;
1009 state->rstate.ra.ptr.p_double[4] = err;
1010 state->rstate.ra.ptr.p_double[5] = maxgrowpow;
1011 return result;
1012 }
1013
1014
1015 /*************************************************************************
1016 ODE solver results
1017
1018 Called after OdeSolverIteration returned False.
1019
1020 INPUT PARAMETERS:
1021 State - algorithm state (used by OdeSolverIteration).
1022
1023 OUTPUT PARAMETERS:
1024 M - number of tabulated values, M>=1
1025 XTbl - array[0..M-1], values of X
1026 YTbl - array[0..M-1,0..N-1], values of Y in X[i]
1027 Rep - solver report:
1028 * Rep.TerminationType completetion code:
1029 * -2 X is not ordered by ascending/descending or
1030 there are non-distinct X[], i.e. X[i]=X[i+1]
1031 * -1 incorrect parameters were specified
1032 * 1 task has been solved
1033 * Rep.NFEV contains number of function calculations
1034
1035 -- ALGLIB --
1036 Copyright 01.09.2009 by Bochkanov Sergey
1037 *************************************************************************/
odesolverresults(odesolverstate * state,ae_int_t * m,ae_vector * xtbl,ae_matrix * ytbl,odesolverreport * rep,ae_state * _state)1038 void odesolverresults(odesolverstate* state,
1039 ae_int_t* m,
1040 /* Real */ ae_vector* xtbl,
1041 /* Real */ ae_matrix* ytbl,
1042 odesolverreport* rep,
1043 ae_state *_state)
1044 {
1045 double v;
1046 ae_int_t i;
1047
1048 *m = 0;
1049 ae_vector_clear(xtbl);
1050 ae_matrix_clear(ytbl);
1051 _odesolverreport_clear(rep);
1052
1053 rep->terminationtype = state->repterminationtype;
1054 if( rep->terminationtype>0 )
1055 {
1056 *m = state->m;
1057 rep->nfev = state->repnfev;
1058 ae_vector_set_length(xtbl, state->m, _state);
1059 v = state->xscale;
1060 ae_v_moved(&xtbl->ptr.p_double[0], 1, &state->xg.ptr.p_double[0], 1, ae_v_len(0,state->m-1), v);
1061 ae_matrix_set_length(ytbl, state->m, state->n, _state);
1062 for(i=0; i<=state->m-1; i++)
1063 {
1064 ae_v_move(&ytbl->ptr.pp_double[i][0], 1, &state->ytbl.ptr.pp_double[i][0], 1, ae_v_len(0,state->n-1));
1065 }
1066 }
1067 else
1068 {
1069 rep->nfev = 0;
1070 }
1071 }
1072
1073
1074 /*************************************************************************
1075 Internal initialization subroutine
1076 *************************************************************************/
odesolver_odesolverinit(ae_int_t solvertype,ae_vector * y,ae_int_t n,ae_vector * x,ae_int_t m,double eps,double h,odesolverstate * state,ae_state * _state)1077 static void odesolver_odesolverinit(ae_int_t solvertype,
1078 /* Real */ ae_vector* y,
1079 ae_int_t n,
1080 /* Real */ ae_vector* x,
1081 ae_int_t m,
1082 double eps,
1083 double h,
1084 odesolverstate* state,
1085 ae_state *_state)
1086 {
1087 ae_int_t i;
1088 double v;
1089
1090 _odesolverstate_clear(state);
1091
1092
1093 /*
1094 * Prepare RComm
1095 */
1096 ae_vector_set_length(&state->rstate.ia, 5+1, _state);
1097 ae_vector_set_length(&state->rstate.ba, 0+1, _state);
1098 ae_vector_set_length(&state->rstate.ra, 5+1, _state);
1099 state->rstate.stage = -1;
1100 state->needdy = ae_false;
1101
1102 /*
1103 * check parameters.
1104 */
1105 if( (n<=0||m<1)||ae_fp_eq(eps,(double)(0)) )
1106 {
1107 state->repterminationtype = -1;
1108 return;
1109 }
1110 if( ae_fp_less(h,(double)(0)) )
1111 {
1112 h = -h;
1113 }
1114
1115 /*
1116 * quick exit if necessary.
1117 * after this block we assume that M>1
1118 */
1119 if( m==1 )
1120 {
1121 state->repnfev = 0;
1122 state->repterminationtype = 1;
1123 ae_matrix_set_length(&state->ytbl, 1, n, _state);
1124 ae_v_move(&state->ytbl.ptr.pp_double[0][0], 1, &y->ptr.p_double[0], 1, ae_v_len(0,n-1));
1125 ae_vector_set_length(&state->xg, m, _state);
1126 ae_v_move(&state->xg.ptr.p_double[0], 1, &x->ptr.p_double[0], 1, ae_v_len(0,m-1));
1127 return;
1128 }
1129
1130 /*
1131 * check again: correct order of X[]
1132 */
1133 if( ae_fp_eq(x->ptr.p_double[1],x->ptr.p_double[0]) )
1134 {
1135 state->repterminationtype = -2;
1136 return;
1137 }
1138 for(i=1; i<=m-1; i++)
1139 {
1140 if( (ae_fp_greater(x->ptr.p_double[1],x->ptr.p_double[0])&&ae_fp_less_eq(x->ptr.p_double[i],x->ptr.p_double[i-1]))||(ae_fp_less(x->ptr.p_double[1],x->ptr.p_double[0])&&ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i-1])) )
1141 {
1142 state->repterminationtype = -2;
1143 return;
1144 }
1145 }
1146
1147 /*
1148 * auto-select H if necessary
1149 */
1150 if( ae_fp_eq(h,(double)(0)) )
1151 {
1152 v = ae_fabs(x->ptr.p_double[1]-x->ptr.p_double[0], _state);
1153 for(i=2; i<=m-1; i++)
1154 {
1155 v = ae_minreal(v, ae_fabs(x->ptr.p_double[i]-x->ptr.p_double[i-1], _state), _state);
1156 }
1157 h = 0.001*v;
1158 }
1159
1160 /*
1161 * store parameters
1162 */
1163 state->n = n;
1164 state->m = m;
1165 state->h = h;
1166 state->eps = ae_fabs(eps, _state);
1167 state->fraceps = ae_fp_less(eps,(double)(0));
1168 ae_vector_set_length(&state->xg, m, _state);
1169 ae_v_move(&state->xg.ptr.p_double[0], 1, &x->ptr.p_double[0], 1, ae_v_len(0,m-1));
1170 if( ae_fp_greater(x->ptr.p_double[1],x->ptr.p_double[0]) )
1171 {
1172 state->xscale = (double)(1);
1173 }
1174 else
1175 {
1176 state->xscale = (double)(-1);
1177 ae_v_muld(&state->xg.ptr.p_double[0], 1, ae_v_len(0,m-1), -1);
1178 }
1179 ae_vector_set_length(&state->yc, n, _state);
1180 ae_v_move(&state->yc.ptr.p_double[0], 1, &y->ptr.p_double[0], 1, ae_v_len(0,n-1));
1181 state->solvertype = solvertype;
1182 state->repterminationtype = 0;
1183
1184 /*
1185 * Allocate arrays
1186 */
1187 ae_vector_set_length(&state->y, n, _state);
1188 ae_vector_set_length(&state->dy, n, _state);
1189 }
1190
1191
_odesolverstate_init(void * _p,ae_state * _state,ae_bool make_automatic)1192 void _odesolverstate_init(void* _p, ae_state *_state, ae_bool make_automatic)
1193 {
1194 odesolverstate *p = (odesolverstate*)_p;
1195 ae_touch_ptr((void*)p);
1196 ae_vector_init(&p->yc, 0, DT_REAL, _state, make_automatic);
1197 ae_vector_init(&p->escale, 0, DT_REAL, _state, make_automatic);
1198 ae_vector_init(&p->xg, 0, DT_REAL, _state, make_automatic);
1199 ae_vector_init(&p->y, 0, DT_REAL, _state, make_automatic);
1200 ae_vector_init(&p->dy, 0, DT_REAL, _state, make_automatic);
1201 ae_matrix_init(&p->ytbl, 0, 0, DT_REAL, _state, make_automatic);
1202 ae_vector_init(&p->yn, 0, DT_REAL, _state, make_automatic);
1203 ae_vector_init(&p->yns, 0, DT_REAL, _state, make_automatic);
1204 ae_vector_init(&p->rka, 0, DT_REAL, _state, make_automatic);
1205 ae_vector_init(&p->rkc, 0, DT_REAL, _state, make_automatic);
1206 ae_vector_init(&p->rkcs, 0, DT_REAL, _state, make_automatic);
1207 ae_matrix_init(&p->rkb, 0, 0, DT_REAL, _state, make_automatic);
1208 ae_matrix_init(&p->rkk, 0, 0, DT_REAL, _state, make_automatic);
1209 _rcommstate_init(&p->rstate, _state, make_automatic);
1210 }
1211
1212
_odesolverstate_init_copy(void * _dst,void * _src,ae_state * _state,ae_bool make_automatic)1213 void _odesolverstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic)
1214 {
1215 odesolverstate *dst = (odesolverstate*)_dst;
1216 odesolverstate *src = (odesolverstate*)_src;
1217 dst->n = src->n;
1218 dst->m = src->m;
1219 dst->xscale = src->xscale;
1220 dst->h = src->h;
1221 dst->eps = src->eps;
1222 dst->fraceps = src->fraceps;
1223 ae_vector_init_copy(&dst->yc, &src->yc, _state, make_automatic);
1224 ae_vector_init_copy(&dst->escale, &src->escale, _state, make_automatic);
1225 ae_vector_init_copy(&dst->xg, &src->xg, _state, make_automatic);
1226 dst->solvertype = src->solvertype;
1227 dst->needdy = src->needdy;
1228 dst->x = src->x;
1229 ae_vector_init_copy(&dst->y, &src->y, _state, make_automatic);
1230 ae_vector_init_copy(&dst->dy, &src->dy, _state, make_automatic);
1231 ae_matrix_init_copy(&dst->ytbl, &src->ytbl, _state, make_automatic);
1232 dst->repterminationtype = src->repterminationtype;
1233 dst->repnfev = src->repnfev;
1234 ae_vector_init_copy(&dst->yn, &src->yn, _state, make_automatic);
1235 ae_vector_init_copy(&dst->yns, &src->yns, _state, make_automatic);
1236 ae_vector_init_copy(&dst->rka, &src->rka, _state, make_automatic);
1237 ae_vector_init_copy(&dst->rkc, &src->rkc, _state, make_automatic);
1238 ae_vector_init_copy(&dst->rkcs, &src->rkcs, _state, make_automatic);
1239 ae_matrix_init_copy(&dst->rkb, &src->rkb, _state, make_automatic);
1240 ae_matrix_init_copy(&dst->rkk, &src->rkk, _state, make_automatic);
1241 _rcommstate_init_copy(&dst->rstate, &src->rstate, _state, make_automatic);
1242 }
1243
1244
_odesolverstate_clear(void * _p)1245 void _odesolverstate_clear(void* _p)
1246 {
1247 odesolverstate *p = (odesolverstate*)_p;
1248 ae_touch_ptr((void*)p);
1249 ae_vector_clear(&p->yc);
1250 ae_vector_clear(&p->escale);
1251 ae_vector_clear(&p->xg);
1252 ae_vector_clear(&p->y);
1253 ae_vector_clear(&p->dy);
1254 ae_matrix_clear(&p->ytbl);
1255 ae_vector_clear(&p->yn);
1256 ae_vector_clear(&p->yns);
1257 ae_vector_clear(&p->rka);
1258 ae_vector_clear(&p->rkc);
1259 ae_vector_clear(&p->rkcs);
1260 ae_matrix_clear(&p->rkb);
1261 ae_matrix_clear(&p->rkk);
1262 _rcommstate_clear(&p->rstate);
1263 }
1264
1265
_odesolverstate_destroy(void * _p)1266 void _odesolverstate_destroy(void* _p)
1267 {
1268 odesolverstate *p = (odesolverstate*)_p;
1269 ae_touch_ptr((void*)p);
1270 ae_vector_destroy(&p->yc);
1271 ae_vector_destroy(&p->escale);
1272 ae_vector_destroy(&p->xg);
1273 ae_vector_destroy(&p->y);
1274 ae_vector_destroy(&p->dy);
1275 ae_matrix_destroy(&p->ytbl);
1276 ae_vector_destroy(&p->yn);
1277 ae_vector_destroy(&p->yns);
1278 ae_vector_destroy(&p->rka);
1279 ae_vector_destroy(&p->rkc);
1280 ae_vector_destroy(&p->rkcs);
1281 ae_matrix_destroy(&p->rkb);
1282 ae_matrix_destroy(&p->rkk);
1283 _rcommstate_destroy(&p->rstate);
1284 }
1285
1286
_odesolverreport_init(void * _p,ae_state * _state,ae_bool make_automatic)1287 void _odesolverreport_init(void* _p, ae_state *_state, ae_bool make_automatic)
1288 {
1289 odesolverreport *p = (odesolverreport*)_p;
1290 ae_touch_ptr((void*)p);
1291 }
1292
1293
_odesolverreport_init_copy(void * _dst,void * _src,ae_state * _state,ae_bool make_automatic)1294 void _odesolverreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic)
1295 {
1296 odesolverreport *dst = (odesolverreport*)_dst;
1297 odesolverreport *src = (odesolverreport*)_src;
1298 dst->nfev = src->nfev;
1299 dst->terminationtype = src->terminationtype;
1300 }
1301
1302
_odesolverreport_clear(void * _p)1303 void _odesolverreport_clear(void* _p)
1304 {
1305 odesolverreport *p = (odesolverreport*)_p;
1306 ae_touch_ptr((void*)p);
1307 }
1308
1309
_odesolverreport_destroy(void * _p)1310 void _odesolverreport_destroy(void* _p)
1311 {
1312 odesolverreport *p = (odesolverreport*)_p;
1313 ae_touch_ptr((void*)p);
1314 }
1315
1316
1317 #endif
1318
1319 }
1320
1321