/[ascend]/branches/ksenija2/solvers/ida/idaanalyse.c
ViewVC logotype

Contents of /branches/ksenija2/solvers/ida/idaanalyse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2859 - (show annotations) (download) (as text)
Sun Mar 22 01:40:44 2015 UTC (7 years, 3 months ago) by jpye
File MIME type: text/x-csrc
File size: 35771 byte(s)
added debug output for relations. seems that in solardynamics.a4c, some
expected relations are not considered active; raises question about the
'analysis' routine before ida_analyse gets started?

1 /* ASCEND modelling environment
2 Copyright (C) 2006-2011 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 *//** @file
17 Analysis routines for the ASCEND wrapper of the IDA integrator.
18 These functions perform sorting of variables and relations and create
19 additional lists of variables as required for use by our ida.c code.
20 */
21 #include "idaanalyse.h"
22 #include "idatypes.h"
23 #include "idaio.h"
24
25 #include <ascend/general/panic.h>
26 #include <ascend/utilities/error.h>
27
28 #include <ascend/linear/linsolqr.h>
29
30 #include <ascend/system/diffvars.h>
31 #include <ascend/system/slv_stdcalls.h>
32 #include <ascend/system/block.h>
33 #include <ascend/system/diffvars.h>
34 #include <ascend/system/diffvars_impl.h>
35 #include <ascend/system/jacobian.h>
36 #include <ascend/system/cond_config.h>
37 #include <ascend/solver/slvDOF.h>
38
39 #define ANALYSE_DEBUG
40
41 /*
42 define DERIV_WITHOUT_DIFF to enable experimental handling of derivatives
43 for which corresponding differential vars were not found to be incident.
44
45 (this is important! it means that if we have eg "Qdot = der(Q)" in our model,
46 but Q itself isn't referred to anywhere. We want IDA to include Q in our
47 system in that case, because we want to calculate the integral of Qdot, Q,
48 as time marches forward. -- JP)
49 */
50 #define DERIV_WITHOUT_DIFF
51
52 #define VARMSG(MSG) \
53 varname = var_make_name(integ->system,v); \
54 CONSOLE_DEBUG(MSG,varname); \
55 ASC_FREE(varname)
56
57 const var_filter_t integrator_ida_nonderiv = {
58 VAR_SVAR | VAR_ACTIVE | VAR_FIXED | VAR_DERIV | VAR_NONBASIC,
59 VAR_SVAR | VAR_ACTIVE | 0 | 0 | 0
60 };
61
62 const var_filter_t integrator_ida_deriv = {
63 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | VAR_FIXED | VAR_DERIV | VAR_NONBASIC,
64 VAR_SVAR | VAR_INCIDENT | VAR_ACTIVE | 0 | VAR_DERIV | 0
65 };
66
67 const rel_filter_t integrator_ida_rel = {
68 REL_INCLUDED | REL_EQUALITY | REL_ACTIVE,
69 REL_INCLUDED | REL_EQUALITY | REL_ACTIVE
70 };
71
72
73 //static int integrator_ida_check_partitioning(IntegratorSystem *integ);
74 static int integrator_ida_check_diffindex(IntegratorSystem *integ);
75 /* static int integrator_ida_rebuild_diffindex(IntegratorSystem *integ); */
76 static int check_dups(IntegratorSystem *integ, struct var_variable **list,int n,int allownull);
77 static int integrator_ida_check_diffindex(IntegratorSystem *integ);
78 int integrator_ida_diffindex(const IntegratorSystem *integ, const struct var_variable *deriv);
79 static int integrator_ida_check_vars(IntegratorSystem *integ);
80 static int integrator_ida_flag_rels(IntegratorSystem *integ);
81 static int integrator_ida_sort_rels_and_vars(IntegratorSystem *integ);
82 static int integrator_ida_create_lists(IntegratorSystem *integ);
83 static int integrator_ida_check_index(IntegratorSystem *integ);
84
85 static int integrator_ida_vars_debug(IntegratorSystem *integ);
86 static int integrator_ida_rels_debug(IntegratorSystem *integ);
87
88 /*------------------------------------------------------------------------------
89 ANALYSIS ROUTINE (new implementation)
90 */
91
92 /**
93 Perform additional problem analysis to prepare problem for integration with
94 IDA.
95
96 We assume that the analyse_generate_diffvars routine has been called, so
97 that we just need to call slv_get_diffvars to access the derivative
98 chains.
99
100 We can also assume that the independent variable has been found.
101
102 See mailing list, ~Jan 2007.
103
104 Note, the stuff for identifying the static and output sub-problems should
105 be part of integrator.c, not this file. We will assume this is handled
106
107 @return 0 on success
108 @see integrator_analyse
109 */
110 int integrator_ida_analyse(IntegratorSystem *integ){
111 int res;
112 const SolverDiffVarCollection *diffvars;
113 int i;
114 #ifdef ANALYSE_DEBUG
115 char *varname;
116 #endif
117
118 asc_assert(integ->engine==INTEG_IDA);
119
120 #ifdef ANALYSE_DEBUG
121 CONSOLE_DEBUG("System contains a total of %d bnds and %d rels"
122 ,slv_get_num_solvers_bnds(integ->system)
123 ,slv_get_num_solvers_rels(integ->system)
124 );
125
126 /* set the active flags on variables depending on the state of WHENs */
127 CONSOLE_DEBUG("Currently %d rels active",slv_count_solvers_rels(integ->system, &integrator_ida_rel));
128 #endif
129
130 reanalyze_solver_lists_cont(integ->system);
131
132 #ifdef ANALYSE_DEBUG
133 CONSOLE_DEBUG("After analysing WHENs, there are %d rels active"
134 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
135 );
136 #endif
137
138
139 #ifdef ANALYSE_DEBUG
140 CONSOLE_DEBUG("Starting IDA analysis");
141 #endif
142
143 #if 1
144 {
145 CONSOLE_DEBUG("Setting 'nonbasic' flags...");
146 /* set the VAR_NONBASIC flag only for the indep var */
147 struct var_variable **list = slv_get_solvers_var_list(integ->system);
148 int n = slv_get_num_solvers_vars(integ->system);
149 for(i=0;i<n;++i){
150 var_set_nonbasic(list[i],0);
151 }
152 var_set_nonbasic(integ->x,1);
153 }
154 #endif
155
156 /* set the flags on differential and derivative and algebraic vars */
157 res = integrator_ida_check_vars(integ);
158 if(res){
159 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem with vars");
160 return 1;
161 }
162
163 res = integrator_ida_flag_rels(integ);
164 if(res){
165 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem classifying differential equations");
166 return 1;
167 }
168
169 #ifdef ANALYSE_DEBUG
170 CONSOLE_DEBUG("Sorting rels and vars");
171 #endif
172
173 res = integrator_ida_sort_rels_and_vars(integ);
174 if(res){
175 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem sorting rels and vars");
176 return 1;
177 }
178
179 #ifdef ANALYSE_DEBUG
180 CONSOLE_DEBUG("Creating lists");
181
182 CONSOLE_DEBUG("BEFORE MAKING LISTS");
183 integrator_ida_debug(integ,stderr);
184 #endif
185
186 res = integrator_ida_create_lists(integ);
187 if(res){
188 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem creating lists");
189 return 1;
190 }
191
192 #ifdef ANALYSE_DEBUG
193 CONSOLE_DEBUG("After ida_create_lists, there are %d rels active"
194 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
195 );
196 CONSOLE_DEBUG("Checking lists");
197 #endif
198
199 #ifdef ANALYSE_DEBUG
200 asc_assert(integ->y);
201 asc_assert(integ->ydot);
202 asc_assert(integ->y_id);
203
204 integrator_ida_vars_debug(integ);
205 integrator_ida_debug(integ,stderr);
206 #endif
207
208 if(integrator_ida_check_diffindex(integ)){
209 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error with diffindex");
210 return 360;
211 }
212
213
214 #ifdef ANALYSE_DEBUG
215 CONSOLE_DEBUG("After ida_check_diffindex, there are %d rels active"
216 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
217 );
218 #endif
219
220 /* the following causes TestIDA.testincidence3 to fail:
221 if(integ->n_y != slv_get_num_solvers_rels(integ->system)){
222 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Problem is not square");
223 return 2;
224 }*/
225
226 #if 0
227 /* WARNING: uncommenting this section seems to break the solver!!!! -- JP */
228 /* check structural singularity for the two IDACalcIC scenarios */
229
230 /* ...(1) FIX the derivatives */
231 CONSOLE_DEBUG("Checking system with derivatives fixed...");
232 for(i=0;i<integ->n_y;++i){
233 if(integ->ydot[i])var_set_fixed(integ->ydot[i],1);
234 }
235
236 slv_block_partition(integ->system);
237 res = integrator_ida_block_check(integ);
238 if(res)return 100 + res;
239
240 /* ...(2) FREE the derivatives, FIX the diffvars */
241 CONSOLE_DEBUG("Checking system with differential variables fixed...");
242 for(i=0;i<integ->n_y;++i){
243 if(integ->ydot[i]){
244 var_set_fixed(integ->ydot[i],0);
245 var_set_fixed(integ->y[i],1);
246 }
247 }
248 slv_block_partition(integ->system);
249 res = integrator_ida_block_check(integ);
250 if(res)return 200 + res;
251
252 /* ...(3) restore as it was, FREE the diffvars */
253 for(i=0;i<integ->n_y;++i){
254 if(integ->ydot[i]){
255 var_set_fixed(integ->y[i],0);
256 }
257 }
258 #endif
259
260 res = integrator_ida_check_index(integ);
261 if(res){
262 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Your DAE system has an index problem");
263 return 100 + res;
264 }
265
266 #ifdef ANALYSE_DEBUG
267 CONSOLE_DEBUG("After ida_check_index, there are %d rels active"
268 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
269 );
270 #endif
271
272
273 /* get the the dervative chains from the system */
274 diffvars = system_get_diffvars(integ->system);
275
276 /* check the indep var is same as was located earlier */
277 if(diffvars->nindep>1){
278 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Multiple variables specified as independent (ode_type=-1)");
279 return 3;
280 }else if(diffvars->nindep<1){
281 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Independent var not set (ode_type=-1)");
282 return 4;
283 }else if(diffvars->indep[0]!=integ->x){
284 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Indep var doesn't match");
285 return 5;
286 }
287
288 #ifdef ANALYSE_DEBUG
289 CONSOLE_DEBUG("Collecting observed variables");
290 #endif
291
292 /* get the observations */
293 /** @TODO this should use a 'problem provider' API instead of static data
294 from the system object */
295 integ->n_obs = diffvars->nobs;
296 integ->obs = ASC_NEW_ARRAY(struct var_variable *,integ->n_obs);
297 for(i=0;i<integ->n_obs;++i){
298 /* we get them all, regardless of flags etc */
299 integ->obs[i] = diffvars->obs[i];
300 #ifdef ANALYSE_DEBUG
301 varname = var_make_name(integ->system,integ->obs[i]);
302 CONSOLE_DEBUG("'%s' is observation",varname);
303 ASC_FREE(varname);
304 #endif
305 }
306
307 integ->n_dobs = diffvars->ndobs;
308 integ->dobs = ASC_NEW_ARRAY(struct dis_discrete *,integ->n_dobs);
309 for(i=0;i<integ->n_dobs;++i){
310 /* we get them all, regardless of flags etc */
311 integ->dobs[i] = diffvars->dobs[i];
312 #ifdef ANALYSE_DEBUG
313 varname = dis_make_name(integ->system,integ->dobs[i]);
314 CONSOLE_DEBUG("'%s' is observation",varname);
315 ASC_FREE(varname);
316 #endif
317 }
318
319 #ifdef ANALYSE_DEBUG
320 CONSOLE_DEBUG("rels matchbits: 0x%x",integrator_ida_rel.matchbits);
321 CONSOLE_DEBUG("rels matchvalue: 0x%x",integrator_ida_rel.matchvalue);
322
323 integrator_ida_rels_debug(integ);
324
325 CONSOLE_DEBUG("At the end of ida_analyse, there are %d rels active"
326 ,slv_count_solvers_rels(integ->system, &integrator_ida_rel)
327 );
328 #endif
329
330 return 0;
331 }
332
333
334 /*------------------------------------------------------------------------------
335 ADD ANY MISSING VARIABLES FOR WHICH DERIVATIVES ARE PRESENT
336 */
337
338 /**
339 This is the first step in the DAE analysis process. We inspect the
340 derivative chains from the solver's analyse() routine.
341 Check derivative chains: if a var or derivative is inelligible,
342 follow down the chain fixing & setting zero any derivatives.
343 */
344 static int integrator_ida_check_vars(IntegratorSystem *integ){
345 const SolverDiffVarCollection *diffvars;
346 char *varname;
347 int n_y = 0;
348 int i, j, k, n;
349 struct var_variable *v;
350 SolverDiffVarSequence seq;
351 int vok;
352 struct var_variable **vlist, **oldvl;
353 int var_found;
354
355 #ifdef ANALYSE_DEBUG
356 CONSOLE_DEBUG("BEFORE CHECKING VARS");
357 system_diffvars_debug(integ->system,stderr);
358 #endif
359
360 /* we shouldn't have allocated these yet: just be sure */
361 asc_assert(integ->y==NULL);
362 asc_assert(integ->ydot==NULL);
363 asc_assert(integ->n_y==0);
364
365 /* get the the dervative chains from the system */
366 diffvars = system_get_diffvars(integ->system);
367 if(diffvars==NULL){
368 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
369 return 1;
370 }
371
372 #ifdef ANALYSE_DEBUG
373 CONSOLE_DEBUG("AFTER system_get_diffvars...");
374 system_var_list_debug(integ->system);
375 #endif
376
377 /* add the variables from the derivative chains */
378 for(i=0; i<diffvars->nseqs; ++i){
379 seq = diffvars->seqs[i];
380 asc_assert(seq.n >= 1);
381 v = seq.vars[0];
382
383 /* update the var_fixed flag */
384 (void)var_fixed(v);
385 vok = var_apply_filter(v,&integrator_ida_nonderiv);
386 if(!var_incident(v) || !var_active(v)){
387 /* VARMSG("good var '%s' is not incident"); */
388 /* var meets our filter, but perhaps it's not incident? */
389 if(seq.n == 1 || var_apply_filter(seq.vars[1],&integrator_ida_nonderiv)){
390 /* VARMSG("DEACTIVATING NON-INCIDENT VAR '%s' (NO DERIVATIVE)"); */
391 var_set_active(v,0);
392 vok = 0;
393 }else{
394 var_found = 0;
395 #ifdef DERIV_WITHOUT_DIFF
396 VARMSG("'%s' has a derivative present, so needs to be included in the system");
397 #ifdef ANALYSE_DEBUG
398 CONSOLE_DEBUG("That var %s active",(var_active(v) ? "is" : "is NOT"));
399 #endif
400 var_set_active(v,1);
401 vok = var_apply_filter(v,&integrator_ida_nonderiv);
402 if (vok) {
403 if (!var_incident(v)) var_set_incident(v,1);
404 n = slv_get_num_solvers_vars(integ->system);
405 vlist = ASC_NEW_ARRAY(struct var_variable*,n + 2);
406 oldvl = slv_get_solvers_var_list(integ->system);
407 for(k = 0; k < n; k++) {
408 vlist[k] = oldvl[k];
409 if (vlist[k]==v) var_found = 1;
410 }
411 vlist[n] = v;
412 vlist[n+1] = NULL;
413 if (var_found) {
414 slv_set_solvers_var_list(integ->system,oldvl,n);
415 ascfree(vlist);
416 }
417 else {
418 slv_set_solvers_var_list(integ->system,vlist,n+1);
419 ascfree(oldvl);
420 }
421 }
422 #else
423 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Non-incident var with an incident derivative. ASCEND can't handle this case at the moment, but we hope to fix it.");
424 return 1;
425 #endif
426 }
427 }
428
429 if(!vok){
430 #ifdef ANALYSE_DEBUG
431 if(var_fixed(v)){
432 VARMSG("Fixed variable '%s' fails non-deriv filter");
433 }else{
434 VARMSG("'%s' fails non-deriv filter");
435 }
436 CONSOLE_DEBUG("passes nonderiv? %s (flags = 0x%x)"
437 , (var_apply_filter(v,&integrator_ida_nonderiv) ? "TRUE" : "false")
438 , var_flags(v)
439 );
440 #endif
441 for(j=1;j<seq.n;++j){
442 v = seq.vars[j];
443 var_set_active(v,FALSE);
444 var_set_value(v,0);
445 #ifdef ANALYSE_DEBUG
446 VARMSG("Derivative '%s' SET ZERO AND INACTIVE");
447 #endif
448 }
449 continue;
450 }
451
452 /* VARMSG("We will use var '%s'"); */
453 if(seq.n > 2){
454 varname = var_make_name(integ->system,v);
455 ERROR_REPORTER_HERE(ASC_USER_ERROR,"Too-long derivative chain with var '%s'",varname);
456 ASC_FREE(varname);
457 return 2;
458 }else if(seq.n==2){
459 /* diff var */
460 if(var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
461 /* add the diff & deriv vars to the lists */
462 n_y++;
463 /* VARMSG("Added diff var '%s'");
464 v = seq.vars[1]; VARMSG("and its derivative '%s'"); */
465 continue;
466 }
467 /* VARMSG("Diff var '%s' being converted to alg var...");
468 v = seq.vars[1]; VARMSG("...because deriv var '%s' fails filter"); */
469 /* fall through */
470 }
471
472 /* VARMSG("Adding '%s' to algebraic"); */
473 n_y++;
474 }
475
476 #ifdef ANALYSE_DEBUG
477 system_var_list_debug(integ->system);
478 #endif
479
480 /* we assert that all vars in y meet the integrator_ida_nonderiv filter */
481 /* we assert that all vars in ydot meet the integrator_ida_deriv filter */
482
483 #ifdef ANALYSE_DEBUG
484 CONSOLE_DEBUG("Found %d good non-derivative vars", n_y);
485 #endif
486 integ->n_y = n_y;
487
488 return 0;
489 }
490
491 /*------------------------------------------------------------------------------
492 FLAG THE RELATIONS CONTAINING DERIVATIVES
493 */
494 /**
495 Flag relations that contain derivatives as 'REL_DIFFERENTIAL'
496 in
497 @TODO what to do about relations that make reference to 't'?
498 */
499 static int integrator_ida_flag_rels(IntegratorSystem *integ){
500 int i, n, c, nd=0;
501 struct rel_relation **rels;
502 n = slv_get_num_solvers_rels(integ->system);
503 rels = slv_get_solvers_rel_list(integ->system);
504 for(i=0;i<n;++i){
505 c = rel_classify_differential(rels[i]);
506 if(c){
507 nd++;
508 /* CONSOLE_DEBUG("Rel %d is DIFFERENTIAL", i); */
509 }
510 }
511 #ifdef ANALYSE_DEBUG
512 CONSOLE_DEBUG("Found %d differential equations (so %d algebraic)",nd, n - nd);
513 #endif
514 integ->n_diffeqs = nd;
515 return 0;
516 }
517
518 /*------------------------------------------------------------------------------
519 SORT REL AND VAR LISTS
520 */
521 /**
522 Sort the lists. First we will put the non-derivative vars, then we will
523 put the derivative vars. Then we will put all the others.
524
525 See http://ascend4.org/File:Diffvars.png
526 */
527 static int integrator_ida_sort_rels_and_vars(IntegratorSystem *integ){
528 int ny1, nydot, nr;
529
530 #ifdef ANALYSE_DEBUG
531 CONSOLE_DEBUG("BEFORE SORTING RELS AND VARS");
532 system_diffvars_debug(integ->system,stderr);
533 #endif
534
535 /* we should not have allocated y or ydot yet */
536 asc_assert(integ->y==NULL && integ->ydot==NULL);
537
538 /* but we should have found some variables (and know how many) */
539 asc_assert(integ->n_y);
540
541 if(system_cut_vars(integ->system, 0, &integrator_ida_nonderiv, &ny1)){
542 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting non-derivs");
543 return 1;
544 }
545
546 #ifdef ANALYSE_DEBUG
547 CONSOLE_DEBUG("Cut %d non-derivative vars to start of list. cf integ->n_y = %d",ny1,integ->n_y);
548 integrator_ida_vars_debug(integ);
549 #endif
550 asc_assert(ny1 == integ->n_y);
551
552 ERROR_REPORTER_HERE(ASC_USER_NOTE,"moving derivs to start of remainder\n");
553 if(system_cut_vars(integ->system, ny1, &integrator_ida_deriv, &nydot)){
554 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
555 return 1;
556 }
557 #ifdef ANALYSE_DEBUG
558 CONSOLE_DEBUG("Cut %d derivatives to just after the non-derivatives",nydot);
559 #endif
560
561 if(system_cut_rels(integ->system, 0, &integrator_ida_rel, &nr)){
562 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem cutting derivs");
563 return 1;
564 }
565
566 #if 0
567 /* TODO we need to work out if it's reasonable to return an error now, or not... */
568 if(ny1 != nr){
569 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Problem is not square (ny = %d, nr = %d)",ny1,nr);
570 return 2;
571 }
572 #endif
573
574 return 0;
575 }
576
577 /*------------------------------------------------------------------------------
578 CREATE 'YDOT' AND 'Y_ID' LISTS
579 */
580 /**
581 Build up the lists 'ydot' and 'y_id.
582 'ydot' allows us to find the derivative of a variable given its sindex.
583 'y_id' allows us to locate the variable of a derivative given its sindex.
584
585 Hence
586 y_id[var_sindex(derivvar)-integ->n_y] = var_sindex(diffvar)
587 ydot[var_sindex(diffvar)] = derivvar (or NULL if diffvar is algebraic)
588
589 Note that there is 'shifting' required when looking up y_id.
590
591 Note that we want to get rid of 'y' but for the moment we are also assigning
592 that. We want to be rid of it because if the vars are ordered correctly
593 it shouldn't be needed.
594 */
595 static int integrator_ida_create_lists(IntegratorSystem *integ){
596 const SolverDiffVarCollection *diffvars;
597 int i, j;
598 struct var_variable *v;
599
600 SolverDiffVarSequence seq;
601
602 asc_assert(integ->y==NULL);
603 asc_assert(integ->ydot==NULL);
604 asc_assert(integ->y_id== NULL);
605
606 /* get the the dervative chains from the system */
607 diffvars = system_get_diffvars(integ->system);
608 if(diffvars==NULL){
609 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Derivative structure is empty");
610 return 1;
611 }
612
613 /* allocate space (more than enough) */
614 integ->y = ASC_NEW_ARRAY(struct var_variable *,integ->n_y);
615 integ->ydot = ASC_NEW_ARRAY_CLEAR(struct var_variable *,integ->n_y);
616 integ->n_ydot = 0;
617 for(i=0;i<integ->n_y;++i){
618 asc_assert(integ->ydot[i] == 0);
619 }
620
621 #ifdef ANALYSE_DEBUG
622 CONSOLE_DEBUG("Passing through chains...");
623 #endif
624 /* create the lists y and ydot, ignoring 'bad' vars */
625 for(i=0; i<diffvars->nseqs; ++i){
626 /* CONSOLE_DEBUG("i = %d",i); */
627
628 seq = diffvars->seqs[i];
629 asc_assert(seq.n >= 1);
630 v = seq.vars[0];
631 j = var_sindex(v);
632
633 if(!var_apply_filter(v,&integrator_ida_nonderiv)){
634 continue;
635 }
636
637 integ->y[j] = v;
638 /* VARMSG("'%s' is good non-deriv"); */
639
640 if(seq.n > 1 && var_apply_filter(seq.vars[1],&integrator_ida_deriv)){
641 v = seq.vars[1];
642 asc_assert(var_sindex(v) >= integ->n_y);
643 asc_assert(var_sindex(v)-integ->n_y < integ->n_y);
644
645 integ->ydot[j] = v;
646 integ->n_ydot++;
647 /* VARMSG("'%s' is good deriv"); */
648 }else{
649 asc_assert(integ->ydot[j]==NULL);
650 }
651 }
652
653 #ifdef ANALYSE_DEBUG
654 CONSOLE_DEBUG("Found %d good non-derivs",j);
655 #endif
656 /* create the list y_id by looking at non-NULLs from ydot */
657 integ->y_id = ASC_NEW_ARRAY(int,integ->n_ydot);
658 for(i=0,j=0; i < integ->n_y; ++i){
659 if(integ->ydot[i]==NULL)continue;
660 /* v = integ->ydot[i]; VARMSG("deriv '%s'..."); */
661 /* v = integ->y[i]; VARMSG("diff '%s'..."); */
662 integ->y_id[var_sindex(integ->ydot[i]) - integ->n_y] = i;
663 j++;
664 }
665
666 asc_assert(j==integ->n_ydot);
667
668 return 0;
669 }
670
671
672 /*------------------------------------------------------------------------------
673 CHECKING THE Y, YDOT AND Y_ID LISTS AGAINST THE SOLVERS_VAR_LISTS
674 */
675
676 /*
677 We are going to check that
678
679 - n_y in range (0,n_vars)
680 - no duplicates anywhere in the varlist
681 - no duplicates in non-NULL elements of ydot
682
683 - first n_y vars in solver's var list match those in vector y and match the non-deriv filter.
684 - var_sindex for first n_y vars match their position the the solver's var list
685
686 - ydot contains n_ydot non-NULL elements, each match the deriv filter.
687 - vars in solver's var list positions [n_y,n_y+n_ydot) all match deriv filter
688
689 - y_id contains n_ydot elements, with int values in range [0,n_y)
690 - var_sindex(ydot[y_id[i]]) - n_y == i for i in [0,n_ydot)
691
692 - all the vars [n_y+n_ydot,n_vars) fail the non-deriv filter and fail the deriv filter.
693 */
694 static int integrator_ida_check_diffindex(IntegratorSystem *integ){
695 int i, n_vars, n_ydot=0;
696 struct var_variable **list, *v;
697 char *varname;
698 const char *msg;
699
700 #ifdef ANALYSE_DEBUG
701 CONSOLE_DEBUG("Checking diffindex vector");
702 #endif
703
704 if(integ->y_id == NULL || integ->y == NULL || integ->ydot == NULL){
705 ERROR_REPORTER_HERE(ASC_PROG_ERR,"list(s) NULL");
706 return 1;
707 }
708
709 list = slv_get_solvers_var_list(integ->system);
710 n_vars = slv_get_num_solvers_vars(integ->system);
711
712 if(check_dups(integ, list, n_vars,FALSE)){
713 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates or NULLs in solver's var list");
714 return 1;
715 }
716
717 if(check_dups(integ, integ->ydot, integ->n_y,TRUE)){
718 ERROR_REPORTER_HERE(ASC_PROG_ERR,"duplicates in ydot vector");
719 return 1;
720 }
721
722 /* check n_y in range */
723 if(integ->n_y <= 0 || integ->n_y >= n_vars){
724 ERROR_REPORTER_HERE(ASC_PROG_ERR,"n_y = %d invalid (n_vars = %d)",integ->n_y, n_vars);
725 return 1;
726 }
727
728 /* check first n_y vars */
729 for(i=0; i < integ->n_y; ++i){
730 v = list[i];
731 if(!var_apply_filter(v, &integrator_ida_nonderiv)){
732 msg = "'%s' (in first n_y vars) fails non-deriv filter"; goto finish;
733 }else if(var_sindex(v) != i){
734 msg = "'%s' has wrong var_sindex"; goto finish;
735 }else if(v != integ->y[i]){
736 msg = "'%s' not matched in y vector"; goto finish;
737 }
738 /* meets filter, matches in y vector, has correct var_sindex. */
739 }
740
741 /* count non-NULLs in ydot */
742 for(i=0; i < integ->n_y; ++i){
743 v = integ->ydot[i];
744 if(v!=NULL){
745 if(var_sindex(v) < integ->n_y){
746 msg = "'%s' has var_sindex < n_y"; goto finish;
747 }else if(!var_apply_filter(v,&integrator_ida_deriv)){
748 msg = "'%s' (in next n_ydot vars) fails deriv filter"; goto finish;
749 }
750 /* lies beyond first n_y vars, match deriv filter */
751 n_ydot++;
752 }
753 }
754
755 /* check that vars [n_y, n_y+n_ydot) in solver's var list match the deriv filter */
756 for(i=integ->n_y; i< integ->n_y + n_ydot; ++i){
757 v = list[i];
758 if(!var_apply_filter(v,&integrator_ida_deriv)){
759 msg = "'%s', in range [n_y,n_y+n_ydot), fails deriv filter"; goto finish;
760 }
761 }
762
763 /* check values in y_id are ints int range [0,n_y), and that they point to correct vars */
764 for(i=0; i<n_ydot; ++i){
765 if(integ->y_id[i] < 0 || integ->y_id[i] >= integ->n_y){
766 ERROR_REPORTER_HERE(ASC_PROG_ERR,"value of y_id[%d] is out of range",i);
767 return 1;
768 }
769 v = integ->ydot[integ->y_id[i]];
770 if(var_sindex(v) - integ->n_y != i){
771 msg = "'%s' not index correctly from y_id"; goto finish;
772 }
773 }
774
775 /* check remaining vars fail both filters */
776 for(i=integ->n_y + n_ydot; i<n_vars; ++i){
777 v = list[i];
778 if(var_apply_filter(v,&integrator_ida_nonderiv)){
779 msg = "Var '%s' at end meets non-deriv filter, but shouldn't"; goto finish;
780 }
781 if(var_apply_filter(v,&integrator_ida_deriv)){
782 CONSOLE_DEBUG("position = %d",i);
783 msg = "Var '%s' at end meets deriv filter, but shouldn't"; goto finish;
784 }
785 }
786
787 return 0;
788 finish:
789 varname = var_make_name(integ->system,v);
790 ERROR_REPORTER_HERE(ASC_PROG_ERR,msg,varname);
791 ASC_FREE(varname);
792 return 1;
793 }
794
795 /* return 0 on succes */
796 static int check_dups(IntegratorSystem *integ, struct var_variable **list,int n,int allownull){
797 int i,j;
798 struct var_variable *v;
799 #ifdef ANALYSE_DEBUG
800 char *varname;
801 #endif
802 for(i=0; i< n; ++i){
803 v=list[i];
804 if(v==NULL){
805 if(allownull)continue;
806 else return 2;
807 }
808 asc_assert(v!=(void *)0x31);
809 for(j=0; j<i-1;++j){
810 if(list[j]==NULL)continue;
811 if(v==list[j]){
812 #ifdef ANALYSE_DEBUG
813 varname = var_make_name(integ->system,v);
814 if(varname){
815 CONSOLE_DEBUG("Duplicate of '%s' found",varname);
816 ASC_FREE(varname);
817 }else{
818 CONSOLE_DEBUG("Duplicate found (couldn't retrieve name)");
819 }
820 ASC_FREE(varname);
821 #endif
822 return 1;
823 }
824 }
825 }
826 return 0;
827 }
828
829
830 /*------------------------------------------------------------------------------
831 INDEX CHECKING
832 */
833 /**
834 Check for index-1 DAE system. We plan to do this by checking that df/dyd'
835 and dg/dya are both non-singular (and of course square). Valid systems
836 can (we think?) be written that don't meet these requirements but they
837 will have index problems and may not solve with IDA.
838 */
839 static int integrator_ida_check_index(IntegratorSystem *integ){
840 #if 1
841 linsolqr_system_t L;
842 mtx_range_t range;
843 mtx_region_t R;
844 int res, r;
845 struct SystemJacobianStruct df_dydp, dg_dya;
846
847 #ifdef ANALYSE_DEBUG
848 CONSOLE_DEBUG("system has total of %d rels and %d vars"
849 ,slv_get_num_solvers_rels(integ->system)
850 ,slv_get_num_solvers_vars(integ->system)
851 );
852
853 CONSOLE_DEBUG("VAR_DERIV = 0x%x = %d",VAR_DERIV, VAR_DERIV);
854 CONSOLE_DEBUG("system_vfilter_deriv.matchbits = 0x%x",system_vfilter_deriv.matchbits);
855 CONSOLE_DEBUG("system_vfilter_deriv.matchvalue= 0x%x",system_vfilter_deriv.matchvalue);
856 #endif
857
858 asc_assert(system_vfilter_deriv.matchbits & VAR_DERIV);
859 asc_assert(system_vfilter_deriv.matchvalue & VAR_DERIV);
860
861 #ifdef ANALYSE_DEBUG
862 CONSOLE_DEBUG("system has %d vars matching deriv filter",slv_count_solvers_vars(integ->system, &system_vfilter_deriv));
863 #endif
864
865 res = system_jacobian(integ->system
866 , &system_rfilter_diff
867 , &system_vfilter_deriv
868 , 1 /* safe evaluation */
869 , &df_dydp
870 );
871
872 if(res){
873 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating df/dyd'");
874 }
875
876 #ifdef ANALYSE_DEBUG
877 CONSOLE_DEBUG("df/dyd': nr = %d, nv = %d",df_dydp.n_rels,df_dydp.n_vars);
878 #endif
879
880 res = system_jacobian(integ->system
881 , &system_rfilter_algeb
882 , &system_vfilter_algeb
883 , 1 /* safe evaluation */
884 , &dg_dya
885 );
886
887 if(res){
888 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error calculating dg/dya");
889 }
890
891 #ifdef ANALYSE_DEBUG
892 CONSOLE_DEBUG("dg/dya: nr = %d, nv = %d",dg_dya.n_rels,dg_dya.n_vars);
893 #endif
894
895 if((df_dydp.n_rels == 0) ^ (df_dydp.n_vars == 0)){
896 ERROR_REPORTER_HERE(ASC_PROG_ERR,"df/dyd' is a bit ambiguous");
897 }
898
899 if(dg_dya.n_rels <= 0){
900 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No algebraic equations were found in the DAE system!");
901 }else if(dg_dya.n_rels != dg_dya.n_vars){
902 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The algebraic part of the DAE jacobian, dg/dya, is not square!");
903 }else{
904 /* check the rank */
905 range.low = 0; range.high = mtx_order(dg_dya.M) - 1;
906 R.row = range; R.col = range;
907
908 L = linsolqr_create_default();
909 linsolqr_set_matrix(L,dg_dya.M);
910 linsolqr_set_region(L,R);
911 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
912 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
913 linsolqr_factor(L,linsolqr_fmethod(L));
914 r = linsolqr_rank(L);
915
916 linsolqr_destroy(L);
917
918 if(r != dg_dya.n_rels){
919 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix dg/dya is not full rank");
920 }
921 }
922
923 ASC_FREE(dg_dya.vars);
924 ASC_FREE(dg_dya.rels);
925 mtx_destroy(dg_dya.M);
926
927 if(df_dydp.n_rels <= 0){
928 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"No differential equations were found in the DAE system!");
929 }else if(df_dydp.n_rels != df_dydp.n_vars){
930 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"The differential part of the the jacobian dg/dya is not square!");
931 ASC_FREE(df_dydp.vars);
932 ASC_FREE(df_dydp.rels);
933 mtx_destroy(df_dydp.M);
934 return 1;
935 }else{
936 /* check the rank */
937 range.low = 0; range.high = mtx_order(df_dydp.M) - 1;
938 R.row = range; R.col = range;
939
940 L = linsolqr_create_default();
941 linsolqr_set_matrix(L,df_dydp.M);
942 linsolqr_set_region(L,R);
943 linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
944 linsolqr_reorder(L, &R, linsolqr_rmethod(L));
945 linsolqr_factor(L,linsolqr_fmethod(L));
946 r = linsolqr_rank(L);
947
948 linsolqr_destroy(L);
949
950 if(r != df_dydp.n_rels){
951 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Your DAE system has an index problem: the matrix df/dyd' is not full rank");
952 }
953 }
954
955 if(df_dydp.n_rels + dg_dya.n_rels == 0){
956 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Both df/dyd' and dg/dya were empty!");
957 }
958
959 ASC_FREE(df_dydp.vars);
960 ASC_FREE(df_dydp.rels);
961 mtx_destroy(df_dydp.M);
962 return 0;
963 #else
964 ERROR_REPORTER_HERE(ASC_PROG_ERR,"check_index disabled");
965 return 0;
966 #endif
967 }
968
969 /*------------------------------------------------------------------------------
970 CHECKING FOR STRUCTURAL SINGULARITY
971 */
972
973 #if 0 /* disused function, apparently -- JP */
974 static int integrator_ida_check_partitioning(IntegratorSystem *integ){
975 long i, nv;
976 int err=0;
977 char *varname;
978 struct var_variable **vlist, *v;
979 nv = slv_get_num_solvers_vars(integ->system);
980 vlist = slv_get_solvers_var_list(integ->system);
981 var_filter_t vf = {VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|VAR_FIXED
982 ,VAR_SVAR|VAR_ACTIVE|VAR_INCIDENT|0 };
983 for(i=0;i<nv;++i){
984 v=vlist[i];
985 if(!var_apply_filter(v,&vf))continue;
986 varname = var_make_name(integ->system,v);
987 if(!var_deriv(v)){
988 #ifdef ANALYSE_DEBUG
989 fprintf(stderr,"vlist[%ld] = '%s' (nonderiv)\n",i,varname);
990 #endif
991 if(i>=integ->n_y){
992 ERROR_REPORTER_HERE(ASC_PROG_ERR,"non-deriv var '%s' is not at the start",varname);
993 err++;
994 }
995 }else{
996 #ifdef ANALYSE_DEBUG
997 fprintf(stderr,"vlist[%ld] = '%s' (derivative)\n",i,varname);
998 #endif
999 if(i<integ->n_y){
1000 ERROR_REPORTER_HERE(ASC_PROG_ERR,"deriv var '%s' is not at the end (n_y = %d, i = %d)"
1001 ,varname, integ->n_y, i
1002 );
1003 err++;
1004 }
1005 }
1006 ASC_FREE(varname);
1007 }
1008 if(err){
1009 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error in IDA partitioning");
1010 return 1;
1011 }
1012 return 0;
1013 }
1014
1015 int integrator_ida_block_check(IntegratorSystem *integ){
1016 int res;
1017 int dof;
1018 #ifdef ANALYSE_DEBUG
1019 int *vlist, *vp, i, nv, nv_ok;
1020 char *varname;
1021 struct var_variable **solversvars;
1022
1023 var_filter_t vfilt = {
1024 VAR_ACTIVE | VAR_INCIDENT | VAR_FIXED,
1025 VAR_ACTIVE | VAR_INCIDENT | 0
1026 };
1027
1028 nv = slv_get_num_solvers_vars(integ->system);
1029 solversvars = slv_get_solvers_var_list(integ->system);
1030 CONSOLE_DEBUG("-------------- nv = %d -------------",nv);
1031 for(nv_ok=0, i=0; i < nv; ++i){
1032 if(var_apply_filter(solversvars[i],&vfilt)){
1033 varname = var_make_name(integ->system,solversvars[i]);
1034 fprintf(stderr,"%s\n",varname);
1035 ASC_FREE(varname);
1036 nv_ok++;
1037 }
1038 }
1039 CONSOLE_DEBUG("----------- got %d ok -------------",nv_ok);
1040 #endif
1041
1042 if(!slvDOF_status(integ->system, &res, &dof)){
1043 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to determine DOF status");
1044 return -1;
1045 }
1046 switch(res){
1047 case 1: CONSOLE_DEBUG("System is underspecified (%d degrees of freedom)",dof);break;
1048 case 2: CONSOLE_DEBUG("System is square"); return 0; /* all OK */
1049 case 3: CONSOLE_DEBUG("System is structurally singular"); break;
1050 case 4: CONSOLE_DEBUG("System is overspecified"); break;
1051 default:
1052 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unrecognised slfDOF_status");
1053 return -2;
1054 }
1055
1056 #ifdef ANALYSE_DEBUG
1057 /* if it was underspecified, what vars could be fixed? */
1058 if(res==1){
1059 CONSOLE_DEBUG("Need to FIX %d of the following vars:",dof);
1060 solversvars = slv_get_solvers_var_list(integ->system);
1061 if(!slvDOF_eligible(integ->system, &vlist)){
1062 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unable to det slvDOF_eligble list");
1063 return -3;
1064 }
1065 for(vp=vlist;*vp!=-1;++vp){
1066 varname = var_make_name(integ->system, solversvars[*vp]);
1067 CONSOLE_DEBUG("Fixable var: %s",varname);
1068 ASC_FREE(varname);
1069 }
1070 CONSOLE_DEBUG("(Found %d fixable vars)",(int)(vp-vlist));
1071 return 1;
1072 }
1073 #endif
1074
1075 return res;
1076 }
1077 #endif /* disused code */
1078
1079 /*------------------------------------------------------------------------------
1080 DEBUGGING/OUTPUT ROUTINES (should be moved to idaio.c)
1081 */
1082
1083 /**
1084 This function prints out a table of the solvers_var_list with the relevant
1085 variable tags also shown. It is aimed to allow debugging of problems with
1086 the IDA analysis routines, but perhaps it could be useful enough to be
1087 reported back to the user?
1088 */
1089 static int integrator_ida_vars_debug(IntegratorSystem *integ){
1090 struct var_variable **list;
1091 int n, i;
1092 /* get all the var names, get the length of the longest one */
1093 n = slv_get_num_solvers_vars(integ->system);
1094 list = slv_get_solvers_var_list(integ->system);
1095 char *name[n]; // we'll need to free all these pointers later
1096 int lmax = 0;
1097 for(i=0;i<n;++i){
1098 name[i] = var_make_name(integ->system,list[i]);
1099 int l = strlen(name[i]);
1100 if(l>lmax)lmax=l;
1101 }
1102
1103 #define FSTRING " %-10s %-*s %6s %6s %6s %6s %6s"
1104 #define AST(F) ((var_##F(v))?"*":"")
1105 fprintf(stderr,"%6s" FSTRING "\n","sindex","",lmax,"","incid","active","fixed","deriv","indep");
1106 fprintf(stderr,"%6s" FSTRING "\n","------","",lmax,"","-----","------","-----","-----","-----");
1107 struct var_variable *v;
1108 for(i=0;i<integ->n_y;++i){
1109 v = list[i];
1110 fprintf(stderr,"%6d" FSTRING "\n",i,(integ->ydot?(integ->ydot[i]?"diff":"algeb"):"?"),lmax,name[i]
1111 ,AST(incident),AST(active),AST(fixed),AST(deriv),AST(nonbasic)
1112 );
1113 }
1114 fprintf(stderr,"%6s %-10s\n","","-------");
1115 for(i=integ->n_y;i<integ->n_y+integ->n_ydot;++i){
1116 v = list[i];
1117 fprintf(stderr,"%6d" FSTRING " (of %s)\n",i,"deriv",lmax,name[i]
1118 ,AST(incident),AST(active),AST(fixed),AST(deriv),AST(nonbasic)
1119 ,integ->y_id?(name[integ->y_id[i - integ->n_y]]):"?"
1120 );
1121 }
1122 fprintf(stderr,"%6s %-10s\n","","-------");
1123 for(i=integ->n_y+integ->n_ydot;i<n;++i){
1124 v = list[i];
1125 fprintf(stderr,"%6d" FSTRING "\n",i,"other",lmax,name[i]
1126 ,AST(incident),AST(active),AST(fixed),AST(deriv),AST(nonbasic)
1127 );
1128 }
1129 fprintf(stderr,"%6s %-10s\n","","=======");
1130 #undef FSTRING
1131 #undef AST
1132 for(i=0;i<n;++i){
1133 ASC_FREE(name[i]);
1134 }
1135
1136 return 0;
1137 }
1138
1139
1140 static int integrator_ida_rels_debug(IntegratorSystem *integ){
1141 struct rel_relation **list;
1142 int n, i;
1143 /* get all the rel names, get the length of the longest one */
1144 n = slv_get_num_solvers_rels(integ->system);
1145 list = slv_get_solvers_rel_list(integ->system);
1146 char *name[n]; // we'll need to free all these pointers later
1147 int lmax = 0;
1148 for(i=0;i<n;++i){
1149 name[i] = rel_make_name(integ->system,list[i]);
1150 int l = strlen(name[i]);
1151 if(l>lmax)lmax=l;
1152 }
1153
1154 #define FSTRING " %-10s %-*s %4s %6s %4s %5s %4s %4s %4s"
1155 #define AST(F) ((rel_##F(r))?"*":"")
1156 fprintf(stderr,"%6s" FSTRING "\n","sindex","",lmax,"","incl","active","diff","cond","event","when","bbox");
1157 fprintf(stderr,"%6s" FSTRING "\n","------","",lmax,"","----","------","----","----","-----","----","----");
1158 struct rel_relation *r;
1159 for(i=0;i<n;++i){
1160 r = list[i];
1161 fprintf(stderr,"%6d" FSTRING "\n",i,(rel_differential(r)?"diff":"algeb")
1162 ,lmax,name[i]
1163 ,AST(included),AST(active),AST(differential),AST(conditional),AST(in_when),AST(in_event),AST(blackbox)
1164 );
1165 }
1166 fprintf(stderr,"%6s %-10s\n","","=======");
1167 #undef FSTRING
1168 #undef AST
1169 for(i=0;i<n;++i){
1170 ASC_FREE(name[i]);
1171 }
1172 }
1173
1174 /*------------------------------------------------------------------------------
1175 UTILITY FUNCTIONS (used elsewhere, not only in this file)
1176 */
1177
1178 /**
1179 Given a derivative variable, return the index of its corresponding differential
1180 variable in the y vector (and equivalently the var_sindex of the diff var)
1181 */
1182 int integrator_ida_diffindex(const IntegratorSystem *integ, const struct var_variable *deriv){
1183 asc_assert(var_sindex(deriv) >= integ->n_y);
1184 asc_assert(var_sindex(deriv) < integ->n_y + integ->n_ydot);
1185 return integ->y_id[var_sindex(deriv) - integ->n_y];
1186 }
1187
1188 /**
1189 A less assertive version of diffindex...
1190 */
1191 int integrator_ida_diffindex1(const IntegratorSystem *integ, const struct var_variable *deriv){
1192 if(var_sindex(deriv) >= integ->n_y)return -1;
1193 if(var_sindex(deriv) < integ->n_y + integ->n_ydot)return -2;
1194 return integ->y_id[var_sindex(deriv) - integ->n_y];
1195 }
1196
1197

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22