Mercurial > octave

comparison liboctave/LSODE.cc @ 258:1c9a678906fb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-12-14 06:24:04 by jwe]
author jwe
date Tue, 14 Dec 1993 06:24:21 +0000
parents 780cbbc57b7c
children c23f50e61c58
comparison
equal deleted inserted replaced
257:126791334c68 258:1c9a678906fb
29 #include "f77-uscore.h" 29 #include "f77-uscore.h"
30 #include "lo-error.h" 30 #include "lo-error.h"
31 31
32 extern "C" 32 extern "C"
33 { 33 {
34 int F77_FCN (lsode) (int (*)(), int *, double *, double *, double *, 34 int F77_FCN (lsode) (int (*)(int*, double*, double*, double*, int*),
35 int *, double *, double *, double *,
35 int *, double *, double *, int *, int *, int *, 36 int *, double *, double *, int *, int *, int *,
36 double *, int *, int *, int *, int (*)(), int *); 37 double *, int *, int *, int *,
38 int (*)(int*, double*, double*, int*, int*,
39 double*, int*), int *);
37 } 40 }
38 41
39 static ColumnVector (*user_fun) (ColumnVector&, double); 42 static ColumnVector (*user_fun) (ColumnVector&, double);
40 static Matrix (*user_jac) (ColumnVector&, double); 43 static Matrix (*user_jac) (ColumnVector&, double);
41 static ColumnVector *tmp_x; 44 static ColumnVector *tmp_x;
49 relative_tolerance = 1.0e-6; 52 relative_tolerance = 1.0e-6;
50 53
51 stop_time_set = 0; 54 stop_time_set = 0;
52 stop_time = 0.0; 55 stop_time = 0.0;
53 56
57 integration_error = 0;
54 restart = 1; 58 restart = 1;
55 59
56 istate = 1; 60 istate = 1;
57 itol = 1; 61 itol = 1;
58 itask = 1; 62 itask = 1;
82 relative_tolerance = 1.0e-6; 86 relative_tolerance = 1.0e-6;
83 87
84 stop_time_set = 0; 88 stop_time_set = 0;
85 stop_time = 0.0; 89 stop_time = 0.0;
86 90
91 integration_error = 0;
87 restart = 1; 92 restart = 1;
88 93
89 istate = 1; 94 istate = 1;
90 itol = 1; 95 itol = 1;
91 itask = 1; 96 itask = 1;
116 relative_tolerance = 1.0e-6; 121 relative_tolerance = 1.0e-6;
117 122
118 stop_time_set = 0; 123 stop_time_set = 0;
119 stop_time = 0.0; 124 stop_time = 0.0;
120 125
126 integration_error = 0;
121 restart = 1; 127 restart = 1;
122 128
123 istate = 1; 129 istate = 1;
124 itol = 1; 130 itol = 1;
125 itask = 1; 131 itask = 1;
145 delete [] rwork; 151 delete [] rwork;
146 delete [] iwork; 152 delete [] iwork;
147 } 153 }
148 154
149 int 155 int
150 lsode_f (int *neq, double *time, double *state, double *deriv) 156 lsode_f (int *neq, double *time, double *state, double *deriv, int *ierr)
151 { 157 {
152 int nn = *neq; 158 int nn = *neq;
153 ColumnVector tmp_deriv (nn); 159 ColumnVector tmp_deriv (nn);
154 160
155 /* 161 /*
157 * In that case we have to create a temporary vector object 163 * In that case we have to create a temporary vector object
158 * and copy. 164 * and copy.
159 */ 165 */
160 tmp_deriv = (*user_fun) (*tmp_x, *time); 166 tmp_deriv = (*user_fun) (*tmp_x, *time);
161 167
162 for (int i = 0; i < nn; i++) 168 if (tmp_deriv.length () == 0)
163 deriv [i] = tmp_deriv.elem (i); 169 *ierr = -1;
170 else
171 {
172 for (int i = 0; i < nn; i++)
173 deriv [i] = tmp_deriv.elem (i);
174 }
164 175
165 return 0; 176 return 0;
166 } 177 }
167 178
168 int 179 int
192 if (jac == NULL) 203 if (jac == NULL)
193 method_flag = 22; 204 method_flag = 22;
194 else 205 else
195 method_flag = 21; 206 method_flag = 21;
196 207
208 integration_error = 0;
209
197 double *xp = x.fortran_vec (); 210 double *xp = x.fortran_vec ();
198 211
199 // NOTE: this won't work if LSODE passes copies of the state vector. 212 // NOTE: this won't work if LSODE passes copies of the state vector.
200 // In that case we have to create a temporary vector object 213 // In that case we have to create a temporary vector object
201 // and copy. 214 // and copy.
234 &itask, &istate, &iopt, rwork, &lrw, iwork, 247 &itask, &istate, &iopt, rwork, &lrw, iwork,
235 &liw, lsode_j, &method_flag); 248 &liw, lsode_j, &method_flag);
236 249
237 switch (istate) 250 switch (istate)
238 { 251 {
252 case -13: // Return requested in user-supplied function.
239 case -6: // error weight became zero during problem. (solution 253 case -6: // error weight became zero during problem. (solution
240 // component i vanished, and atol or atol(i) = 0.) 254 // component i vanished, and atol or atol(i) = 0.)
241 break;
242 case -5: // repeated convergence failures (perhaps bad jacobian 255 case -5: // repeated convergence failures (perhaps bad jacobian
243 // supplied or wrong choice of mf or tolerances). 256 // supplied or wrong choice of mf or tolerances).
244 break;
245 case -4: // repeated error test failures (check all inputs). 257 case -4: // repeated error test failures (check all inputs).
246 break;
247 case -3: // illegal input detected (see printed message). 258 case -3: // illegal input detected (see printed message).
248 break;
249 case -2: // excess accuracy requested (tolerances too small). 259 case -2: // excess accuracy requested (tolerances too small).
260 integration_error = 1;
261 return ColumnVector ();
250 break; 262 break;
251 case -1: // excess work done on this call (perhaps wrong mf). 263 case -1: // excess work done on this call (perhaps wrong mf).
252 working_too_hard++; 264 working_too_hard++;
253 if (working_too_hard > 20) 265 if (working_too_hard > 20)
254 { 266 {
255 (*current_liboctave_error_handler) 267 (*current_liboctave_error_handler)
256 ("giving up after more than %d steps attempted in lsode", 268 ("giving up after more than %d steps attempted in lsode",
257 iwork[5] * 20); 269 iwork[5] * 20);
270 integration_error = 1;
258 return ColumnVector (); 271 return ColumnVector ();
259 } 272 }
260 else 273 else
261 { 274 {
262 istate = 2; 275 istate = 2;
315 retval.elem (0, i) = x.elem (i); 328 retval.elem (0, i) = x.elem (i);
316 329
317 for (int j = 1; j < n_out; j++) 330 for (int j = 1; j < n_out; j++)
318 { 331 {
319 ColumnVector x_next = integrate (tout.elem (j)); 332 ColumnVector x_next = integrate (tout.elem (j));
333
334 if (integration_error)
335 return retval;
336
320 for (i = 0; i < n; i++) 337 for (i = 0; i < n; i++)
321 retval.elem (j, i) = x_next.elem (i); 338 retval.elem (j, i) = x_next.elem (i);
322 } 339 }
323 } 340 }
324 341
392 i_out++; 409 i_out++;
393 } 410 }
394 411
395 ColumnVector x_next = integrate (t_out); 412 ColumnVector x_next = integrate (t_out);
396 413
414 if (integration_error)
415 return retval;
416
397 if (save_output) 417 if (save_output)
398 { 418 {
399 for (i = 0; i < n; i++) 419 for (i = 0; i < n; i++)
400 retval.elem (i_out-1, i) = x_next.elem (i); 420 retval.elem (i_out-1, i) = x_next.elem (i);
401 } 421 }
403 if (do_restart) 423 if (do_restart)
404 force_restart (); 424 force_restart ();
405 } 425 }
406 } 426 }
407 else 427 else
408 retval = integrate (tout); 428 {
429 retval = integrate (tout);
430
431 if (integration_error)
432 return retval;
433 }
409 } 434 }
410 435
411 return retval; 436 return retval;
412 } 437 }
413 438