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