Mercurial > hg > octave-nkf

comparison liboctave/DASSL.cc @ 1842:0574a1f3a273

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-02-03 11:44:02 by jwe]
author jwe
date Sat, 03 Feb 1996 11:44:02 +0000
parents a272c4056bab
children 2ffe49eb95a5
comparison
equal deleted inserted replaced
1841:fc5667a20dd2 1842:0574a1f3a273
1 // DAE.cc -*- C++ -*- 1 // DASSL.cc -*- C++ -*-
2 /* 2 /*
3 3
4 Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton 4 Copyright (C) 1996 John W. Eaton
5 5
6 This file is part of Octave. 6 This file is part of Octave.
7 7
8 Octave is free software; you can redistribute it and/or modify it 8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the 9 under the terms of the GNU General Public License as published by the
27 27
28 #ifdef HAVE_CONFIG_H 28 #ifdef HAVE_CONFIG_H
29 #include <config.h> 29 #include <config.h>
30 #endif 30 #endif
31 31
32 #include "DAE.h" 32 #include <cfloat>
33 #include <cmath>
34
35 #include "DASSL.h"
33 #include "f77-uscore.h" 36 #include "f77-uscore.h"
34 #include "lo-error.h" 37 #include "lo-error.h"
35 38
36 extern "C" 39 extern "C"
37 { 40 {
50 53
51 static DAEFunc::DAERHSFunc user_fun; 54 static DAEFunc::DAERHSFunc user_fun;
52 static DAEFunc::DAEJacFunc user_jac; 55 static DAEFunc::DAEJacFunc user_jac;
53 static int nn; 56 static int nn;
54 57
55 DAE::DAE (void) 58 DASSL::DASSL (void) : DAE ()
56 { 59 {
57 n = 0;
58 t = 0.0;
59
60 stop_time_set = 0; 60 stop_time_set = 0;
61 stop_time = 0.0; 61 stop_time = 0.0;
62
63 integration_error = 0;
64 restart = 1;
65
66 DAEFunc::set_function (0);
67 DAEFunc::set_jacobian_function (0);
68 62
69 liw = 0; 63 liw = 0;
70 lrw = 0; 64 lrw = 0;
71 65
72 info = new int [15]; 66 info = new int [15];
75 69
76 for (int i = 0; i < 15; i++) 70 for (int i = 0; i < 15; i++)
77 info [i] = 0; 71 info [i] = 0;
78 } 72 }
79 73
80 DAE::DAE (int size) 74 DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f)
81 { 75 : DAE (state, time, f)
82 n = size; 76 {
83 t = 0.0; 77 n = size ();
84 78
85 stop_time_set = 0; 79 stop_time_set = 0;
86 stop_time = 0.0; 80 stop_time = 0.0;
87
88 integration_error = 0;
89 restart = 1;
90
91 DAEFunc::set_function (0);
92 DAEFunc::set_jacobian_function (0);
93 81
94 liw = 20 + n; 82 liw = 20 + n;
95 lrw = 40 + 9*n + n*n; 83 lrw = 40 + 9*n + n*n;
96 84
97 info = new int [15]; 85 info = new int [15];
100 88
101 for (int i = 0; i < 15; i++) 89 for (int i = 0; i < 15; i++)
102 info [i] = 0; 90 info [i] = 0;
103 } 91 }
104 92
105 DAE::DAE (const ColumnVector& state, double time, DAEFunc& f) 93 DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv,
106 { 94 double time, DAEFunc& f)
107 n = state.capacity (); 95 : DAE (state, deriv, time, f)
108 t = time; 96 {
109 x = state; 97 n = size ();
110 xdot.resize (n, 0.0);
111 98
112 stop_time_set = 0; 99 stop_time_set = 0;
113 stop_time = 0.0; 100 stop_time = 0.0;
114
115 integration_error = 0;
116 restart = 1;
117 101
118 DAEFunc::set_function (f.function ()); 102 DAEFunc::set_function (f.function ());
119 DAEFunc::set_jacobian_function (f.jacobian_function ()); 103 DAEFunc::set_jacobian_function (f.jacobian_function ());
120 104
121 liw = 20 + n; 105 liw = 20 + n;
127 111
128 for (int i = 0; i < 15; i++) 112 for (int i = 0; i < 15; i++)
129 info [i] = 0; 113 info [i] = 0;
130 } 114 }
131 115
132 DAE::DAE (const ColumnVector& state, const ColumnVector& deriv, 116 DASSL::~DASSL (void)
133 double time, DAEFunc& f)
134 {
135 if (deriv.capacity () != state.capacity ())
136 {
137 (*current_liboctave_error_handler)
138 ("x, xdot size mismatch in DAE constructor");
139 n = 0;
140 t = 0.0;
141 return;
142 }
143
144 n = state.capacity ();
145 t = time;
146 xdot = deriv;
147 x = state;
148
149 stop_time_set = 0;
150 stop_time = 0.0;
151
152 DAEFunc::set_function (f.function ());
153 DAEFunc::set_jacobian_function (f.jacobian_function ());
154
155 liw = 20 + n;
156 lrw = 40 + 9*n + n*n;
157
158 info = new int [15];
159 iwork = new int [liw];
160 rwork = new double [lrw];
161
162 for (int i = 0; i < 15; i++)
163 info [i] = 0;
164 }
165
166 DAE::~DAE (void)
167 { 117 {
168 delete info; 118 delete info;
169 delete rwork; 119 delete rwork;
170 delete iwork; 120 delete iwork;
171 } 121 }
172 122
173 ColumnVector 123 void
174 DAE::deriv (void) 124 DASSL::force_restart (void)
175 { 125 {
176 return xdot; 126 restart = 1;
177 }
178
179 void
180 DAE::initialize (const ColumnVector& state, double time)
181 {
182 integration_error = 0; 127 integration_error = 0;
183 restart = 1; 128 }
184 x = state; 129
185 int nx = x.capacity (); 130 void
186 xdot.resize (nx, 0.0); 131 DASSL::set_stop_time (double t)
187 t = time; 132 {
188 } 133 stop_time_set = 1;
189 134 stop_time = t;
190 void 135 }
191 DAE::initialize (const ColumnVector& state, 136
192 const ColumnVector& deriv, double time) 137 void
193 { 138 DASSL::clear_stop_time (void)
194 integration_error = 0; 139 {
195 restart = 1; 140 stop_time_set = 0;
196 xdot = deriv;
197 x = state;
198 t = time;
199 } 141 }
200 142
201 int 143 int
202 ddassl_f (const double& time, double *state, double *deriv, 144 ddassl_f (const double& time, double *state, double *deriv,
203 double *delta, int& ires, double *, int *) 145 double *delta, int& ires, double *, int *)
253 195
254 return 0; 196 return 0;
255 } 197 }
256 198
257 ColumnVector 199 ColumnVector
258 DAE::integrate (double tout) 200 DASSL::do_integrate (double tout)
259 { 201 {
260 integration_error = 0; 202 integration_error = 0;
261 203
262 if (DAEFunc::jacobian_function ()) 204 if (DAEFunc::jacobian_function ())
263 iwork [4] = 1; 205 iwork [4] = 1;
361 303
362 return x; 304 return x;
363 } 305 }
364 306
365 Matrix 307 Matrix
366 DAE::integrate (const ColumnVector& tout, Matrix& xdot_out) 308 DASSL::do_integrate (const ColumnVector& tout)
309 {
310 Matrix dummy;
311 return integrate (tout, dummy);
312 }
313
314 Matrix
315 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
367 { 316 {
368 Matrix retval; 317 Matrix retval;
369 int n_out = tout.capacity (); 318 int n_out = tout.capacity ();
370 319
371 if (n_out > 0 && n > 0) 320 if (n_out > 0 && n > 0)
379 xdot_out.elem (0, i) = xdot.elem (i); 328 xdot_out.elem (0, i) = xdot.elem (i);
380 } 329 }
381 330
382 for (int j = 1; j < n_out; j++) 331 for (int j = 1; j < n_out; j++)
383 { 332 {
384 ColumnVector x_next = integrate (tout.elem (j)); 333 ColumnVector x_next = do_integrate (tout.elem (j));
385 334
386 if (integration_error) 335 if (integration_error)
387 return retval; 336 return retval;
388 337
389 for (int i = 0; i < n; i++) 338 for (int i = 0; i < n; i++)
396 345
397 return retval; 346 return retval;
398 } 347 }
399 348
400 Matrix 349 Matrix
401 DAE::integrate (const ColumnVector& tout, Matrix& xdot_out, 350 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out,
402 const ColumnVector& tcrit) 351 const ColumnVector& tcrit)
403 { 352 {
404 Matrix retval; 353 Matrix retval;
405 int n_out = tout.capacity (); 354 int n_out = tout.capacity ();
406 355
407 if (n_out > 0 && n > 0) 356 if (n_out > 0 && n > 0)
467 t_out = next_out; 416 t_out = next_out;
468 save_output = 1; 417 save_output = 1;
469 i_out++; 418 i_out++;
470 } 419 }
471 420
472 ColumnVector x_next = integrate (t_out); 421 ColumnVector x_next = do_integrate (t_out);
473 422
474 if (integration_error) 423 if (integration_error)
475 return retval; 424 return retval;
476 425
477 if (save_output) 426 if (save_output)
497 } 446 }
498 447
499 return retval; 448 return retval;
500 } 449 }
501 450
451 DASSL_options::DASSL_options (void)
452 {
453 init ();
454 }
455
456 DASSL_options::DASSL_options (const DASSL_options& opt)
457 {
458 copy (opt);
459 }
460
461 DASSL_options&
462 DASSL_options::operator = (const DASSL_options& opt)
463 {
464 if (this != &opt)
465 copy (opt);
466
467 return *this;
468 }
469
470 DASSL_options::~DASSL_options (void)
471 {
472 }
473
474 void
475 DASSL_options::init (void)
476 {
477 double sqrt_eps = sqrt (DBL_EPSILON);
478 x_absolute_tolerance = sqrt_eps;
479 x_initial_step_size = -1.0;
480 x_maximum_step_size = -1.0;
481 x_minimum_step_size = 0.0;
482 x_relative_tolerance = sqrt_eps;
483 }
484
485 void
486 DASSL_options::copy (const DASSL_options& opt)
487 {
488 x_absolute_tolerance = opt.x_absolute_tolerance;
489 x_initial_step_size = opt.x_initial_step_size;
490 x_maximum_step_size = opt.x_maximum_step_size;
491 x_minimum_step_size = opt.x_minimum_step_size;
492 x_relative_tolerance = opt.x_relative_tolerance;
493 }
494
495 void
496 DASSL_options::set_default_options (void)
497 {
498 init ();
499 }
500
501 void
502 DASSL_options::set_absolute_tolerance (double val)
503 {
504 x_absolute_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
505 }
506
507 void
508 DASSL_options::set_initial_step_size (double val)
509 {
510 x_initial_step_size = (val >= 0.0) ? val : -1.0;
511 }
512
513 void
514 DASSL_options::set_maximum_step_size (double val)
515 {
516 x_maximum_step_size = (val >= 0.0) ? val : -1.0;
517 }
518
519 void
520 DASSL_options::set_minimum_step_size (double val)
521 {
522 x_minimum_step_size = (val >= 0.0) ? val : 0.0;
523 }
524
525 void
526 DASSL_options::set_relative_tolerance (double val)
527 {
528 x_relative_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
529 }
530
531 double
532 DASSL_options::absolute_tolerance (void)
533 {
534 return x_absolute_tolerance;
535 }
536
537 double
538 DASSL_options::initial_step_size (void)
539 {
540 return x_initial_step_size;
541 }
542
543 double
544 DASSL_options::maximum_step_size (void)
545 {
546 return x_maximum_step_size;
547 }
548
549 double
550 DASSL_options::minimum_step_size (void)
551 {
552 return x_minimum_step_size;
553 }
554
555 double
556 DASSL_options::relative_tolerance (void)
557 {
558 return x_relative_tolerance;
559 }
560
502 /* 561 /*
503 ;;; Local Variables: *** 562 ;;; Local Variables: ***
504 ;;; mode: C++ *** 563 ;;; mode: C++ ***
505 ;;; page-delimiter: "^/\\*" *** 564 ;;; page-delimiter: "^/\\*" ***
506 ;;; End: *** 565 ;;; End: ***