Mercurial > octave

comparison liboctave/LSODE.cc @ 5275:23b37da9fd5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-08 16:07:35 by jwe]
author jwe
date Fri, 08 Apr 2005 16:07:37 +0000
parents e35b034d3523
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5274:eae7b40388e9 5275:23b37da9fd5b
31 #include "f77-fcn.h" 31 #include "f77-fcn.h"
32 #include "lo-error.h" 32 #include "lo-error.h"
33 #include "lo-sstream.h" 33 #include "lo-sstream.h"
34 #include "quit.h" 34 #include "quit.h"
35 35
36 typedef int (*lsode_fcn_ptr) (const int&, const double&, double*, 36 typedef octave_idx_type (*lsode_fcn_ptr) (const octave_idx_type&, const double&, double*,
37 double*, int&); 37 double*, octave_idx_type&);
38 38
39 typedef int (*lsode_jac_ptr) (const int&, const double&, double*, 39 typedef octave_idx_type (*lsode_jac_ptr) (const octave_idx_type&, const double&, double*,
40 const int&, const int&, double*, const 40 const octave_idx_type&, const octave_idx_type&, double*, const
41 int&); 41 octave_idx_type&);
42 42
43 extern "C" 43 extern "C"
44 { 44 {
45 F77_RET_T 45 F77_RET_T
46 F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, int&, double*, double&, 46 F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, octave_idx_type&, double*, double&,
47 double&, int&, double&, const double*, int&, 47 double&, octave_idx_type&, double&, const double*, octave_idx_type&,
48 int&, int&, double*, int&, int*, int&, 48 octave_idx_type&, octave_idx_type&, double*, octave_idx_type&, octave_idx_type*, octave_idx_type&,
49 lsode_jac_ptr, int&); 49 lsode_jac_ptr, octave_idx_type&);
50 } 50 }
51 51
52 static ODEFunc::ODERHSFunc user_fun; 52 static ODEFunc::ODERHSFunc user_fun;
53 static ODEFunc::ODEJacFunc user_jac; 53 static ODEFunc::ODEJacFunc user_jac;
54 static ColumnVector *tmp_x; 54 static ColumnVector *tmp_x;
55 55
56 static int 56 static octave_idx_type
57 lsode_f (const int& neq, const double& time, double *, 57 lsode_f (const octave_idx_type& neq, const double& time, double *,
58 double *deriv, int& ierr) 58 double *deriv, octave_idx_type& ierr)
59 { 59 {
60 BEGIN_INTERRUPT_WITH_EXCEPTIONS; 60 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
61 61
62 ColumnVector tmp_deriv; 62 ColumnVector tmp_deriv;
63 63
69 69
70 if (tmp_deriv.length () == 0) 70 if (tmp_deriv.length () == 0)
71 ierr = -1; 71 ierr = -1;
72 else 72 else
73 { 73 {
74 for (int i = 0; i < neq; i++) 74 for (octave_idx_type i = 0; i < neq; i++)
75 deriv [i] = tmp_deriv.elem (i); 75 deriv [i] = tmp_deriv.elem (i);
76 } 76 }
77 77
78 END_INTERRUPT_WITH_EXCEPTIONS; 78 END_INTERRUPT_WITH_EXCEPTIONS;
79 79
80 return 0; 80 return 0;
81 } 81 }
82 82
83 static int 83 static octave_idx_type
84 lsode_j (const int& neq, const double& time, double *, 84 lsode_j (const octave_idx_type& neq, const double& time, double *,
85 const int&, const int&, double *pd, const int& nrowpd) 85 const octave_idx_type&, const octave_idx_type&, double *pd, const octave_idx_type& nrowpd)
86 { 86 {
87 BEGIN_INTERRUPT_WITH_EXCEPTIONS; 87 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
88 88
89 Matrix tmp_jac (neq, neq); 89 Matrix tmp_jac (neq, neq);
90 90
92 // In that case we have to create a temporary vector object 92 // In that case we have to create a temporary vector object
93 // and copy. 93 // and copy.
94 94
95 tmp_jac = (*user_jac) (*tmp_x, time); 95 tmp_jac = (*user_jac) (*tmp_x, time);
96 96
97 for (int j = 0; j < neq; j++) 97 for (octave_idx_type j = 0; j < neq; j++)
98 for (int i = 0; i < neq; i++) 98 for (octave_idx_type i = 0; i < neq; i++)
99 pd [nrowpd * j + i] = tmp_jac (i, j); 99 pd [nrowpd * j + i] = tmp_jac (i, j);
100 100
101 END_INTERRUPT_WITH_EXCEPTIONS; 101 END_INTERRUPT_WITH_EXCEPTIONS;
102 102
103 return 0; 103 return 0;
106 ColumnVector 106 ColumnVector
107 LSODE::do_integrate (double tout) 107 LSODE::do_integrate (double tout)
108 { 108 {
109 ColumnVector retval; 109 ColumnVector retval;
110 110
111 static int nn = 0; 111 static octave_idx_type nn = 0;
112 112
113 if (! initialized || restart || ODEFunc::reset || LSODE_options::reset) 113 if (! initialized || restart || ODEFunc::reset || LSODE_options::reset)
114 { 114 {
115 integration_error = false; 115 integration_error = false;
116 116
117 initialized = true; 117 initialized = true;
118 118
119 istate = 1; 119 istate = 1;
120 120
121 int n = size (); 121 octave_idx_type n = size ();
122 122
123 nn = n; 123 nn = n;
124 124
125 int max_maxord = 0; 125 octave_idx_type max_maxord = 0;
126 126
127 if (integration_method () == "stiff") 127 if (integration_method () == "stiff")
128 { 128 {
129 max_maxord = 5; 129 max_maxord = 5;
130 130
164 } 164 }
165 } 165 }
166 166
167 iwork.resize (liw); 167 iwork.resize (liw);
168 168
169 for (int i = 4; i < 9; i++) 169 for (octave_idx_type i = 4; i < 9; i++)
170 iwork(i) = 0; 170 iwork(i) = 0;
171 171
172 rwork.resize (lrw); 172 rwork.resize (lrw);
173 173
174 for (int i = 4; i < 9; i++) 174 for (octave_idx_type i = 4; i < 9; i++)
175 rwork(i) = 0; 175 rwork(i) = 0;
176 176
177 if (stop_time_set) 177 if (stop_time_set)
178 { 178 {
179 itask = 4; 179 itask = 4;
219 // LSODE_options 219 // LSODE_options
220 220
221 rel_tol = relative_tolerance (); 221 rel_tol = relative_tolerance ();
222 abs_tol = absolute_tolerance (); 222 abs_tol = absolute_tolerance ();
223 223
224 int abs_tol_len = abs_tol.length (); 224 octave_idx_type abs_tol_len = abs_tol.length ();
225 225
226 if (abs_tol_len == 1) 226 if (abs_tol_len == 1)
227 itol = 1; 227 itol = 1;
228 else if (abs_tol_len == n) 228 else if (abs_tol_len == n)
229 itol = 2; 229 itol = 2;
255 { 255 {
256 rwork(6) = minss; 256 rwork(6) = minss;
257 iopt = 1; 257 iopt = 1;
258 } 258 }
259 259
260 int sl = step_limit (); 260 octave_idx_type sl = step_limit ();
261 if (sl > 0) 261 if (sl > 0)
262 { 262 {
263 iwork(5) = sl; 263 iwork(5) = sl;
264 iopt = 1; 264 iopt = 1;
265 } 265 }
382 Matrix 382 Matrix
383 LSODE::do_integrate (const ColumnVector& tout) 383 LSODE::do_integrate (const ColumnVector& tout)
384 { 384 {
385 Matrix retval; 385 Matrix retval;
386 386
387 int n_out = tout.capacity (); 387 octave_idx_type n_out = tout.capacity ();
388 int n = size (); 388 octave_idx_type n = size ();
389 389
390 if (n_out > 0 && n > 0) 390 if (n_out > 0 && n > 0)
391 { 391 {
392 retval.resize (n_out, n); 392 retval.resize (n_out, n);
393 393
394 for (int i = 0; i < n; i++) 394 for (octave_idx_type i = 0; i < n; i++)
395 retval.elem (0, i) = x.elem (i); 395 retval.elem (0, i) = x.elem (i);
396 396
397 for (int j = 1; j < n_out; j++) 397 for (octave_idx_type j = 1; j < n_out; j++)
398 { 398 {
399 ColumnVector x_next = do_integrate (tout.elem (j)); 399 ColumnVector x_next = do_integrate (tout.elem (j));
400 400
401 if (integration_error) 401 if (integration_error)
402 return retval; 402 return retval;
403 403
404 for (int i = 0; i < n; i++) 404 for (octave_idx_type i = 0; i < n; i++)
405 retval.elem (j, i) = x_next.elem (i); 405 retval.elem (j, i) = x_next.elem (i);
406 } 406 }
407 } 407 }
408 408
409 return retval; 409 return retval;
412 Matrix 412 Matrix
413 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit) 413 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
414 { 414 {
415 Matrix retval; 415 Matrix retval;
416 416
417 int n_out = tout.capacity (); 417 octave_idx_type n_out = tout.capacity ();
418 int n = size (); 418 octave_idx_type n = size ();
419 419
420 if (n_out > 0 && n > 0) 420 if (n_out > 0 && n > 0)
421 { 421 {
422 retval.resize (n_out, n); 422 retval.resize (n_out, n);
423 423
424 for (int i = 0; i < n; i++) 424 for (octave_idx_type i = 0; i < n; i++)
425 retval.elem (0, i) = x.elem (i); 425 retval.elem (0, i) = x.elem (i);
426 426
427 int n_crit = tcrit.capacity (); 427 octave_idx_type n_crit = tcrit.capacity ();
428 428
429 if (n_crit > 0) 429 if (n_crit > 0)
430 { 430 {
431 int i_crit = 0; 431 octave_idx_type i_crit = 0;
432 int i_out = 1; 432 octave_idx_type i_out = 1;
433 double next_crit = tcrit.elem (0); 433 double next_crit = tcrit.elem (0);
434 double next_out; 434 double next_out;
435 while (i_out < n_out) 435 while (i_out < n_out)
436 { 436 {
437 bool do_restart = false; 437 bool do_restart = false;
438 438
439 next_out = tout.elem (i_out); 439 next_out = tout.elem (i_out);
440 if (i_crit < n_crit) 440 if (i_crit < n_crit)
441 next_crit = tcrit.elem (i_crit); 441 next_crit = tcrit.elem (i_crit);
442 442
443 int save_output; 443 octave_idx_type save_output;
444 double t_out; 444 double t_out;
445 445
446 if (next_crit == next_out) 446 if (next_crit == next_out)
447 { 447 {
448 set_stop_time (next_crit); 448 set_stop_time (next_crit);
483 if (integration_error) 483 if (integration_error)
484 return retval; 484 return retval;
485 485
486 if (save_output) 486 if (save_output)
487 { 487 {
488 for (int i = 0; i < n; i++) 488 for (octave_idx_type i = 0; i < n; i++)
489 retval.elem (i_out-1, i) = x_next.elem (i); 489 retval.elem (i_out-1, i) = x_next.elem (i);
490 } 490 }
491 491
492 if (do_restart) 492 if (do_restart)
493 force_restart (); 493 force_restart ();