Mercurial > hg > octave-lyh

comparison liboctave/DASSL.cc @ 256:e592734b002b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-12-08 23:36:33 by jwe]
author jwe
date Wed, 08 Dec 1993 23:36:43 +0000
parents 780cbbc57b7c
children c23f50e61c58
comparison
equal deleted inserted replaced
255:98246fedc941 256:e592734b002b
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 (ddassl) (int (*)(), const int*, double*, double*, 34 int F77_FCN (ddassl) (int (*)(double*, double*, double*, double*,
35 double*, double*, const int*, const double*, 35 int*, double*, int*),
36 const int*, double*, double*, double*,
37 double*, const int*, const double*,
36 const double*, int*, double*, const int*, 38 const double*, int*, double*, const int*,
37 int*, const int*, const double*, const int*, 39 int*, const int*, const double*, const int*,
38 int (*)()); 40 int (*)(double*, double*, double*, double*,
41 double*, double*, int*));
39 } 42 }
40 43
41 static DAERHSFunc user_fun; 44 static DAERHSFunc user_fun;
42 static DAEJacFunc user_jac; 45 static DAEJacFunc user_jac;
43 static int nn; 46 static int nn;
48 t = 0.0; 51 t = 0.0;
49 52
50 stop_time_set = 0; 53 stop_time_set = 0;
51 stop_time = 0.0; 54 stop_time = 0.0;
52 55
56 integration_error = 0;
53 restart = 1; 57 restart = 1;
54 58
55 DAEFunc::set_function (NULL); 59 DAEFunc::set_function (NULL);
56 DAEFunc::set_jacobian_function (NULL); 60 DAEFunc::set_jacobian_function (NULL);
57 61
75 relative_tolerance = 1.0e-6; 79 relative_tolerance = 1.0e-6;
76 80
77 stop_time_set = 0; 81 stop_time_set = 0;
78 stop_time = 0.0; 82 stop_time = 0.0;
79 83
84 integration_error = 0;
80 restart = 1; 85 restart = 1;
81 86
82 DAEFunc::set_function (NULL); 87 DAEFunc::set_function (NULL);
83 DAEFunc::set_jacobian_function (NULL); 88 DAEFunc::set_jacobian_function (NULL);
84 89
104 relative_tolerance = 1.0e-6; 109 relative_tolerance = 1.0e-6;
105 110
106 stop_time_set = 0; 111 stop_time_set = 0;
107 stop_time = 0.0; 112 stop_time = 0.0;
108 113
114 integration_error = 0;
109 restart = 1; 115 restart = 1;
110 116
111 DAEFunc::set_function (f.function ()); 117 DAEFunc::set_function (f.function ());
112 DAEFunc::set_jacobian_function (f.jacobian_function ()); 118 DAEFunc::set_jacobian_function (f.jacobian_function ());
113 119
172 } 178 }
173 179
174 void 180 void
175 DAE::initialize (Vector& state, double time) 181 DAE::initialize (Vector& state, double time)
176 { 182 {
183 integration_error = 0;
177 restart = 1; 184 restart = 1;
178 x = state; 185 x = state;
179 int nx = x.capacity (); 186 int nx = x.capacity ();
180 xdot.resize (nx, 0.0); 187 xdot.resize (nx, 0.0);
181 t = time; 188 t = time;
182 } 189 }
183 190
184 void 191 void
185 DAE::initialize (Vector& state, Vector& deriv, double time) 192 DAE::initialize (Vector& state, Vector& deriv, double time)
186 { 193 {
194 integration_error = 0;
187 restart = 1; 195 restart = 1;
188 xdot = deriv; 196 xdot = deriv;
189 x = state; 197 x = state;
190 t = time; 198 t = time;
191 } 199 }
204 tmp_state.elem (i) = state [i]; 212 tmp_state.elem (i) = state [i];
205 } 213 }
206 214
207 tmp_delta = user_fun (tmp_state, tmp_deriv, *time); 215 tmp_delta = user_fun (tmp_state, tmp_deriv, *time);
208 216
209 for (i = 0; i < nn; i++) 217 if (tmp_delta.length () == 0)
210 delta [i] = tmp_delta.elem (i); 218 *ires = -2;
219 else
220 {
221 for (i = 0; i < nn; i++)
222 delta [i] = tmp_delta.elem (i);
223 }
211 224
212 return 0; 225 return 0;
213 } 226 }
214 227
215 int 228 int
242 } 255 }
243 256
244 Vector 257 Vector
245 DAE::integrate (double tout) 258 DAE::integrate (double tout)
246 { 259 {
260 integration_error = 0;
261
247 if (DAEFunc::jac == NULL) 262 if (DAEFunc::jac == NULL)
248 iwork [4] = 0; 263 iwork [4] = 0;
249 else 264 else
250 iwork [4] = 1; 265 iwork [4] = 1;
251 266
292 case 3: // The integration to TOUT was successfully completed 307 case 3: // The integration to TOUT was successfully completed
293 // (T=TOUT) by stepping past TOUT. Y(*) is obtained by 308 // (T=TOUT) by stepping past TOUT. Y(*) is obtained by
294 // interpolation. YPRIME(*) is obtained by interpolation. 309 // interpolation. YPRIME(*) is obtained by interpolation.
295 break; 310 break;
296 case -1: // A large amount of work has been expended. (About 500 steps). 311 case -1: // A large amount of work has been expended. (About 500 steps).
297 break;
298 case -2: // The error tolerances are too stringent. 312 case -2: // The error tolerances are too stringent.
299 break;
300 case -3: // The local error test cannot be satisfied because you 313 case -3: // The local error test cannot be satisfied because you
301 // specified a zero component in ATOL and the 314 // specified a zero component in ATOL and the
302 // corresponding computed solution component is zero. 315 // corresponding computed solution component is zero.
303 // Thus, a pure relative error test is impossible for 316 // Thus, a pure relative error test is impossible for
304 // this component. 317 // this component.
305 break;
306 case -6: // DDASSL had repeated error test failures on the last 318 case -6: // DDASSL had repeated error test failures on the last
307 // attempted step. 319 // attempted step.
308 break;
309 case -7: // The corrector could not converge. 320 case -7: // The corrector could not converge.
310 break;
311 case -8: // The matrix of partial derivatives is singular. 321 case -8: // The matrix of partial derivatives is singular.
312 break;
313 case -9: // The corrector could not converge. There were repeated 322 case -9: // The corrector could not converge. There were repeated
314 // error test failures in this step. 323 // error test failures in this step.
315 break;
316 case -10: // The corrector could not converge because IRES was 324 case -10: // The corrector could not converge because IRES was
317 // equal to minus one. 325 // equal to minus one.
318 break;
319 case -11: // IRES equal to -2 was encountered and control is being 326 case -11: // IRES equal to -2 was encountered and control is being
320 // returned to the calling program. 327 // returned to the calling program.
321 break;
322 case -12: // DDASSL failed to compute the initial YPRIME. 328 case -12: // DDASSL failed to compute the initial YPRIME.
323 break;
324 case -33: // The code has encountered trouble from which it cannot 329 case -33: // The code has encountered trouble from which it cannot
325 // recover. A message is printed explaining the trouble 330 // recover. A message is printed explaining the trouble
326 // and control is returned to the calling program. For 331 // and control is returned to the calling program. For
327 // example, this occurs when invalid input is detected. 332 // example, this occurs when invalid input is detected.
328 break;
329 default: 333 default:
330 // Error? 334 integration_error = 1;
331 break; 335 break;
332 } 336 }
333 337
334 t = tout; 338 t = tout;
335 339
354 } 358 }
355 359
356 for (int j = 1; j < n_out; j++) 360 for (int j = 1; j < n_out; j++)
357 { 361 {
358 ColumnVector x_next = integrate (tout.elem (j)); 362 ColumnVector x_next = integrate (tout.elem (j));
363
364 if (integration_error)
365 return retval;
366
359 for (i = 0; i < n; i++) 367 for (i = 0; i < n; i++)
360 { 368 {
361 retval.elem (j, i) = x_next.elem (i); 369 retval.elem (j, i) = x_next.elem (i);
362 xdot_out.elem (j, i) = xdot.elem (i); 370 xdot_out.elem (j, i) = xdot.elem (i);
363 } 371 }
438 i_out++; 446 i_out++;
439 } 447 }
440 448
441 ColumnVector x_next = integrate (t_out); 449 ColumnVector x_next = integrate (t_out);
442 450
451 if (integration_error)
452 return retval;
453
443 if (save_output) 454 if (save_output)
444 { 455 {
445 for (i = 0; i < n; i++) 456 for (i = 0; i < n; i++)
446 { 457 {
447 retval.elem (i_out-1, i) = x_next.elem (i); 458 retval.elem (i_out-1, i) = x_next.elem (i);
452 if (do_restart) 463 if (do_restart)
453 force_restart (); 464 force_restart ();
454 } 465 }
455 } 466 }
456 else 467 else
457 retval = integrate (tout, xdot_out); 468 {
469 retval = integrate (tout, xdot_out);
470
471 if (integration_error)
472 return retval;
473 }
458 } 474 }
459 475
460 return retval; 476 return retval;
461 } 477 }