Mercurial > octave-nkf

comparison liboctave/DASRT.cc @ 3994:a41827ec5677

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-07-16 23:57:09 by jwe]
author jwe
date Tue, 16 Jul 2002 23:57:09 +0000
parents f23bc69132cc
children ee0304212be0
comparison
equal deleted inserted replaced
3993:f23bc69132cc 3994:a41827ec5677
86 86
87 ColumnVector tmp_deriv (nn); 87 ColumnVector tmp_deriv (nn);
88 for (int i = 0; i < nn; i++) 88 for (int i = 0; i < nn; i++)
89 tmp_deriv(i) = deriv[i]; 89 tmp_deriv(i) = deriv[i];
90 90
91 ColumnVector tmp_fval = user_fsub (tmp_state, tmp_deriv, t, ires); 91 ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires);
92 92
93 if (tmp_fval.length () == 0) 93 if (tmp_fval.length () == 0)
94 ires = -2; 94 ires = -2;
95 else 95 else
96 { 96 {
99 } 99 }
100 100
101 return 0; 101 return 0;
102 } 102 }
103 103
104 //typedef int (*efptr) (const double& t, const int& n, double *state,
105 // double *ework, double *rpar, int *ipar,
106 // const int& ieform, int& ires);
107
108 //static efptr e_fun;
109
110 int 104 int
111 ddasrt_j (const double& time, const double *state, const double *deriv, 105 ddasrt_j (const double& time, const double *state, const double *deriv,
112 double *pd, const double& cj, double *, int *) 106 double *pd, const double& cj, double *, int *)
113 { 107 {
114 // XXX FIXME XXX -- would be nice to avoid copying the data. 108 // XXX FIXME XXX -- would be nice to avoid copying the data.
120 { 114 {
121 tmp_deriv.elem (i) = deriv [i]; 115 tmp_deriv.elem (i) = deriv [i];
122 tmp_state.elem (i) = state [i]; 116 tmp_state.elem (i) = state [i];
123 } 117 }
124 118
125 Matrix tmp_pd = user_jsub (tmp_state, tmp_deriv, time, cj); 119 Matrix tmp_pd = (*user_jsub) (tmp_state, tmp_deriv, time, cj);
126 120
127 for (int j = 0; j < nn; j++) 121 for (int j = 0; j < nn; j++)
128 for (int i = 0; i < nn; i++) 122 for (int i = 0; i < nn; i++)
129 pd [nn * j + i] = tmp_pd.elem (i, j); 123 pd [nn * j + i] = tmp_pd.elem (i, j);
130 124
139 133
140 ColumnVector tmp_state (n); 134 ColumnVector tmp_state (n);
141 for (int i = 0; i < n; i++) 135 for (int i = 0; i < n; i++)
142 tmp_state(i) = state[i]; 136 tmp_state(i) = state[i];
143 137
144 ColumnVector tmp_fval = user_csub (tmp_state, t); 138 ColumnVector tmp_fval = (*user_csub) (tmp_state, t);
145 139
146 for (int i = 0; i < ng; i++) 140 for (int i = 0; i < ng; i++)
147 gout[i] = tmp_fval(i); 141 gout[i] = tmp_fval(i);
148 142
149 return 0; 143 return 0;
150 } 144 }
151 145
152
153 DASRT::DASRT (void) 146 DASRT::DASRT (void)
154 : DAERT () 147 : DAERT ()
155 { 148 {
156 initialized = false; 149 initialized = false;
157 restart = false; 150 restart = false;
161 154
162 sanity_checked = false; 155 sanity_checked = false;
163 156
164 info.resize (30, 0); 157 info.resize (30, 0);
165 158
166 npar = 0;
167 ng = 0; 159 ng = 0;
168 160
169 liw = 0; 161 liw = 0;
170 lrw = 0; 162 lrw = 0;
163 }
164
165 DASRT::DASRT (const ColumnVector& state, double time, DAERTFunc& f)
166 : DAERT (state, time, f)
167 {
168 n = size ();
169
170 initialized = false;
171 restart = false;
172
173 stop_time_set = false;
174 stop_time = 0.0;
175
176 liw = 20 + n;
177 lrw = 50 + 9*n + n*n;
178
179 sanity_checked = false;
180
181 info.resize (15, 0);
182
183 DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function ();
184
185 if (tmp_csub)
186 {
187 ColumnVector tmp = tmp_csub (state, time);
188 ng = tmp.length ();
189 }
190 else
191 ng = 0;
192
193 jroot.resize (ng, 1);
171 } 194 }
172 195
173 DASRT::DASRT (const ColumnVector& state, const ColumnVector& deriv, 196 DASRT::DASRT (const ColumnVector& state, const ColumnVector& deriv,
174 double time, DAERTFunc& f) 197 double time, DAERTFunc& f)
175 : DAERT (state, deriv, time, f) 198 : DAERT (state, deriv, time, f)
183 stop_time = 0.0; 206 stop_time = 0.0;
184 207
185 sanity_checked = false; 208 sanity_checked = false;
186 209
187 info.resize (30, 0); 210 info.resize (30, 0);
188 jroot.resize (ng, 1);
189
190 npar = 0;
191 211
192 DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function (); 212 DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function ();
193 213
194 if (tmp_csub) 214 if (tmp_csub)
195 { 215 {
197 ng = tmp.length (); 217 ng = tmp.length ();
198 } 218 }
199 else 219 else
200 ng = 0; 220 ng = 0;
201 221
202 rpar.resize (npar+1); 222 liw = 20 + n + 1000;
203 ipar.resize (npar+1); 223 lrw = 50 + 9*n + n*n + 1000;
204 224
205 info(11) = npar; 225 jroot.resize (ng, 1);
206
207 // Also store it here, for communication with user-supplied
208 // subroutines.
209 ipar(0) = npar;
210
211 y.resize (n, 1, 0.0);
212 ydot.resize (n, 1, 0.0);
213 }
214
215 void
216 DASRT::init_work_size (int info_zero)
217 {
218 double t;
219 double *py = y.fortran_vec ();
220 double *pydot = ydot.fortran_vec ();
221 double rel_tol = relative_tolerance ();
222 double abs_tol = absolute_tolerance ();
223 int *pinfo = info.fortran_vec ();
224 double *prpar = rpar.fortran_vec ();
225 int *pipar = ipar.fortran_vec ();
226 int *pjroot = jroot.fortran_vec ();
227 int idid;
228
229 // We do not have to lie.
230 rwork.resize (5000+9*n+n*n, 0.0);
231 iwork.resize (n+20, 0);
232
233 liw = n+20;
234 lrw = 5000+9*n+n*n;
235
236 double *prwork = rwork.fortran_vec ();
237 int *piwork = iwork.fortran_vec ();
238
239 F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo,
240 &rel_tol, &abs_tol, idid, prwork, lrw,
241 piwork, liw, prpar, pipar, ddasrt_j,
242 ddasrt_g, ng, pjroot);
243
244 int iwadd = iwork(18);
245
246 if (iwadd > 0)
247 liw += iwadd;
248
249 info(0) = 0;
250
251 iwork.resize (liw, 0);
252
253 piwork = iwork.fortran_vec ();
254
255 F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo,
256 &rel_tol, &abs_tol, idid, prwork, lrw,
257 piwork, liw, prpar, pipar, ddasrt_j,
258 ddasrt_g, ng, pjroot);
259
260 int rwadd = iwork(19);
261
262 if (rwadd > 0)
263 lrw += rwadd;
264
265 rwork.resize (lrw, 0.0);
266
267 info(0) = info_zero;
268 } 226 }
269 227
270 void 228 void
271 DASRT::force_restart (void) 229 DASRT::force_restart (void)
272 { 230 {
294 252
295 if (! initialized) 253 if (! initialized)
296 { 254 {
297 info(0) = 0; 255 info(0) = 0;
298 256
299 for (int i = 0; i < n; i++)
300 {
301 y(i,0) = x(i);
302 ydot(i,0) = xdot(i);
303 }
304
305 integration_error = false; 257 integration_error = false;
306
307 nn = n;
308 258
309 user_fsub = DAEFunc::function (); 259 user_fsub = DAEFunc::function ();
310 user_jsub = DAEFunc::jacobian_function (); 260 user_jsub = DAEFunc::jacobian_function ();
311 user_csub = DAERTFunc::constraint_function (); 261 user_csub = DAERTFunc::constraint_function ();
312 262
313 if (user_jsub) 263 if (user_jsub)
314 info(4) = 1; 264 info(4) = 1;
315 else 265 else
316 info(4) = 0; 266 info(4) = 0;
317 267
268 px = x.fortran_vec ();
269 pxdot = xdot.fortran_vec ();
270
271 nn = n;
272
318 if (! sanity_checked) 273 if (! sanity_checked)
319 { 274 {
320 int ires = 0; 275 int ires = 0;
321 276
322 ColumnVector fval = user_fsub (x, xdot, t, ires); 277 ColumnVector fval = (*user_fsub) (x, xdot, t, ires);
323 278
324 if (fval.length () != x.length ()) 279 if (fval.length () != x.length ())
325 { 280 {
326 (*current_liboctave_error_handler) 281 (*current_liboctave_error_handler)
327 ("dassl: inconsistent sizes for state and residual vectors"); 282 ("dassl: inconsistent sizes for state and residual vectors");
331 } 286 }
332 287
333 sanity_checked = true; 288 sanity_checked = true;
334 } 289 }
335 290
336 init_work_size (info(0));
337
338 if (iwork.length () != liw) 291 if (iwork.length () != liw)
339 iwork.resize (liw); 292 iwork.resize (liw);
340 293
341 if (rwork.length () != lrw) 294 if (rwork.length () != lrw)
342 rwork.resize (lrw); 295 rwork.resize (lrw);
366 info(6) = 1; 319 info(6) = 1;
367 } 320 }
368 else 321 else
369 info(6) = 0; 322 info(6) = 0;
370 323
371 py = y.fortran_vec ();
372 pydot = ydot.fortran_vec ();
373 pinfo = info.fortran_vec (); 324 pinfo = info.fortran_vec ();
374 piwork = iwork.fortran_vec (); 325 piwork = iwork.fortran_vec ();
375 prwork = rwork.fortran_vec (); 326 prwork = rwork.fortran_vec ();
376 prpar = rpar.fortran_vec ();
377 pipar = ipar.fortran_vec ();
378 pjroot = jroot.fortran_vec (); 327 pjroot = jroot.fortran_vec ();
379 328
380 info(5) = 0; 329 info(5) = 0;
381 info(8) = 0; 330 info(8) = 0;
382 initialized = true; 331 initialized = true;
393 } 342 }
394 else 343 else
395 info(3) = 0; 344 info(3) = 0;
396 } 345 }
397 346
398 F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, py, pydot, tout, pinfo, 347 double *dummy = 0;
348 int *idummy = 0;
349
350 F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, px, pxdot, tout, pinfo,
399 &rel_tol, &abs_tol, idid, prwork, lrw, 351 &rel_tol, &abs_tol, idid, prwork, lrw,
400 piwork, liw, prpar, pipar, ddasrt_j, 352 piwork, liw, dummy, idummy, ddasrt_j,
401 ddasrt_g, ng, pjroot)); 353 ddasrt_g, ng, pjroot));
402 354
403 if (f77_exception_encountered) 355 if (f77_exception_encountered)
404 { 356 {
405 integration_error = true; 357 integration_error = true;
418 // (T=TOUT) by stepping past TOUT. Y(*) is obtained by 370 // (T=TOUT) by stepping past TOUT. Y(*) is obtained by
419 // interpolation. YPRIME(*) is obtained by interpolation. 371 // interpolation. YPRIME(*) is obtained by interpolation.
420 case 5: // The integration to TSTOP was successfully completed 372 case 5: // The integration to TSTOP was successfully completed
421 // (T=TSTOP) by stepping to TSTOP within the 373 // (T=TSTOP) by stepping to TSTOP within the
422 // tolerance. Must restart to continue. 374 // tolerance. Must restart to continue.
423 for (int i = 0; i < n; i++)
424 x(i) = y(i,0);
425 t = tout; 375 t = tout;
426 break; 376 break;
427 377
428 case 4: // We've hit the stopping condition. 378 case 4: // We've hit the stopping condition.
429 for (int i = 0; i < n; i++)
430 x(i) = y(i,0);
431 break; 379 break;
432 380
433 case -1: // A large amount of work has been expended. (~500 steps). 381 case -1: // A large amount of work has been expended. (~500 steps).
434 case -2: // The error tolerances are too stringent. 382 case -2: // The error tolerances are too stringent.
435 case -3: // The local error test cannot be satisfied because you 383 case -3: // The local error test cannot be satisfied because you
466 { 414 {
467 DASRT_result retval; 415 DASRT_result retval;
468 416
469 Matrix x_out; 417 Matrix x_out;
470 Matrix xdot_out; 418 Matrix xdot_out;
471 ColumnVector t_out; 419 ColumnVector t_out = tout;
472 420
473 int oldj = 0; 421 int oldj = 0;
474 422
475 int n_out = tout.capacity (); 423 int n_out = tout.capacity ();
476 424
477 if (n_out > 0 && n > 0) 425 if (n_out > 0 && n > 0)
478 { 426 {
479 x_out.resize (n_out, n); 427 x_out.resize (n_out, n);
480 xdot_out.resize (n_out, n); 428 xdot_out.resize (n_out, n);
481 t_out.resize (n_out); 429
482 430 for (int i = 0; i < n; i++)
483 for (int j = 0; j < n_out; j++) 431 {
432 x_out(0,i) = x(i);
433 xdot_out(0,i) = xdot(i);
434 }
435
436 for (int j = 1; j < n_out; j++)
484 { 437 {
485 integrate (tout(j)); 438 integrate (tout(j));
439
486 if (integration_error) 440 if (integration_error)
487 { 441 {
488 retval = DASRT_result (x_out, xdot_out, t_out); 442 retval = DASRT_result (x_out, xdot_out, t_out);
489 return retval; 443 return retval;
490 } 444 }
494 else 448 else
495 t_out(j) = tout(j); 449 t_out(j) = tout(j);
496 450
497 for (int i = 0; i < n; i++) 451 for (int i = 0; i < n; i++)
498 { 452 {
499 x_out(j,i) = y(i,0); 453 x_out(j,i) = x(i);
500 xdot_out(j,i) = ydot(i,0); 454 xdot_out(j,i) = xdot(i);
501 } 455 }
502 456
503 if (idid == 4) 457 if (idid == 4)
504 { 458 {
505 oldj = j; 459 x_out.resize (j+1, n);
506 j = n_out; 460 xdot_out.resize (j+1, n);
507 x_out.resize (oldj+1, n); 461 t_out.resize (j+1);
508 xdot_out.resize (oldj+1, n); 462 break;
509 t_out.resize (oldj+1);
510 } 463 }
511 } 464 }
512 } 465 }
513 466
514 retval = DASRT_result (x_out, xdot_out, t_out); 467 retval = DASRT_result (x_out, xdot_out, t_out);
521 { 474 {
522 DASRT_result retval; 475 DASRT_result retval;
523 476
524 Matrix x_out; 477 Matrix x_out;
525 Matrix xdot_out; 478 Matrix xdot_out;
526 ColumnVector t_outs; 479 ColumnVector t_outs = tout;
527 480
528 int n_out = tout.capacity (); 481 int n_out = tout.capacity ();
529 482
530 if (n_out > 0 && n > 0) 483 if (n_out > 0 && n > 0)
531 { 484 {
532 x_out.resize (n_out, n); 485 x_out.resize (n_out, n);
533 xdot_out.resize (n_out, n); 486 xdot_out.resize (n_out, n);
534 t_outs.resize (n_out);
535 487
536 int n_crit = tcrit.capacity (); 488 int n_crit = tcrit.capacity ();
537 489
538 if (n_crit > 0) 490 if (n_crit > 0)
539 { 491 {
540 int i_crit = 0; 492 int i_crit = 0;
541 int i_out = 0; 493 int i_out = 1;
542 double next_crit = tcrit(0); 494 double next_crit = tcrit(0);
543 double next_out; 495 double next_out;
544 while (i_out < n_out) 496 while (i_out < n_out)
545 { 497 {
546 bool do_restart = false; 498 bool do_restart = false;
600 552
601 if (save_output) 553 if (save_output)
602 { 554 {
603 for (int i = 0; i < n; i++) 555 for (int i = 0; i < n; i++)
604 { 556 {
605 x_out(i_out-1,i) = y(i,0); 557 x_out(i_out-1,i) = x(i);
606 xdot_out(i_out-1,i) = ydot(i,0); 558 xdot_out(i_out-1,i) = xdot(i);
607 } 559 }
560
608 t_outs(i_out-1) = t_out; 561 t_outs(i_out-1) = t_out;
609 if (idid ==4) 562
563 if (idid == 4)
610 { 564 {
611 x_out.resize (i_out, n); 565 x_out.resize (i_out, n);
612 xdot_out.resize (i_out, n); 566 xdot_out.resize (i_out, n);
613 t_outs.resize (i_out); 567 t_outs.resize (i_out);
614 i_out = n_out; 568 i_out = n_out;