Mercurial > hg > octave-nkf

comparison liboctave/DASRT.cc @ 4133:402d7b86a0a2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-10-29 23:57:34 by jwe]
author jwe
date Tue, 29 Oct 2002 23:57:34 +0000
parents 47d3baea432d
children 84fe3ca3a246
comparison
equal deleted inserted replaced
4132:87eb044020ae 4133:402d7b86a0a2
62 static int 62 static int
63 ddasrt_f (const double& t, const double *state, const double *deriv, 63 ddasrt_f (const double& t, const double *state, const double *deriv,
64 double *delta, int& ires, double *rpar, int *ipar) 64 double *delta, int& ires, double *rpar, int *ipar)
65 { 65 {
66 ColumnVector tmp_state (nn); 66 ColumnVector tmp_state (nn);
67 ColumnVector tmp_deriv (nn);
68
67 for (int i = 0; i < nn; i++) 69 for (int i = 0; i < nn; i++)
68 tmp_state(i) = state[i]; 70 {
69 71 tmp_state(i) = state[i];
70 ColumnVector tmp_deriv (nn); 72 tmp_deriv(i) = deriv[i];
71 for (int i = 0; i < nn; i++) 73 }
72 tmp_deriv(i) = deriv[i];
73 74
74 ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires); 75 ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires);
75 76
76 if (tmp_fval.length () == 0) 77 if (tmp_fval.length () == 0)
77 ires = -2; 78 ires = -2;
151 152
152 int n = size (); 153 int n = size ();
153 154
154 nn = n; 155 nn = n;
155 156
156 liw = 20 + n;
157 lrw = 50 + 9*n + n*n;
158
159 iwork.resize (liw);
160 rwork.resize (lrw);
161
162 info(0) = 0;
163
164 if (stop_time_set)
165 {
166 info(3) = 1;
167 rwork(0) = stop_time;
168 }
169 else
170 info(3) = 0;
171
172 px = x.fortran_vec ();
173 pxdot = xdot.fortran_vec ();
174
175 piwork = iwork.fortran_vec ();
176 prwork = rwork.fortran_vec ();
177
178 restart = false;
179
180 // DAEFunc
181
182 user_fsub = DAEFunc::function ();
183 user_jsub = DAEFunc::jacobian_function ();
184
185 if (user_fsub)
186 {
187 int ires = 0;
188
189 ColumnVector fval = (*user_fsub) (x, xdot, t, ires);
190
191 if (fval.length () != x.length ())
192 {
193 (*current_liboctave_error_handler)
194 ("dasrt: inconsistent sizes for state and residual vectors");
195
196 integration_error = true;
197 return;
198 }
199 }
200 else
201 {
202 (*current_liboctave_error_handler)
203 ("dasrt: no user supplied RHS subroutine!");
204
205 integration_error = true;
206 return;
207 }
208
209 info(4) = user_jsub ? 1 : 0;
210
211 DAEFunc::reset = false;
212
213 // DAERTFunc 157 // DAERTFunc
214 158
215 user_csub = DAERTFunc::constraint_function (); 159 user_csub = DAERTFunc::constraint_function ();
216 160
217 if (user_csub) 161 if (user_csub)
219 ColumnVector tmp = (*user_csub) (x, t); 163 ColumnVector tmp = (*user_csub) (x, t);
220 ng = tmp.length (); 164 ng = tmp.length ();
221 } 165 }
222 else 166 else
223 ng = 0; 167 ng = 0;
224
225 jroot.resize (ng, 1);
226
227 pjroot = jroot.fortran_vec ();
228
229 DAERTFunc::reset = false;
230
231 // DASRT_options
232
233 double mss = maximum_step_size ();
234 if (mss >= 0.0)
235 {
236 rwork(1) = mss;
237 info(6) = 1;
238 }
239 else
240 info(6) = 0;
241
242 double iss = initial_step_size ();
243 if (iss >= 0.0)
244 {
245 rwork(2) = iss;
246 info(7) = 1;
247 }
248 else
249 info(7) = 0;
250 168
251 int maxord = maximum_order (); 169 int maxord = maximum_order ();
252 if (maxord >= 0) 170 if (maxord >= 0)
253 { 171 {
254 if (maxord > 0 && maxord < 6) 172 if (maxord > 0 && maxord < 6)
262 ("dassl: invalid value for maximum order"); 180 ("dassl: invalid value for maximum order");
263 integration_error = true; 181 integration_error = true;
264 return; 182 return;
265 } 183 }
266 } 184 }
185
186 liw = 20 + n;
187 lrw = 50 + 9*n + n*n + 3*ng;
188
189 iwork.resize (liw);
190 rwork.resize (lrw);
191
192 info(0) = 0;
193
194 if (stop_time_set)
195 {
196 info(3) = 1;
197 rwork(0) = stop_time;
198 }
199 else
200 info(3) = 0;
201
202 px = x.fortran_vec ();
203 pxdot = xdot.fortran_vec ();
204
205 piwork = iwork.fortran_vec ();
206 prwork = rwork.fortran_vec ();
207
208 restart = false;
209
210 // DAEFunc
211
212 user_fsub = DAEFunc::function ();
213 user_jsub = DAEFunc::jacobian_function ();
214
215 if (user_fsub)
216 {
217 int ires = 0;
218
219 ColumnVector fval = (*user_fsub) (x, xdot, t, ires);
220
221 if (fval.length () != x.length ())
222 {
223 (*current_liboctave_error_handler)
224 ("dasrt: inconsistent sizes for state and residual vectors");
225
226 integration_error = true;
227 return;
228 }
229 }
230 else
231 {
232 (*current_liboctave_error_handler)
233 ("dasrt: no user supplied RHS subroutine!");
234
235 integration_error = true;
236 return;
237 }
238
239 info(4) = user_jsub ? 1 : 0;
240
241 DAEFunc::reset = false;
242
243 jroot.resize (ng, 1);
244
245 pjroot = jroot.fortran_vec ();
246
247 DAERTFunc::reset = false;
248
249 // DASRT_options
250
251 double mss = maximum_step_size ();
252 if (mss >= 0.0)
253 {
254 rwork(1) = mss;
255 info(6) = 1;
256 }
257 else
258 info(6) = 0;
259
260 double iss = initial_step_size ();
261 if (iss >= 0.0)
262 {
263 rwork(2) = iss;
264 info(7) = 1;
265 }
266 else
267 info(7) = 0;
267 268
268 if (step_limit () >= 0) 269 if (step_limit () >= 0)
269 { 270 {
270 info(11) = 1; 271 info(11) = 1;
271 iwork(18) = step_limit (); 272 iwork(18) = step_limit ();