Mercurial > hg > octave-nkf

comparison liboctave/DASSL.cc @ 4049:a35a3c5d4740

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-08-16 08:54:31 by jwe]
author jwe
date Fri, 16 Aug 2002 08:54:31 +0000
parents 7b0c139ac8af
children b79da8779a0e
comparison
equal deleted inserted replaced
4048:c9991c59cf5c 4049:a35a3c5d4740
54 static DAEFunc::DAERHSFunc user_fun; 54 static DAEFunc::DAERHSFunc user_fun;
55 static DAEFunc::DAEJacFunc user_jac; 55 static DAEFunc::DAEJacFunc user_jac;
56 56
57 static int nn; 57 static int nn;
58 58
59 DASSL::DASSL (void) : DAE () 59 static int
60 {
61 liw = 0;
62 lrw = 0;
63
64 sanity_checked = false;
65
66 info.resize (15);
67
68 for (int i = 0; i < 15; i++)
69 info.elem (i) = 0;
70 }
71
72 DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f)
73 : DAE (state, time, f)
74 {
75 n = size ();
76
77 liw = 20 + n;
78 lrw = 40 + 9*n + n*n;
79
80 sanity_checked = false;
81
82 info.resize (15, 0);
83 }
84
85 DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv,
86 double time, DAEFunc& f)
87 : DAE (state, deriv, time, f)
88 {
89 n = size ();
90
91 DAEFunc::set_function (f.function ());
92 DAEFunc::set_jacobian_function (f.jacobian_function ());
93
94 liw = 20 + n;
95 lrw = 40 + 9*n + n*n;
96
97 sanity_checked = false;
98
99 info.resize (15);
100
101 for (int i = 0; i < 15; i++)
102 info.elem (i) = 0;
103 }
104
105 int
106 ddassl_f (const double& time, const double *state, const double *deriv, 60 ddassl_f (const double& time, const double *state, const double *deriv,
107 double *delta, int& ires, double *, int *) 61 double *delta, int& ires, double *, int *)
108 { 62 {
109 // XXX FIXME XXX -- would be nice to avoid copying the data. 63 // XXX FIXME XXX -- would be nice to avoid copying the data.
110 64
132 } 86 }
133 87
134 return 0; 88 return 0;
135 } 89 }
136 90
137 int 91 static int
138 ddassl_j (const double& time, const double *state, const double *deriv, 92 ddassl_j (const double& time, const double *state, const double *deriv,
139 double *pd, const double& cj, double *, int *) 93 double *pd, const double& cj, double *, int *)
140 { 94 {
141 // XXX FIXME XXX -- would be nice to avoid copying the data. 95 // XXX FIXME XXX -- would be nice to avoid copying the data.
142 96
159 } 113 }
160 114
161 ColumnVector 115 ColumnVector
162 DASSL::do_integrate (double tout) 116 DASSL::do_integrate (double tout)
163 { 117 {
164 // XXX FIXME XXX -- should handle all this option stuff just once
165 // for each new problem.
166
167 ColumnVector retval; 118 ColumnVector retval;
168 119
169 if (restart) 120 if (! initialized || restart || DAEFunc::reset|| DASSL_options::reset)
170 { 121 {
122 integration_error = false;
123
124 initialized = true;
125
126 info.resize (15);
127
128 for (int i = 0; i < 15; i++)
129 info(i) = 0;
130
131 pinfo = info.fortran_vec ();
132
133 int n = size ();
134
135 liw = 20 + n;
136 lrw = 40 + 9*n + n*n;
137
138 nn = n;
139
140 iwork.resize (liw);
141 rwork.resize (lrw);
142
143 info(0) = 0;
144
145 if (stop_time_set)
146 {
147 rwork(0) = stop_time;
148 info(3) = 1;
149 }
150 else
151 info(3) = 0;
152
153 px = x.fortran_vec ();
154 pxdot = xdot.fortran_vec ();
155
156 piwork = iwork.fortran_vec ();
157 prwork = rwork.fortran_vec ();
158
171 restart = false; 159 restart = false;
172 info.elem (0) = 0; 160
173 } 161 // DAEFunc
174 162
175 if (iwork.length () != liw) 163 user_fun = DAEFunc::function ();
176 iwork.resize (liw); 164 user_jac = DAEFunc::jacobian_function ();
177 165
178 if (rwork.length () != lrw) 166 if (user_fun)
179 rwork.resize (lrw); 167 {
180 168 int ires = 0;
181 integration_error = false; 169
182 170 ColumnVector res = (*user_fun) (x, xdot, t, ires);
183 user_fun = DAEFunc::fun; 171
184 user_jac = DAEFunc::jac; 172 if (res.length () != x.length ())
185 173 {
186 if (user_jac) 174 (*current_liboctave_error_handler)
187 info.elem (4) = 1; 175 ("dassl: inconsistent sizes for state and residual vectors");
188 else 176
189 info.elem (4) = 0; 177 integration_error = true;
190 178 return retval;
191 double *px = x.fortran_vec (); 179 }
192 double *pxdot = xdot.fortran_vec (); 180 }
193 181 else
194 nn = n;
195
196 if (! sanity_checked)
197 {
198 int ires = 0;
199
200 ColumnVector res = (*user_fun) (x, xdot, t, ires);
201
202 if (res.length () != x.length ())
203 { 182 {
204 (*current_liboctave_error_handler) 183 (*current_liboctave_error_handler)
205 ("dassl: inconsistent sizes for state and residual vectors"); 184 ("dassl: no user supplied RHS subroutine!");
206 185
207 integration_error = true; 186 integration_error = true;
208 return retval; 187 return retval;
209 } 188 }
210 189
211 sanity_checked = true; 190 info(4) = user_jac ? 1 : 0;
212 }
213 191
214 if (stop_time_set) 192 DAEFunc::reset = false;
215 { 193
216 rwork.elem (0) = stop_time; 194 // DASSL_options
217 info.elem (3) = 1; 195
218 } 196 double hmax = maximum_step_size ();
219 else 197 if (hmax >= 0.0)
220 info.elem (3) = 0; 198 {
221 199 rwork(1) = hmax;
222 Array<double> abs_tol = absolute_tolerance (); 200 info(6) = 1;
223 Array<double> rel_tol = relative_tolerance (); 201 }
224 202 else
225 int abs_tol_len = abs_tol.length (); 203 info(6) = 0;
226 int rel_tol_len = rel_tol.length (); 204
227 205 double h0 = initial_step_size ();
228 if (abs_tol_len == 1 && rel_tol_len == 1) 206 if (h0 >= 0.0)
229 { 207 {
230 info.elem (1) = 0; 208 rwork(2) = h0;
231 } 209 info(7) = 1;
232 else if (abs_tol_len == n && rel_tol_len == n) 210 }
233 { 211 else
234 info.elem (1) = 1; 212 info(7) = 0;
235 } 213
236 else 214 int maxord = maximum_order ();
237 { 215 if (maxord >= 0)
238 (*current_liboctave_error_handler) 216 {
239 ("dassl: inconsistent sizes for tolerance arrays"); 217 if (maxord > 0 && maxord < 6)
240 218 {
241 integration_error = true; 219 info(8) = 1;
242 return retval; 220 iwork(2) = maxord;
243 } 221 }
244 222 else
245 double hmax = maximum_step_size (); 223 {
246 if (hmax >= 0.0) 224 (*current_liboctave_error_handler)
247 { 225 ("dassl: invalid value for maximum order");
248 rwork.elem (1) = hmax; 226 integration_error = true;
249 info.elem (6) = 1; 227 return retval;
250 } 228 }
251 else 229 }
252 info.elem (6) = 0; 230
253 231 int enc = enforce_nonnegativity_constraints ();
254 double h0 = initial_step_size (); 232 info(9) = enc ? 1 : 0;
255 if (h0 >= 0.0) 233
256 { 234 int ccic = compute_consistent_initial_condition ();
257 rwork.elem (2) = h0; 235 info(10) = ccic ? 1 : 0;
258 info.elem (7) = 1; 236
259 } 237 abs_tol = absolute_tolerance ();
260 else 238 rel_tol = relative_tolerance ();
261 info.elem (7) = 0; 239
262 240 int abs_tol_len = abs_tol.length ();
263 int maxord = maximum_order (); 241 int rel_tol_len = rel_tol.length ();
264 if (maxord >= 0) 242
265 { 243 if (abs_tol_len == 1 && rel_tol_len == 1)
266 if (maxord > 0 && maxord < 6) 244 {
267 { 245 info(1) = 0;
268 info(8) = 1; 246 }
269 iwork(2) = maxord; 247 else if (abs_tol_len == n && rel_tol_len == n)
248 {
249 info(1) = 1;
270 } 250 }
271 else 251 else
272 { 252 {
273 (*current_liboctave_error_handler) 253 (*current_liboctave_error_handler)
274 ("dassl: invalid value for maximum order"); 254 ("dassl: inconsistent sizes for tolerance arrays");
255
275 integration_error = true; 256 integration_error = true;
276 return retval; 257 return retval;
277 } 258 }
278 } 259
279 260 pabs_tol = abs_tol.fortran_vec ();
280 int enc = enforce_nonnegativity_constraints (); 261 prel_tol = rel_tol.fortran_vec ();
281 info (9) = enc ? 1 : 0; 262
282 263 DASSL_options::reset = false;
283 int ccic = compute_consistent_initial_condition (); 264 }
284 info(10) = ccic ? 1 : 0; 265
285 266 static double *dummy = 0;
286 double *dummy = 0; 267 static int *idummy = 0;
287 int *idummy = 0; 268
288 269 F77_XFCN (ddassl, DDASSL, (ddassl_f, nn, t, px, pxdot, tout, pinfo,
289 int *pinfo = info.fortran_vec ();
290 int *piwork = iwork.fortran_vec ();
291 double *prwork = rwork.fortran_vec ();
292 double *pabs_tol = abs_tol.fortran_vec ();
293 double *prel_tol = rel_tol.fortran_vec ();
294
295 // again:
296
297 F77_XFCN (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo,
298 prel_tol, pabs_tol, istate, prwork, lrw, 270 prel_tol, pabs_tol, istate, prwork, lrw,
299 piwork, liw, dummy, idummy, ddassl_j)); 271 piwork, liw, dummy, idummy, ddassl_j));
300 272
301 if (f77_exception_encountered) 273 if (f77_exception_encountered)
302 { 274 {
364 336
365 Matrix 337 Matrix
366 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out) 338 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
367 { 339 {
368 Matrix retval; 340 Matrix retval;
341
369 int n_out = tout.capacity (); 342 int n_out = tout.capacity ();
343 int n = size ();
370 344
371 if (n_out > 0 && n > 0) 345 if (n_out > 0 && n > 0)
372 { 346 {
373 retval.resize (n_out, n); 347 retval.resize (n_out, n);
374 xdot_out.resize (n_out, n); 348 xdot_out.resize (n_out, n);
407 Matrix 381 Matrix
408 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out, 382 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out,
409 const ColumnVector& tcrit) 383 const ColumnVector& tcrit)
410 { 384 {
411 Matrix retval; 385 Matrix retval;
386
412 int n_out = tout.capacity (); 387 int n_out = tout.capacity ();
388 int n = size ();
413 389
414 if (n_out > 0 && n > 0) 390 if (n_out > 0 && n > 0)
415 { 391 {
416 retval.resize (n_out, n); 392 retval.resize (n_out, n);
417 xdot_out.resize (n_out, n); 393 xdot_out.resize (n_out, n);