/[ascend]/trunk/pygtk/simulation.cpp
ViewVC logotype

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 921 - (show annotations) (download) (as text)
Mon Nov 6 07:49:06 2006 UTC (15 years, 7 months ago) by johnpye
File MIME type: text/x-c++src
File size: 22582 byte(s)
Fixing up some malloc/free problems with  integrator.c, removing some debug output from slv3.
Working on fixing IDA for systems with inactive variables (ongoing).
1 /* ASCEND modelling environment
2 Copyright (C) 2006 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, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 */
19 #include <iostream>
20 #include <iomanip>
21 #include <stdexcept>
22 #include <sstream>
23 using namespace std;
24
25 #include "config.h"
26
27 extern "C"{
28 #include <utilities/ascConfig.h>
29 #include <utilities/error.h>
30 #include <utilities/ascSignal.h>
31 #include <utilities/ascMalloc.h>
32 #include <general/dstring.h>
33 #include <general/tm_time.h>
34 #include <compiler/instance_enum.h>
35 #include <compiler/fractions.h>
36 #include <compiler/compiler.h>
37 #include <compiler/dimen.h>
38 #include <compiler/symtab.h>
39 #include <compiler/instance_io.h>
40 #include <compiler/instantiate.h>
41 #include <compiler/bintoken.h>
42 #include <compiler/instance_enum.h>
43 #include <compiler/instquery.h>
44 #include <compiler/check.h>
45 #include <compiler/name.h>
46 #include <compiler/pending.h>
47 #include <compiler/importhandler.h>
48
49 #include <utilities/readln.h>
50 #include <solver/mtx.h>
51 #include <solver/slv_types.h>
52 #include <solver/var.h>
53 #include <solver/rel.h>
54 #include <solver/discrete.h>
55 #include <solver/conditional.h>
56 #include <solver/logrel.h>
57 #include <solver/bnd.h>
58 #include <solver/calc.h>
59 #include <solver/relman.h>
60 #include <solver/slv_common.h>
61 #include <solver/linsol.h>
62 #include <solver/linsolqr.h>
63 #include <solver/slv_client.h>
64 #include <solver/system.h>
65 #include <solver/slv_interface.h>
66 #include <solver/slvDOF.h>
67 #include <solver/slv3.h>
68 #include <solver/slv_stdcalls.h>
69 #include <solver/slv_server.h>
70 }
71
72 #include "simulation.h"
73 #include "solver.h"
74 #include "solverparameters.h"
75 #include "name.h"
76 #include "incidencematrix.h"
77 #include "variable.h"
78 #include "solverstatus.h"
79 #include "solverreporter.h"
80
81 /**
82 Create an instance of a type (call compiler etc)
83
84 @TODO fix mutex on compile command filenames
85 */
86 Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
87 sys = NULL;
88 //is_built = false;
89 // Create an Instance object for the 'simulation root' (we'll call
90 // it the 'simulation model') and it can be fetched using 'getModel()'
91 // any time later.
92 //simroot = Instanc(GetSimulationRoot(i),name);
93 }
94
95 Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
96 //is_built = old.is_built;
97 sys = old.sys;
98 bin_srcname = old.bin_srcname;
99 bin_objname = old.bin_objname;
100 bin_libname = old.bin_libname;
101 bin_cmd = old.bin_cmd;
102 bin_rm = old.bin_rm;
103 sing = NULL;
104 }
105
106 Simulation::~Simulation(){
107 //CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
108 }
109
110 Instanc &
111 Simulation::getModel(){
112 if(!simroot.getInternalType()){
113 throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
114 }
115 if(InstanceKind(simroot.getInternalType())!=MODEL_INST){
116 throw runtime_error("Simulation::getModel: simroot is not a MODEL instance");
117 }
118 return simroot;
119 }
120
121
122 slv_system_structure *
123 Simulation::getSystem(){
124 if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
125 return sys;
126 }
127
128
129 const string
130 Simulation::getInstanceName(const Instanc &i) const{
131 char *n;
132 n = WriteInstanceNameString(i.getInternalType(),simroot.getInternalType());
133 string s(n);
134 ascfree(n);
135 return s;
136 }
137
138 const int
139 Simulation::getNumVars(){
140 return slv_get_num_solvers_vars(getSystem());
141 }
142
143
144 void
145 Simulation::write(){
146 simroot.write();
147 }
148
149 //------------------------------------------------------------------------------
150 // RUNNING MODEL 'METHODS'
151
152 void
153 Simulation::run(const Method &method){
154 Instanc &model = getModel();
155 this->run(method,model);
156 }
157
158 void
159 Simulation::runDefaultMethod(){
160 Method m = getType().getMethod(SymChar("on_load"));
161 run(m);
162 }
163
164 void
165 Simulation::run(const Method &method, Instanc &model){
166
167 // set the 'sim' pointer to our local variable...
168 CONSOLE_DEBUG("Setting shared pointer 'sim'");
169 importhandler_setsharedpointer("sim",this);
170
171 /*if(not is_built){
172 CONSOLE_DEBUG("WARNING, SIMULATION NOT YET BUILT");
173 }*/
174
175 CONSOLE_DEBUG("Running method %s...", method.getName());
176
177 Nam name = Nam(method.getSym());
178 //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
179
180
181 Proc_enum pe;
182 pe = Initialize(
183 &*(model.getInternalType()) ,name.getInternalType(), "__not_named__"
184 ,ASCERR
185 ,0, NULL, NULL
186 );
187
188 // clear out the 'sim' pointer (soon it will be invalid)
189 importhandler_setsharedpointer("sim",NULL);
190 CONSOLE_DEBUG("Cleared shared pointer 'sim'");
191
192 if(pe == Proc_all_ok){
193 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Method '%s' was run (check above for errors)\n",method.getName());
194 //cerr << "METHOD " << method.getName() << " COMPLETED OK" << endl;
195 }else{
196 stringstream ss;
197 ss << "Simulation::run: Method '" << method.getName() << "' returned error: ";
198 switch(pe){
199 case Proc_CallOK: ss << "Call OK"; break;
200 case Proc_CallError: ss << "Error occurred in call"; break;
201 case Proc_CallReturn: ss << "Request that caller return (OK)"; break;
202 case Proc_CallBreak: ss << "Break out of enclosing loop"; break;
203 case Proc_CallContinue: ss << "Skip to next iteration"; break;
204
205 case Proc_break: ss << "Break"; break;
206 case Proc_continue: ss << "Continue"; break;
207 case Proc_fallthru: ss << "Fall-through"; break;
208 case Proc_return: ss << "Return"; break;
209 case Proc_stop: ss << "Stop"; break;
210 case Proc_stack_exceeded: ss << "Stack exceeded"; break;
211 case Proc_stack_exceeded_this_frame: ss << "Stack exceeded this frame"; break;
212 case Proc_case_matched: ss << "Case matched"; break;
213 case Proc_case_unmatched: ss << "Case unmatched"; break;
214
215 case Proc_case_undefined_value: ss << "Undefined value in case"; break;
216 case Proc_case_boolean_mismatch: ss << "Boolean mismatch in case"; break;
217 case Proc_case_integer_mismatch: ss << "Integer mismatch in case"; break;
218 case Proc_case_symbol_mismatch: ss << "Symbol mismatch in case"; break;
219 case Proc_case_wrong_index: ss << "Wrong index in case"; break;
220 case Proc_case_wrong_value: ss << "Wrong value in case"; break;
221 case Proc_case_extra_values: ss << "Extra values in case"; break;
222 case Proc_bad_statement: ss << "Bad statement"; break;
223 case Proc_bad_name: ss << "Bad name"; break;
224 case Proc_for_duplicate_index: ss << "Duplicate index"; break;
225 case Proc_for_set_err: ss << "For set error"; break;
226 case Proc_for_not_set: ss << "For not set"; break;
227 case Proc_illegal_name_use: ss << "Illegal name use"; break;
228 case Proc_name_not_found: ss << "Name not found"; break;
229 case Proc_instance_not_found: ss << "Instance not found"; break;
230 case Proc_type_not_found: ss << "Type not found"; break;
231 case Proc_illegal_type_use: ss << "Illegal use"; break;
232 case Proc_proc_not_found: ss << "Method not found"; break;
233 case Proc_if_expr_error_typeconflict: ss << "Type conflict in 'if' expression"; break;
234 case Proc_if_expr_error_nameunfound: ss << "Name not found in 'if' expression"; break;
235 case Proc_if_expr_error_incorrectname: ss << "Incorrect name in 'if' expression"; break;
236 case Proc_if_expr_error_undefinedvalue: ss << "Undefined value in 'if' expression"; break;
237 case Proc_if_expr_error_dimensionconflict: ss << "Dimension conflict in 'if' expression"; break;
238 case Proc_if_expr_error_emptychoice: ss << "Empty choice in 'if' expression"; break;
239 case Proc_if_expr_error_emptyintersection: ss << "Empty intersection in 'if' expression"; break;
240 case Proc_if_expr_error_confused: ss << "Confused in 'if' expression"; break;
241 case Proc_if_real_expr: ss << "Real-valued result in 'if' expression"; break;
242 case Proc_if_integer_expr: ss << "Integeter-valued result in 'if' expression"; break;
243 case Proc_if_symbol_expr: ss << "Symbol-valued result in 'if' expression"; break;
244 case Proc_if_set_expr: ss << "Set-valued result in 'if' expression"; break;
245 case Proc_if_not_logical: ss << "If expression is not logical"; break;
246 case Proc_user_interrupt: ss << "User interrupt"; break;
247 case Proc_infinite_loop: ss << "Infinite loop"; break;
248 case Proc_declarative_constant_assignment: ss << "Declarative constant assignment"; break;
249 case Proc_nonsense_assignment: ss << "Nonsense assginment (bogus)"; break;
250 case Proc_nonconsistent_assignment: ss << "Inconsistent assignment"; break;
251 case Proc_nonatom_assignment: ss << "Non-atom assignment"; break;
252 case Proc_nonboolean_assignment: ss << "Non-boolean assignment"; break;
253 case Proc_noninteger_assignment: ss << "Non-integer assignment"; break;
254 case Proc_nonreal_assignment: ss << "Non-real assignment"; break;
255 case Proc_nonsymbol_assignment: ss << "Non-symbol assignment"; break;
256 case Proc_lhs_error: ss << "Left-hand-side error"; break;
257 case Proc_rhs_error: ss << "Right-hand-side error"; break;
258 case Proc_unknown_error: ss << "Unknown error"; break;
259 default:
260 ss << "Invalid error code";
261 }
262
263 ss << " (" << int(pe) << ")";
264 throw runtime_error(ss.str());
265 }
266 }
267
268 //-----------------------------------------------------------------------------
269 // CHECKING METHODS
270
271 /**
272 Check that all the analysis went OK: solver lists are all there, etc...?
273
274 Can't return anything here because of limitations in the C API
275
276 @TODO there's something wrong with this at the moment: even after 'FIX'
277 methods are run, check shows them as not fixed, up until the point that 'SOLVE'
278 successfully completes. Something's not being synchronised properly...
279 */
280 void
281 Simulation::checkInstance(){
282 //cerr << "CHECKING SIMULATION INSTANCE" << endl;
283 /*if(!is_built){
284 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Simulation has not been built");
285 return;
286 }*/
287 Instance *i1 = getModel().getInternalType();
288 CheckInstance(stderr, &*i1);
289 //cerr << "DONE CHECKING INSTANCE" << endl;
290 }
291
292 /**
293 @return 1 = underspecified, 2 = square, 3 = structurally singular, 4 = overspecified
294 */
295 enum StructuralStatus
296 Simulation::checkDoF() const{
297 int dof, status;
298
299 if(!sys){
300 throw runtime_error("System is not built");
301 }
302
303 /*if(!is_built){
304 throw runtime_error("System not yet built");
305 }*/
306 CONSOLE_DEBUG("Calling slvDOF_status...");
307 slvDOF_status(sys, &status, &dof);
308 switch(status){
309 case ASCXX_DOF_UNDERSPECIFIED:
310 case ASCXX_DOF_SQUARE:
311 case ASCXX_DOF_OVERSPECIFIED:
312 case ASCXX_DOF_STRUCT_SINGULAR:
313 return (enum StructuralStatus)status;
314 case 5:
315 throw runtime_error("Unable to resolve degrees of freedom"); break;
316 default:
317 throw runtime_error("Invalid return status from slvDOF_status");
318 }
319 }
320
321 /**
322 Check consistency
323
324 @TODO what is the difference between this and checkStructuralSingularity?
325
326 @return list of freeable variables. List will be empty if sys is consistent.
327 */
328 vector<Variable>
329 Simulation::getFreeableVariables(){
330 vector<Variable> v;
331
332 //cerr << "CHECKING CONSISTENCY..." << endl;
333 int *fixedarrayptr=NULL;
334
335 if(!sys){
336 throw runtime_error("System not yet built");
337 }
338
339 int res = consistency_analysis(sys, &fixedarrayptr);
340
341 if(res==1){
342 cerr << "STRUCTURALLY CONSISTENT" << endl;
343 }else{
344 if(fixedarrayptr ==NULL){
345 ERROR_REPORTER_HERE(ASC_USER_ERROR,"STRUCTURALLY INCONSISTENT");
346 throw runtime_error("Invalid consistency analysis result returned!");
347 }
348
349 struct var_variable **vp = slv_get_master_var_list(sys);
350 for(int i=0; fixedarrayptr[i]!=-1; ++i){
351 v.push_back( Variable(this, vp[fixedarrayptr[i]]) );
352 }
353 }
354 return v;
355 }
356
357 /** Returns TRUE if all is OK (not singular) */
358 bool
359 Simulation::checkStructuralSingularity(){
360 int *vil;
361 int *ril;
362 int *fil;
363
364 cerr << "RETRIEVING slfDOF_structsing INFO" << endl;
365
366 int res = slvDOF_structsing(sys, mtx_FIRST, &vil, &ril, &fil);
367 struct var_variable **varlist = slv_get_solvers_var_list(sys);
368 struct rel_relation **rellist = slv_get_solvers_rel_list(sys);
369
370 if(this->sing){
371 cerr << "DELETING OLD SINGULATING INFO" << endl;
372 delete this->sing;
373 this->sing = NULL;
374 }
375
376 if(res==1){
377 CONSOLE_DEBUG("processing singularity data...");
378 sing = new SingularityInfo();
379
380 // pull in the lists of vars and rels, and the freeable vars:
381 for(int i=0; ril[i]!=-1; ++i){
382 sing->rels.push_back( Relation(this, rellist[ril[i]]) );
383 }
384
385 for(int i=0; vil[i]!=-1; ++i){
386 sing->vars.push_back( Variable(this, varlist[vil[i]]) );
387 }
388
389 for(int i=0; fil[i]!=-1; ++i){
390 sing->freeablevars.push_back( Variable(this, varlist[fil[i]]) );
391 }
392
393 // we're done with those lists now
394 ASC_FREE(vil);
395 ASC_FREE(ril);
396 ASC_FREE(fil);
397
398 if(sing->isSingular()){
399 CONSOLE_DEBUG("singularity found");
400 this->sing = sing;
401 return FALSE;
402 }
403 CONSOLE_DEBUG("no singularity");
404 delete sing;
405 return TRUE;
406 }else{
407 if(res==0){
408 throw runtime_error("Unable to determine singularity lists");
409 }else{
410 throw runtime_error("Invalid return from slvDOF_structsing.");
411 }
412 }
413 }
414
415 /**
416 If the checkStructuralSingularity analysis has been done,
417 this funciton will let you access the SingularityInfo data that was
418 stored.
419 */
420 const SingularityInfo &
421 Simulation::getSingularityInfo() const{
422 if(sing==NULL){
423 throw runtime_error("No singularity info present");
424 }
425 return *sing;
426 }
427
428 //------------------------------------------
429 // ASSIGNING SOLVER TO SIMULATION
430
431 void
432 Simulation::setSolver(Solver &solver){
433 /* CONSOLE_DEBUG("Setting solver on sim %p, root inst %p",this,this->simroot.getInternalType()); */
434
435 try{
436 build();
437 }catch(runtime_error &e){
438 stringstream ss;
439 ss << "Couldn't prepare system for solving:";
440 ss << e.what();
441 throw runtime_error(ss.str());
442 }
443
444 CONSOLE_DEBUG("Calling slv_select_solver...");
445 int selected = slv_select_solver(sys, solver.getIndex());
446 //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
447
448 if(selected<0){
449 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
450 throw runtime_error("Failed to select solver");
451 }
452
453 if(selected!=solver.getIndex()){
454 solver = Solver(slv_solver_name(selected));
455 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
456 }
457
458 if( slv_eligible_solver(sys) <= 0){
459 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
460 throw runtime_error("Inelegible solver");
461 }
462 }
463
464 const Solver
465 Simulation::getSolver() const{
466 int index = slv_get_selected_solver(sys);
467 //cerr << "Simulation::getSolver: index = " << index << endl;
468 if(index<0)throw runtime_error("No solver selected");
469
470 return Solver(slv_solver_name(index));
471 }
472
473 //------------------------------------------------------------------------------
474 // BUILD THE SYSTEM
475
476 /**
477 Build the system (send it to the solver)
478 */
479 void
480 Simulation::build(){
481 if(sys){
482 CONSOLE_DEBUG("System is already built");
483 return;
484 }
485
486 if(simroot.getKind() != MODEL_INST){
487 throw runtime_error("Simulation does not contain a MODEL_INST");
488 }
489
490 if(NumberPendingInstances(simroot.getInternalType())){
491 throw runtime_error("System has pending instances; can't yet send to solver.");
492 }
493
494 sys = system_build(simroot.getInternalType());
495 if(!sys){
496 throw runtime_error("Unable to build system");
497 }
498
499 CONSOLE_DEBUG("System built OK");
500 }
501
502
503 //------------------------------------------------------------------------------
504 // SOLVER CONFIGURATION PARAMETERS
505
506 /**
507 Get solver parameters struct wrapped up as a SolverParameters class.
508 */
509 SolverParameters
510 Simulation::getSolverParameters() const{
511 //if(!is_built)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
512 if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has no 'sys' assigned.");
513
514 slv_parameters_t p;
515 slv_get_parameters(sys,&p);
516 return SolverParameters(p);
517 }
518
519 /**
520 Update the solver parameters by passing a new set back
521 */
522 void
523 Simulation::setSolverParameters(SolverParameters &P){
524 if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
525 CONSOLE_DEBUG("Calling slv_set_parameters");
526 slv_set_parameters(sys, &(P.getInternalType()));
527 }
528
529 //------------------------------------------------------------------------------
530 // PRE-SOLVE DIAGNOSTICS
531
532 /**
533 Get a list of variables to fix to make an underspecified system
534 become square. Also seems to return stuff when you have a structurally
535 singuler system.
536 */
537 vector<Variable>
538 Simulation::getFixableVariables(){
539 //cerr << "GETTING FIXABLE VARIABLES..." << endl;
540 vector<Variable> vars;
541
542 if(!sys){
543 throw runtime_error("Simulation system not yet built");
544 }
545
546 int32 *vip; /** TODO ensure 32 bit integers are used */
547
548 // Get IDs of elegible variables in array at vip...
549 CONSOLE_DEBUG("Calling slvDOF_eligible");
550 if(!slvDOF_eligible(sys,&vip)){
551 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
552 }else{
553 struct var_variable **vp = slv_get_solvers_var_list(sys);
554
555 if(vp==NULL){
556 throw runtime_error("Simulation variable list is null");
557 }
558
559 // iterate through this list until we find a -1:
560 int i=0;
561 int var_index = vip[i];
562 while(var_index >= 0){
563 struct var_variable *var = vp[var_index];
564 vars.push_back( Variable(this, var) );
565 ++i;
566 var_index = vip[i];
567 }
568 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
569 ascfree(vip);
570 }
571
572 return vars;
573 }
574
575 /**
576 Get the list of variables near their bounds. Helps to indentify why
577 you might be having non-convergence problems.
578 */
579 vector<Variable>
580 Simulation::getVariablesNearBounds(const double &epsilon){
581 //cerr << "GETTING VARIABLES NEAR BOUNDS..." << endl;
582 vector<Variable> vars;
583
584 if(!sys){
585 throw runtime_error("Simulation system not yet built");
586 }
587
588 int *vip;
589 CONSOLE_DEBUG("Calling slv_near_bounds...");
590 if(slv_near_bounds(sys,epsilon,&vip)){
591 struct var_variable **vp = slv_get_solvers_var_list(sys);
592 struct var_variable *var;
593 cerr << "VARS FOUND NEAR BOUNDS" << endl;
594 int nlow = vip[0];
595 int nhigh = vip[1];
596 int lim1 = 2 + nlow;
597 for(int i=2; i<lim1; ++i){
598 var = vp[vip[i]];
599 char *var_name = var_make_name(sys,var);
600 cerr << "AT LOWER BOUND: " << var_name << endl;
601 ascfree(var_name);
602 vars.push_back(Variable(this,var));
603 };
604 int lim2 = lim1 + nhigh;
605 for(int i=lim1; i<lim2; ++i){
606 var = vp[vip[i]];
607 char *var_name = var_make_name(sys,var);
608 cerr << "AT UPPER BOUND: " << var_name << endl;
609 ascfree(var_name);
610 vars.push_back(Variable(this,var));
611 }
612 }
613 ascfree(vip);
614 return vars;
615 }
616
617
618 bool
619 SingularityInfo::isSingular() const{
620 if(vars.size()||rels.size()){
621 return true;
622 }
623 return false;
624 }
625
626 //------------------------------------------------------------------------------
627 // SOLVING
628
629 /**
630 Solve the system through to convergence. This function is hardwired with
631 a maximum of 1000 iterations, but will interrupt itself when the 'stop'
632 condition comes back from the SolverReporter.
633 */
634 void
635 Simulation::solve(Solver solver, SolverReporter &reporter){
636
637 if(!sys){
638 try{
639 build();
640 }catch(runtime_error &e){
641 stringstream ss;
642 ss << "Unable to build system: ";
643 ss << e.what();
644 throw runtime_error(ss.str());
645 }
646 }
647
648 setSolver(solver);
649
650 //cerr << "PRESOLVING SYSTEM...";
651 CONSOLE_DEBUG("Calling slv_presolve...");
652 slv_presolve(sys);
653 //cerr << "DONE" << endl;
654
655 //cerr << "SOLVING SYSTEM..." << endl;
656 // Add some stuff here for cleverer iteration....
657 unsigned niter = 1000;
658 //double updateinterval = 0.02;
659
660 double starttime = tm_cpu_time();
661 //double lastupdate = starttime;
662 SolverStatus status;
663 //int solved_vars=0;
664 bool stop=false;
665
666 status.getSimulationStatus(*this);
667 reporter.report(&status);
668
669 for(unsigned iter = 1; iter <= niter && !stop; ++iter){
670
671 if(status.isReadyToSolve()){
672 /* CONSOLE_DEBUG("Calling slv_iterate..."); */
673 slv_iterate(sys);
674 }
675
676 status.getSimulationStatus(*this);
677
678 if(reporter.report(&status)){
679 stop = true;
680 }
681 }
682
683 double elapsed = tm_cpu_time() - starttime;
684 CONSOLE_DEBUG("Elapsed time: %0.3f", elapsed);
685
686 activeblock = status.getCurrentBlockNum();
687
688 reporter.finalise(&status);
689
690 // Just a little bit of console output:
691
692 if(status.isOK()){
693 //cerr << "... SOLVED, STATUS OK" << endl;
694 }else{
695 cerr << "... SOLVER FAILED" << endl;
696 }
697
698 //cerr << "SOLVER PERFORMED " << status.getIterationNum() << " ITERATIONS IN " << elapsed << "s" << endl;
699
700 // communicate solver variable status back to the instance tree
701 processVarStatus();
702 }
703
704 //------------------------------------------------------------------------------
705 // POST-SOLVE DIAGNOSTICS
706
707 const int
708 Simulation::getActiveBlock() const{
709 return activeblock;
710 }
711
712 /**
713 Return an IncidenceMatrix built from the current state of the solver system.
714
715 This will actually return something meaningful even before solve.
716 */
717 IncidenceMatrix
718 Simulation::getIncidenceMatrix(){
719 return IncidenceMatrix(*this);
720 }
721
722 /**
723 This function looks at all the variables in the solve's list and updates
724 the variable status for the corresponding instances.
725
726 It does this by using the 'interface pointer' in the Instance, see
727 the C-API function GetInterfacePtr.
728
729 This is used to display visually which variables have been solved, which
730 ones have not yet been attempted, and which ones were active when the solver
731 failed (ASCXX_VAR_ACTIVE).
732 */
733 void
734 Simulation::processVarStatus(){
735
736 // this is a cheap function call:
737 const mtx_block_t *bb = slv_get_solvers_blocks(getSystem());
738
739 var_variable **vlist = slv_get_solvers_var_list(getSystem());
740 int nvars = slv_get_num_solvers_vars(getSystem());
741
742 slv_status_t status;
743 slv_get_status(getSystem(), &status);
744
745 if(status.block.number_of == 0){
746 cerr << "Variable statuses can't be set: block structure not yet determined." << endl;
747 return;
748 }
749
750 int activeblock = status.block.current_block;
751 int low = bb->block[activeblock].col.low;
752 int high = bb->block[activeblock].col.high;
753 bool allsolved = status.converged;
754 for(int c=0; c < nvars; ++c){
755 var_variable *v = vlist[c];
756 Instanc i((Instance *)var_instance(v));
757 VarStatus s = ASCXX_VAR_STATUS_UNKNOWN;
758 if(i.isFixed()){
759 s = ASCXX_VAR_FIXED;
760 }else if(var_incident(v) && var_active(v)){
761 if(allsolved || c < low){
762 s = ASCXX_VAR_SOLVED;
763 }else if(c <= high){
764 s = ASCXX_VAR_ACTIVE;
765 }else{
766 s = ASCXX_VAR_UNSOLVED;
767 }
768 }
769 i.setVarStatus(s);
770 }
771 }
772

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