/[ascend]/trunk/solvers/lrslv/slv9a.c
ViewVC logotype

Contents of /trunk/solvers/lrslv/slv9a.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2649 - (show annotations) (download) (as text)
Wed Dec 12 12:39:25 2012 UTC (11 years, 2 months ago) by jpye
File MIME type: text/x-csrc
File size: 34010 byte(s)
Fixing GPL header, removing postal address (rpmlint incorrect-fsf-address)
1 /* ASCEND modelling environment
2 Copyright (C) 1997-2007 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 Logical Relation Solver module, used by CMSlv.
18 *//*
19 Logical Relation Solver
20 by Vicente Rico-Ramirez
21 Created: 04/97
22 */
23
24 #include <math.h>
25
26 #include <ascend/general/ascMalloc.h>
27 #include <ascend/general/tm_time.h>
28 #include <ascend/general/mathmacros.h>
29 #include <ascend/general/mem.h>
30 #include <ascend/general/list.h>
31
32 #include <ascend/system/calc.h>
33 #include <ascend/system/relman.h>
34 #include <ascend/system/logrelman.h>
35 #include <ascend/system/bndman.h>
36 #include <ascend/system/slv_stdcalls.h>
37 #include <ascend/system/cond_config.h>
38
39 #include <ascend/solver/solver.h>
40 #include <ascend/solver/logblock.h>
41
42 # define HAVE_LRSLV 1
43
44 typedef struct slv9a_system_structure *slv9a_system_t;
45
46 ASC_DLLSPEC SolverRegisterFn lrslv_register;
47
48 /* Some constants for use with IDA */
49 #define IDA_TRUE = 2
50 #define IDA_FALSE = 3
51
52 #define SLV9A(s) ((slv9a_system_t)(s))
53 #define SERVER (sys->slv)
54 #define slv9a_PA_SIZE 7 /* MUST INCREMENT WHEN ADDING PARAMETERS */
55 #define SHOW_MORE_IMPT_PTR (sys->parm_array[0])
56 #define SHOW_MORE_IMPT ((*(int32 *)SHOW_MORE_IMPT_PTR))
57 #define SHOW_LESS_IMPT_PTR (sys->parm_array[1])
58 #define SHOW_LESS_IMPT ((*(int32 *)SHOW_LESS_IMPT_PTR))
59 #define AUTO_RESOLVE_PTR (sys->parm_array[2])
60 #define AUTO_RESOLVE ((*(int32 *)AUTO_RESOLVE_PTR))
61 #define TIME_LIMIT_PTR (sys->parm_array[3])
62 #define TIME_LIMIT ((*(int32 *)TIME_LIMIT_PTR))
63 #define ITER_LIMIT_PTR (sys->parm_array[4])
64 #define ITER_LIMIT ((*(int32 *)ITER_LIMIT_PTR))
65 #define PERTURB_BOUNDARY_PTR (sys->parm_array[5])
66 #define PERTURB_BOUNDARY ((*(int32 *)PERTURB_BOUNDARY_PTR))
67 #define WITH_IDA_PTR (sys->parm_array[6])
68 #define WITH_IDA ((*(int32 *)WITH_IDA_PTR))
69
70 /*
71 * auxiliar structures
72 */
73 struct structural_data {
74 mtx_matrix_t mtx; /* For use in structural analysis */
75 unsigned *subregions; /* Set of subregion indeces */
76 dof_t *dofdata; /* dof data pointer from server */
77 mtx_region_t reg; /* Current block region */
78 int32 rank; /* Rank of the matrix */
79 };
80
81
82 /*
83 * solver structure
84 */
85 struct slv9a_system_structure {
86
87 /*
88 * Problem definition
89 */
90 slv_system_t slv; /* slv_system_t back-link */
91 struct dis_discrete **vlist; /* Dis vars list (NULL terminated) */
92 struct logrel_relation **rlist; /* Logrelations list(NULL terminated) */
93 struct bnd_boundary **blist; /* Boundaries. Maybe NULL */
94
95 /*
96 * Solver information
97 */
98 int integrity; /* ? Has the system been created */
99 int32 presolved; /* ? Has the system been presolved */
100 slv_parameters_t p; /* parameters */
101 slv_status_t s; /* Status (as of iteration end) */
102 int32 cap; /* Order of matrix/vectors */
103 int32 rank; /* Symbolic rank of problem */
104 int32 vused; /* Free and incident variables */
105 int32 vtot; /* length of varlist */
106 int32 rused; /* Included relations */
107 int32 rtot; /* length of rellist */
108 double clock; /* CPU time */
109
110 void *parm_array[slv9a_PA_SIZE];
111 struct slv_parameter pa[slv9a_PA_SIZE];
112 /*
113 * Calculated data
114 */
115 struct structural_data S; /* structural information */
116 };
117
118
119 /*
120 * Integrity checks
121 * ----------------
122 * check_system(sys)
123 */
124
125 #define OK ((int32)813029392)
126 #define DESTROYED ((int32)103289182)
127
128 static int check_system(slv9a_system_t sys)
129 /*
130 * Checks sys for NULL and for integrity.
131 */
132 {
133 if( sys == NULL ) {
134 FPRINTF(ASCERR,"ERROR: (slv9a) check_system\n");
135 FPRINTF(ASCERR," NULL system handle.\n");
136 return 1;
137 }
138
139 switch( sys->integrity ) {
140 case OK:
141 return 0;
142 case DESTROYED:
143 FPRINTF(ASCERR,"ERROR: (slv9a) check_system\n");
144 FPRINTF(ASCERR," System was recently destroyed.\n");
145 return 1;
146 default:
147 FPRINTF(ASCERR,"ERROR: (slv9a) check_system\n");
148 FPRINTF(ASCERR," System reused or never allocated.\n");
149 return 1;
150 }
151 }
152
153 /*
154 * General input/output routines (discrete vars and log rels)
155 * ----------------------------------------------------------
156 * print_dis_name(out,sys,dvar)
157 * print_logrel_name(out,sys,lrel)
158 */
159
160 #define print_dis_name(a,b,c) slv_print_dis_name((a),(b)->slv,(c))
161 #define print_logrel_name(a,b,c) slv_print_logrel_name((a),(b)->slv,(c))
162
163 /*
164 * Debug output routines
165 * ---------------------
166 * debug_delimiter(fp)
167 * debug_out_dvar_values(fp,sys)
168 * debug_out_logrel_residuals(fp,sys)
169 * debug_out_structural_mtx(fp,sys)
170 */
171
172 static void debug_delimiter( FILE *fp)
173 /*
174 * Outputs a hyphenated line.
175 */
176 {
177 int i;
178 for( i=0; i<60; i++ ) PUTC('-',fp);
179 PUTC('\n',fp);
180 }
181
182 #if DEBUG
183
184 /*
185 * Outputs all variable values in current block.
186 */
187 static void debug_out_dvar_values( FILE *fp, slv9a_system_t sys)
188 {
189 int32 col;
190 struct dis_discrete *dvar;
191
192 FPRINTF(fp,"Discrete var values --> \n");
193 for( col = sys->S.reg.col.low; col <= sys->S.reg.col.high ; col++ ) {
194 dvar = sys->vlist[mtx_col_to_org(sys->S.mtx,col)];
195 print_dis_name(fp,sys,dvar);
196 FPRINTF(fp, "\nI Value Col \n");
197 FPRINTF(fp,"%d\t%d\t%d\n",dis_sindex(dvar),dis_value(dvar),col);
198 }
199 }
200
201
202 /*
203 * Outputs all logical relation residuals in current block.
204 */
205 static void debug_out_logrel_residuals( FILE *fp, slv9a_system_t sys)
206 {
207 int32 row;
208
209 FPRINTF(fp,"Logical rel residuals --> \n");
210 for( row = sys->S.reg.row.low; row <= sys->S.reg.row.high ; row++ ) {
211 struct logrel_relation *lrel;
212 lrel = sys->rlist[mtx_row_to_org(sys->S.mtx,row)];
213 FPRINTF(fp," %d : ",logrel_residual(lrel));
214 print_logrel_name(fp,sys,lrel);
215 PUTC('\n',fp);
216 }
217 PUTC('\n',fp);
218 }
219
220
221 /*
222 * Outputs permutation of the elements in the structural matrix.
223 */
224 static void debug_out_structural_mtx( FILE *fp, slv9a_system_t sys)
225 {
226 mtx_coord_t nz;
227
228 nz.row = sys->S.reg.row.low;
229 for( ; nz.row <= sys->S.reg.row.high; ++(nz.row) ) {
230 FPRINTF(fp," Row %d (lrel %d)\n", nz.row,
231 mtx_row_to_org(sys->S.mtx,nz.row));
232 }
233 }
234
235 #endif /* DEBUG */
236
237 /*
238 * Array operations
239 * -----------------
240 * destroy_array(p)
241 * create_array(len,type)
242 * create_zero_array(len,type)
243 */
244 #define destroy_array(p) \
245 if( (p) != NULL ) ascfree((p))
246 #define create_array(len,type) \
247 ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
248 #define create_zero_array(len,type) \
249 ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
250
251
252 /*
253 * Block routines
254 * --------------
255 * block_feasible(sys)
256 * move_to_next_block(sys)
257 * find_next_unconverged_block(sys)
258 */
259
260 /*
261 * Returns TRUE if the current block is feasible, FALSE otherwise.
262 * It is assumed that the residuals have been computed.
263 */
264 static boolean block_feasible( slv9a_system_t sys)
265 {
266 int32 row;
267 struct logrel_relation *rel;
268
269 if( !sys->s.calc_ok )
270 return(FALSE);
271
272 for( row = sys->S.reg.row.low; row <= sys->S.reg.row.high; row++ ) {
273 rel = sys->rlist[mtx_row_to_org(sys->S.mtx,row)];
274 if( !logrel_satisfied(rel) ) return FALSE;
275 }
276 return TRUE;
277 }
278
279
280 /*
281 * Moves on to the next block, updating all of the solver information.
282 * To move to the first block, set sys->s.block.current_block to -1 before
283 * calling. If already at the last block, then sys->s.block.current_block
284 * will equal the number of blocks and the system will be declared
285 * converged.
286 */
287 static void move_to_next_block( slv9a_system_t sys)
288 {
289 struct dis_discrete *dvar;
290 struct logrel_relation *lrel;
291 int32 row;
292 int32 col;
293 int32 ci;
294
295 if( sys->s.block.current_block >= 0 ) {
296
297 /* Record cost accounting info here. */
298 ci=sys->s.block.current_block;
299 sys->s.cost[ci].size = sys->s.block.current_size;
300 sys->s.cost[ci].iterations = sys->s.block.iteration;
301 sys->s.cost[ci].funcs = sys->s.block.funcs;
302 sys->s.cost[ci].jacs = sys->s.block.jacs;
303 sys->s.cost[ci].functime = sys->s.block.functime;
304 sys->s.cost[ci].jactime = sys->s.block.jactime;
305 sys->s.cost[ci].time = sys->s.block.cpu_elapsed;
306 sys->s.cost[ci].resid = sys->s.block.residual;
307
308 /* De-initialize previous block */
309 if (SHOW_LESS_IMPT && (sys->s.block.current_size >1)) {
310 FPRINTF(LIF(sys),"Block %d converged.\n",
311 sys->s.block.current_block);
312 }
313
314 for( col=sys->S.reg.col.low; col <= sys->S.reg.col.high; col++ ) {
315 dvar = sys->vlist[mtx_col_to_org(sys->S.mtx,col)];
316 dis_set_in_block(dvar,FALSE);
317 }
318
319 for( row=sys->S.reg.row.low; row <= sys->S.reg.row.high; row++ ) {
320 lrel = sys->rlist[mtx_row_to_org(sys->S.mtx,row)];
321 logrel_set_in_block(lrel,FALSE);
322 }
323 sys->s.block.previous_total_size += sys->s.block.current_size;
324
325 }
326
327 sys->s.block.current_block++;
328 if( sys->s.block.current_block < sys->s.block.number_of ) {
329
330 /* Initialize next block */
331 sys->S.reg =
332 (slv_get_solvers_log_blocks(SERVER))->block[sys->s.block.current_block];
333
334 row = sys->S.reg.row.high - sys->S.reg.row.low + 1;
335 col = sys->S.reg.col.high - sys->S.reg.col.low + 1;
336 sys->s.block.current_size = MAX(row,col);
337
338 sys->s.block.iteration = 0;
339 sys->s.block.cpu_elapsed = 0.0;
340 sys->s.block.functime = 0.0;
341 sys->s.block.jactime = 0.0;
342 sys->s.block.funcs = 0;
343 sys->s.block.jacs = 0;
344
345 if(SHOW_LESS_IMPT && ( sys->s.block.current_size > 1)) {
346 debug_delimiter(LIF(sys));
347 }
348 if(SHOW_LESS_IMPT) {
349 FPRINTF(LIF(sys),"\n%-40s ---> %d in [%d..%d]\n",
350 "Current block number", sys->s.block.current_block,
351 0, sys->s.block.number_of-1);
352 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Current block size",
353 sys->s.block.current_size);
354 }
355 sys->s.calc_ok = TRUE;
356
357 } else {
358 /*
359 * Before we claim convergence, we must check if we left behind
360 * some unassigned logrelations. If and only if they happen to be
361 * satisfied at the current point, convergence has been obtained.
362 */
363 if( sys->s.struct_singular ) {
364 sys->s.block.current_size = sys->rused - sys->rank;
365 if(SHOW_LESS_IMPT) {
366 debug_delimiter(LIF(sys));
367 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Unassigned Log Rels",
368 sys->s.block.current_size);
369 }
370 if( block_feasible(sys) ) {
371 if(SHOW_LESS_IMPT) {
372 FPRINTF(LIF(sys),"\nUnassigned logrels ok. You lucked out.\n");
373 }
374 sys->s.converged = TRUE;
375 } else {
376 if(SHOW_LESS_IMPT) {
377 FPRINTF(LIF(sys),"\nProblem inconsistent: %s.\n",
378 "Unassigned logrels not satisfied");
379 }
380 sys->s.inconsistent = TRUE;
381 }
382 if(SHOW_LESS_IMPT) {
383 debug_delimiter(LIF(sys));
384 }
385 } else {
386 sys->s.converged = TRUE;
387 }
388 }
389 }
390
391
392 /*
393 * Moves to next unconverged block, assuming that the current block has
394 * converged (or is -1, to start).
395 */
396 static void find_next_unconverged_block( slv9a_system_t sys)
397 {
398 do {
399 move_to_next_block(sys);
400 #if DEBUG
401 debug_out_dvar_values(ASCERR,sys);
402 debug_out_logrel_residuals(ASCERR,sys);
403 #endif /* DEBUG */
404 } while( !sys->s.converged && block_feasible(sys) );
405 }
406
407
408 /*
409 * Iteration begin/end routines
410 * ----------------------------
411 * iteration_begins(sys)
412 * iteration_ends(sys)
413 */
414
415 /*
416 * Prepares sys for entering an iteration, increasing the iteration counts
417 * and starting the clock.
418 */
419 static void iteration_begins( slv9a_system_t sys)
420 {
421 sys->clock = tm_cpu_time();
422 ++(sys->s.block.iteration);
423 ++(sys->s.iteration);
424 if(SHOW_LESS_IMPT && (sys->s.block.current_size >1 )) {
425 FPRINTF(LIF(sys),"\n%-40s ---> %d\n",
426 "Iteration", sys->s.block.iteration);
427 FPRINTF(LIF(sys),"%-40s ---> %d\n",
428 "Total iteration", sys->s.iteration);
429 }
430 }
431
432
433 /*
434 * Prepares sys for exiting an iteration, stopping the clock and recording
435 * the cpu time.
436 */
437 static void iteration_ends( slv9a_system_t sys)
438 {
439 double cpu_elapsed; /* elapsed this iteration */
440
441 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
442 sys->s.block.cpu_elapsed += cpu_elapsed;
443 sys->s.cpu_elapsed += cpu_elapsed;
444 if(SHOW_LESS_IMPT && (sys->s.block.current_size >1 )) {
445 FPRINTF(LIF(sys),"%-40s ---> %g\n",
446 "Elapsed time", sys->s.block.cpu_elapsed);
447 FPRINTF(LIF(sys),"%-40s ---> %g\n",
448 "Total elapsed time", sys->s.cpu_elapsed);
449 }
450 }
451
452
453 /*
454 * Updates the solver status.
455 */
456 static void update_status( slv9a_system_t sys)
457 {
458 boolean unsuccessful;
459
460 if( !sys->s.converged ) {
461 sys->s.time_limit_exceeded =
462 (sys->s.block.cpu_elapsed >= TIME_LIMIT);
463 sys->s.iteration_limit_exceeded =
464 (sys->s.block.iteration >= ITER_LIMIT);
465 }
466
467 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
468 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded
469 || (sys->s.block.current_size >1 );
470
471 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
472 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
473 }
474
475
476 static int32 slv9a_get_default_parameters(slv_system_t server,
477 SlvClientToken asys,
478 slv_parameters_t *parameters)
479 {
480 slv9a_system_t sys;
481 union parm_arg lo,hi,val;
482 struct slv_parameter *new_parms = NULL;
483 int32 make_macros = 0;
484
485 if (server != NULL && asys != NULL) {
486 sys = SLV9A(asys);
487 make_macros = 1;
488 }
489
490 if (parameters->parms == NULL) {
491 /* an external client wants our parameter list.
492 * an instance of slv9a_system_structure has this pointer
493 * already set in slv9a_create
494 */
495 new_parms = (struct slv_parameter *)
496 ascmalloc((slv9a_PA_SIZE)*sizeof(struct slv_parameter));
497 if (new_parms == NULL) {
498 return -1;
499 }
500 parameters->parms = new_parms;
501 parameters->dynamic_parms = 1;
502 }
503 parameters->num_parms = 0;
504
505 /* begin defining parameters */
506
507 slv_define_parm(parameters, bool_parm,
508 "showmoreimportant", "showmoreimportant", "showmoreimportant",
509 U_p_bool(val,1),U_p_bool(lo,0),U_p_bool(hi,1),-1);
510 SLV_BPARM_MACRO(SHOW_MORE_IMPT_PTR,parameters);
511
512 slv_define_parm(parameters, bool_parm,
513 "showlessimportant", "detailed solving info",
514 "detailed solving info",
515 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
516 SLV_BPARM_MACRO(SHOW_LESS_IMPT_PTR,parameters);
517
518 slv_define_parm(parameters, bool_parm,
519 "autoresolve", "auto-resolve", "auto-resolve",
520 U_p_bool(val,1),U_p_bool(lo,0),U_p_bool(hi,1), 2);
521 SLV_BPARM_MACRO(AUTO_RESOLVE_PTR,parameters);
522
523 slv_define_parm(parameters, int_parm,
524 "timelimit", "time limit (CPU sec/block)",
525 "time limit (CPU sec/block)",
526 U_p_int(val,1500),U_p_int(lo, 1),U_p_int(hi,20000),1);
527 SLV_IPARM_MACRO(TIME_LIMIT_PTR,parameters);
528
529 slv_define_parm(parameters, int_parm,
530 "iterationlimit", "max iterations/block",
531 "max iterations/block",
532 U_p_int(val, 30),U_p_int(lo, 1),U_p_int(hi,20000),1);
533 SLV_IPARM_MACRO(ITER_LIMIT_PTR,parameters);
534
535 slv_define_parm(parameters, bool_parm,
536 "perturbboundaries", "perturb boundaries",
537 "perturb boundaries",
538 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), -1);
539 SLV_BPARM_MACRO(PERTURB_BOUNDARY_PTR,parameters);
540
541 slv_define_parm(parameters, bool_parm,
542 "withida", "LRSlv called by IDA",
543 "LRSlv called by IDA",
544 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), -1);
545 SLV_BPARM_MACRO(WITH_IDA_PTR,parameters);
546 return 1;
547 }
548
549
550
551 /*
552 * External routines
553 * -----------------
554 * See slv_client.h
555 */
556
557 static SlvClientToken slv9a_create(slv_system_t server, int *statusindex)
558 {
559 slv9a_system_t sys;
560
561 sys = (slv9a_system_t)asccalloc(1, sizeof(struct slv9a_system_structure) );
562 if (sys==NULL) {
563 *statusindex = 1;
564 return sys;
565 }
566 SERVER = server;
567 sys->p.parms = sys->pa;
568 sys->p.dynamic_parms = 0;
569 slv9a_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
570 sys->integrity = OK;
571 sys->presolved = 0;
572 sys->p.output.more_important = stdout;
573 sys->p.output.less_important = stdout;
574 sys->p.whose = (*statusindex);
575
576 sys->s.ok = TRUE;
577 sys->s.calc_ok = TRUE;
578 sys->s.costsize = 0;
579 sys->s.cost = NULL; /*redundant, but sanity preserving */
580 sys->vlist = slv_get_solvers_dvar_list(server);
581 sys->rlist = slv_get_solvers_logrel_list(server);
582 sys->blist = slv_get_solvers_bnd_list(server);
583 if (sys->vlist == NULL) {
584 ascfree(sys);
585 FPRINTF(ASCERR,"LRSlv called with no discrete variables.\n");
586 *statusindex = -2;
587 return NULL;
588 }
589 if (sys->rlist == NULL) {
590 ascfree(sys);
591 FPRINTF(ASCERR,"LRSlv called with no logical relations.\n");
592 *statusindex = -1;
593 return NULL;
594 }
595 slv_check_dvar_initialization(server);
596 *statusindex = 0;
597 return((SlvClientToken)sys);
598 }
599
600
601 static void destroy_matrices( slv9a_system_t sys)
602 {
603 if (sys->S.mtx) {
604 mtx_destroy(sys->S.mtx);
605 sys->S.mtx = NULL;
606 }
607 return;
608 }
609
610
611 static int slv9a_eligible_solver(slv_system_t server)
612 {
613 logrel_filter_t lrfilter;
614
615 lrfilter.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
616 lrfilter.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
617 if (slv_count_solvers_logrels(server,&lrfilter)) {
618 return(TRUE);
619 } else {
620 return(FALSE);
621 }
622 }
623
624
625 static void slv9a_get_parameters(slv_system_t server, SlvClientToken asys,
626 slv_parameters_t *parameters)
627 {
628 slv9a_system_t sys;
629 sys = SLV9A(asys);
630 if (server == NULL || sys==NULL) return;
631 if (check_system(sys)) return;
632 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
633 }
634
635
636 static void slv9a_set_parameters(slv_system_t server, SlvClientToken asys,
637 slv_parameters_t *parameters)
638 {
639 slv9a_system_t sys;
640 sys = SLV9A(asys);
641 if (server == NULL || sys==NULL) return;
642 if (check_system(sys)) return;
643 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
644 }
645
646
647
648 static int slv9a_get_status(slv_system_t server, SlvClientToken asys
649 ,slv_status_t *status
650 ){
651 slv9a_system_t sys;
652 (void) server;
653 sys = SLV9A(asys);
654 if (check_system(sys)) return 1;
655 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
656 return 0;
657 }
658
659 /*
660 * Performs structural analysis on the system, setting the flags in
661 * status. The problem must be set up, the logrelation/dis_var list
662 * must be non-NULL. The structural matrix must be created and have
663 * the correct order (stored in sys->cap).
664 * Everything else will be determined here.
665 * On entry there isn't yet a correspondence between dis_sindex and
666 * structural matrix column. Here we establish that.
667 */
668 static void structural_analysis(slv_system_t server, slv9a_system_t sys)
669 {
670 dis_filter_t dvfilter;
671 logrel_filter_t lrfilter;
672
673 /*
674 * The server has marked incidence flags already.
675 */
676 /* count included equalities */
677 lrfilter.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
678 lrfilter.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
679 sys->rused = slv_count_solvers_logrels(server,&lrfilter);
680
681 /* count free and incident vars */
682 dvfilter.matchbits = (DIS_FIXED | DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
683 dvfilter.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
684 sys->vused = slv_count_solvers_dvars(server,&dvfilter);
685
686 /* Symbolic analysis */
687 sys->rtot = slv_get_num_solvers_logrels(server);
688 sys->vtot = slv_get_num_solvers_dvars(server);
689 if (sys->rtot) {
690 if (slv_log_block_partition(server)) {
691 FPRINTF(ASCERR,
692 "Structural Analysis:Error in slv_log_block_partition\n");
693 return;
694 }
695 }
696 sys->S.dofdata = slv_get_log_dofdata(server);
697 sys->rank = sys->S.dofdata->structural_rank;
698
699 /* Initialize Status */
700 sys->s.over_defined = (sys->rused > sys->vused);
701 sys->s.under_defined = (sys->rused < sys->vused);
702 sys->s.struct_singular = (sys->rank < sys->rused);
703 sys->s.block.number_of = (slv_get_solvers_log_blocks(SERVER))->nblocks;
704 }
705
706
707 /*
708 * Create matrix to perform symbolic analysis
709 */
710 static void create_matrices(slv_system_t server, slv9a_system_t sys)
711 {
712 sys->S.mtx = mtx_create();
713 mtx_set_order(sys->S.mtx,sys->cap);
714 structural_analysis(server,sys);
715 }
716
717
718 /*
719 * Here we will check if any fixed or included flags have
720 * changed since the last presolve.
721 */
722 static int32 slv9a_dof_changed(slv9a_system_t sys)
723 {
724 int32 ind, result = 0;
725 /*
726 * Currently we have two copies of the fixed and included flags
727 * which must be kept in sync. The dis_fixed and logrel_included
728 * functions perform the syncronization and hence must be called
729 * over the whole dvar list and logrel list respectively. When we move
730 * to using only one set of flags (bit flags) this function can
731 * be changed to return 1 at the first indication of a change
732 * in the dof.
733 */
734
735 /* search for dvars that were fixed and are now free */
736 for( ind = sys->vused; ind < sys->vtot; ++ind ) {
737 if( !dis_fixed(sys->vlist[ind]) && dis_active(sys->vlist[ind]) ) {
738 ++result;
739 }
740 }
741 /* search for logrels that were unincluded and are now included */
742 for( ind = sys->rused; ind < sys->rtot; ++ind ) {
743 if( logrel_included(sys->rlist[ind]) && logrel_active(sys->rlist[ind])) {
744 ++result;
745 }
746 }
747 /* search for dvars that were free and are now fixed */
748 for( ind = sys->vused -1; ind >= 0; --ind ) {
749 if( dis_fixed(sys->vlist[ind]) || !dis_active(sys->vlist[ind])) {
750 ++result;
751 }
752 }
753 /* search for logrels that were included and are now unincluded */
754 for( ind = sys->rused -1; ind >= 0; --ind ) {
755 if(!logrel_included(sys->rlist[ind]) || !logrel_active(sys->rlist[ind])){
756 ++result;
757 }
758 }
759 return result;
760 }
761
762
763 static void reset_cost(struct slv_block_cost *cost,int32 costsize){
764 int32 ind;
765 for( ind = 0; ind < costsize; ++ind ) {
766 cost[ind].size = 0;
767 cost[ind].iterations = 0;
768 cost[ind].funcs = 0;
769 cost[ind].jacs = 0;
770 cost[ind].functime = 0;
771 cost[ind].jactime = 0;
772 cost[ind].time = 0;
773 cost[ind].resid = 0;
774 }
775 }
776
777
778 static int slv9a_presolve(slv_system_t server, SlvClientToken asys){
779 struct dis_discrete **dvp;
780 struct logrel_relation **lrp;
781 int32 cap, ind;
782 int32 matrix_creation_needed = 1;
783 slv9a_system_t sys;
784
785 sys = SLV9A(asys);
786 iteration_begins(sys);
787 check_system(sys);
788 if( sys->vlist == NULL ) {
789 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Discrete Variable list was never set.");
790 return 1;
791 }
792 if( sys->rlist == NULL ) {
793 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Logical Relation list was never set.");
794 return 2;
795 }
796
797 if(sys->presolved > 0) { /* system has been presolved before */
798 if(!slv9a_dof_changed(sys) ) { /* no changes in fixed or included flags */
799 #if DEBUG
800 FPRINTF(ASCERR,"Avoiding matrix destruction/creation\n");
801 #endif /* DEBUG */
802 matrix_creation_needed = 0;
803 }
804 }
805
806 lrp=sys->rlist;
807 for( ind = 0; ind < sys->rtot; ++ind ) {
808 logrel_set_satisfied(lrp[ind],FALSE);
809 }
810 if( matrix_creation_needed ) {
811 cap = slv_get_num_solvers_logrels(SERVER);
812 sys->cap = slv_get_num_solvers_dvars(SERVER);
813 sys->cap = MAX(sys->cap,cap);
814 dvp=sys->vlist;
815 for( ind = 0; ind < sys->vtot; ++ind ) {
816 dis_set_in_block(dvp[ind],FALSE);
817 }
818 lrp=sys->rlist;
819 for( ind = 0; ind < sys->rtot; ++ind ) {
820 logrel_set_in_block(lrp[ind],FALSE);
821 logrel_set_satisfied(lrp[ind],FALSE);
822 }
823 sys->presolved = 1; /* full presolve recognized here */
824 destroy_matrices(sys);
825 create_matrices(server,sys);
826 sys->s.block.current_reordered_block = -2;
827 }
828
829 /* Reset status */
830 sys->s.iteration = 0;
831 sys->s.cpu_elapsed = 0.0;
832 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
833 sys->s.block.previous_total_size = 0;
834 sys->s.costsize = 1+sys->s.block.number_of;
835
836 if( matrix_creation_needed ) {
837 destroy_array(sys->s.cost);
838 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
839 for( ind = 0; ind < sys->s.costsize; ++ind ) {
840 sys->s.cost[ind].reorder_method = -1;
841 }
842 } else {
843 reset_cost(sys->s.cost,sys->s.costsize);
844 }
845
846 /* set to go to first unconverged block */
847 sys->s.block.current_block = -1;
848 sys->s.block.current_size = 0;
849 sys->s.calc_ok = TRUE;
850 sys->s.block.iteration = 0;
851
852 update_status(sys);
853 iteration_ends(sys);
854 sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;
855 return 0;
856 }
857
858 static int slv9a_resolve(slv_system_t server, SlvClientToken asys){
859 struct dis_discrete **dvp;
860 struct logrel_relation **lrp;
861 slv9a_system_t sys;
862
863 sys = SLV9A(asys);
864 (void) server;
865 check_system(sys);
866
867 for( dvp = sys->vlist ; *dvp != NULL ; ++dvp ) {
868 dis_set_in_block(*dvp,FALSE);
869 }
870 for( lrp = sys->rlist ; *lrp != NULL ; ++lrp ) {
871 logrel_set_in_block(*lrp,FALSE);
872 logrel_set_satisfied(*lrp,FALSE);
873 }
874
875 /* Reset status */
876 sys->s.iteration = 0;
877 sys->s.cpu_elapsed = 0.0;
878 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
879 sys->s.block.previous_total_size = 0;
880
881 /* go to first unconverged block */
882 sys->s.block.current_block = -1;
883 sys->s.block.current_size = 0;
884 sys->s.calc_ok = TRUE;
885 sys->s.block.iteration = 0;
886
887 update_status(sys);
888 return 0;
889 }
890
891
892 /*
893 * The boundary, relation and logrelation structures in this function
894 * are used to perform the solution of logical relations when changing
895 * the value of truth of the boundary expressions. This is sometimes
896 * required for conditional analysis. The use of the structure instance
897 * in this function is an insanity, but we stick with it by now.
898 */
899 static int slv9a_iterate(slv_system_t server, SlvClientToken asys){
900 slv9a_system_t sys;
901 struct bnd_boundary **blist;
902 struct bnd_boundary *cur_bnd;
903 struct rel_relation *rel;
904 struct logrel_relation *logrel;
905 struct gl_list_t *per_insts;
906 struct Instance *i;
907 bnd_filter_t bfilter;
908 int32 numbnds,numper,nb;
909
910 FILE *mif;
911 FILE *lif;
912 int ds_status=0;
913 int per_value;
914 double time0;
915
916 sys = SLV9A(asys);
917 mif = MIF(sys);
918 lif = LIF(sys);
919 if(server == NULL || sys==NULL)return 1;
920 if(check_system(SLV9A(sys)))return 2;
921 if(!sys->s.ready_to_solve){
922 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not ready to solve.");
923 return 3;
924 }
925
926 /*
927 To change truth values of some boundaries
928 */
929 blist = sys->blist;
930 if (blist == NULL && PERTURB_BOUNDARY) {
931 ERROR_REPORTER_HERE(ASC_USER_ERROR,"No boundaries in the problem."
932 " The solver cannot work in perturbation mode."
933 );
934 sys->s.ready_to_solve = FALSE;
935 iteration_ends(sys);
936 return 4;
937 }
938
939 if(blist == NULL && WITH_IDA) {
940 ERROR_REPORTER_HERE(ASC_USER_ERROR, "No boundaries found when"
941 "LRSlving with IDA");
942 sys->s.ready_to_solve = FALSE;
943 iteration_ends(sys);
944 return 4;
945 }
946
947 /*
948 Solution process begins
949 */
950 if(sys->s.block.current_block==-1) {
951 find_next_unconverged_block(sys);
952 update_status(sys);
953 return 0; /* is that right? */
954 }
955
956 /*
957 finding the list of boundaries to be perturbed
958 */
959 per_insts = NULL;
960 per_value = 0;
961 if(PERTURB_BOUNDARY) {
962 numbnds = slv_get_num_solvers_bnds(server);
963 bfilter.matchbits = (BND_PERTURB);
964 bfilter.matchvalue = (BND_PERTURB);
965 numper = slv_count_solvers_bnds(server,&bfilter);
966 if(numper != 0) {
967 per_insts = gl_create(numper);
968 per_value = 1;
969 for (nb=0; nb <numbnds; nb++){
970 cur_bnd = blist[nb];
971 if(bnd_perturb(cur_bnd)) {
972 if(bnd_kind(cur_bnd) == e_bnd_rel) {
973 rel = bnd_rel(bnd_real_cond(cur_bnd));
974 i = (struct Instance *)rel_instance(rel);
975 gl_append_ptr(per_insts,i);
976 }else{
977 if (bnd_kind(cur_bnd) == e_bnd_logrel) {
978 logrel = bnd_logrel(bnd_log_cond(cur_bnd));
979 i = (struct Instance *)logrel_instance(logrel);
980 gl_append_ptr(per_insts,i);
981 }
982 }
983 }
984 }
985 }
986 }
987
988 /* stick the boundary to be toggled onto per_insts. Should only be 1 instance
989 * in the list, i.e only single boundary crossings (for the moment....) */
990 if (WITH_IDA) {
991 numbnds = slv_get_num_solvers_bnds(server);
992 for (nb = 0; nb < numbnds; nb++) {
993 cur_bnd = blist[nb];
994 if (bnd_ida_crossed(cur_bnd)) {
995 per_insts = gl_create(1);
996 if(bnd_ida_value(cur_bnd)) {
997 per_value = 2;
998 } else {
999 per_value = 3;
1000 }
1001
1002 if (bnd_kind(cur_bnd) == e_bnd_rel) {
1003 rel = bnd_rel(bnd_real_cond(cur_bnd));
1004 i = (struct Instance *) rel_instance(rel);
1005 gl_append_ptr(per_insts, i);
1006 } else {
1007 if (bnd_kind(cur_bnd) == e_bnd_logrel) {
1008 logrel = bnd_logrel(bnd_log_cond(cur_bnd));
1009 i = (struct Instance *) logrel_instance(logrel);
1010 gl_append_ptr(per_insts, i);
1011 }
1012 }
1013 }
1014
1015 }
1016 }
1017
1018 iteration_begins(sys);
1019
1020 /*
1021 Attempt direct solve if appropriate
1022 */
1023 if( sys->s.block.current_size == 1 ) {
1024 struct dis_discrete *dvar;
1025 struct logrel_relation *lrel;
1026 dvar = sys->vlist[mtx_col_to_org(sys->S.mtx,sys->S.reg.col.low)];
1027 lrel = sys->rlist[mtx_row_to_org(sys->S.mtx,sys->S.reg.row.low)];
1028 if (SHOW_LESS_IMPT) {
1029 FPRINTF(lif,"%-40s ---> (%d)", "Singleton relation",
1030 mtx_row_to_org(sys->S.mtx,sys->S.reg.row.low));
1031 print_logrel_name(lif,sys,lrel); PUTC('\n',lif);
1032 FPRINTF(lif,"%-40s ---> (%d)", "Singleton variable",
1033 mtx_col_to_org(sys->S.mtx,sys->S.reg.col.low));
1034 print_dis_name(lif,sys,dvar); PUTC('\n',lif);
1035 }
1036
1037 /* Attempt direct solve */
1038 time0=tm_cpu_time();
1039 if (PERTURB_BOUNDARY && per_insts != NULL) {
1040 ds_status=slv_direct_log_solve(SERVER,lrel,dvar,mif,1,per_insts);
1041 gl_destroy(per_insts);
1042 per_insts = NULL;
1043 } else if (WITH_IDA && per_insts != NULL) {
1044 ds_status = slv_direct_log_solve(SERVER,lrel,dvar,mif,per_value,per_insts);
1045 gl_destroy(per_insts);
1046 per_insts = NULL;
1047 } else {
1048
1049 ds_status=slv_direct_log_solve(SERVER,lrel,dvar,mif,0,NULL);
1050 }
1051 sys->s.block.functime += (tm_cpu_time()-time0);
1052
1053 switch( ds_status ) {
1054 case 0:
1055 if (SHOW_LESS_IMPT) {
1056 FPRINTF(lif,"Unable to directly solve a logical relation.\n");
1057 }
1058 FPRINTF(lif,"Bad discrete variable or logrel\n");
1059 return 5;
1060 case 1:
1061 if (SHOW_LESS_IMPT) {
1062 FPRINTF(lif,"Directly solved.\n");
1063 }
1064 iteration_ends(sys);
1065 find_next_unconverged_block(sys);
1066 update_status(sys);
1067 return 0;
1068 case 2:
1069 if (SHOW_LESS_IMPT) {
1070 FPRINTF(lif,"Directly solved.\n");
1071 }
1072 sys->s.inconsistent = TRUE;
1073 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
1074 FPRINTF(ASCERR,"Multiple solution exists for the discrete variable '");
1075 print_dis_name(ASCERR,sys,dvar);
1076 FPRINTF(ASCERR,"' when solving the logical relation '");
1077 print_logrel_name(ASCERR,sys,lrel);
1078 FPRINTF(ASCERR,"'.");
1079 error_reporter_end_flush();
1080 iteration_ends(sys);
1081 update_status(sys);
1082 return 6;
1083 case -1:
1084 sys->s.inconsistent = TRUE;
1085 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
1086 FPRINTF(ASCERR,"No solution exists for the discrete variable '");
1087 print_dis_name(ASCERR,sys,dvar);
1088 FPRINTF(ASCERR,"' when solving the logical relation '");
1089 print_logrel_name(ASCERR,sys,lrel);
1090 FPRINTF(ASCERR,"'.");
1091 error_reporter_end_flush();
1092 iteration_ends(sys);
1093 update_status(sys);
1094 return 7;
1095 }
1096 } else {
1097 FPRINTF(lif,"block number = %d \n",sys->s.block.current_block);
1098 FPRINTF(lif,"block size = %d \n",sys->s.block.current_size );
1099 FPRINTF(lif,"block iteration = %d \n",sys->s.block.iteration);
1100 FPRINTF(lif,"Multiple logical relations in block.\n");
1101 FPRINTF(lif,"Not supported in this solver.\n");
1102 iteration_ends(sys);
1103 update_status(sys);
1104 }
1105
1106 #if DEBUG
1107 FPRINTF(ASCERR,"***********end of iteration*****************\n");
1108 debug_out_dvar_values(LIF(sys), sys);
1109 debug_out_logrel_residuals(LIF(sys), sys);
1110 FPRINTF(ASCERR,"********************************************\n");
1111 #endif /* DEBUG */
1112 return 0;
1113 }
1114
1115
1116 static int slv9a_solve(slv_system_t server, SlvClientToken asys){
1117 slv9a_system_t sys;
1118 int err = 0;
1119 sys = SLV9A(asys);
1120 if (server == NULL || sys==NULL) return -1;
1121 if (check_system(sys)) return -2;
1122 while( sys->s.ready_to_solve )err |= slv9a_iterate(server,sys);
1123 return err;
1124 }
1125
1126
1127 static mtx_matrix_t slv9a_get_structural_matrix(slv_system_t server,
1128 SlvClientToken sys)
1129 {
1130 if (server == NULL || sys==NULL) return NULL;
1131 if (check_system(SLV9A(sys))) return NULL;
1132 return SLV9A(sys)->S.mtx;
1133 }
1134
1135 static int slv9a_destroy(slv_system_t server, SlvClientToken asys)
1136 {
1137 slv9a_system_t sys;
1138 sys = SLV9A(asys);
1139 if (server == NULL || sys==NULL) return 1;
1140 if (check_system(sys)) return 1;
1141 slv_destroy_parms(&(sys->p));
1142 destroy_matrices(sys);
1143 sys->integrity = DESTROYED;
1144 if (sys->s.cost) ascfree(sys->s.cost);
1145 ascfree( (POINTER)asys );
1146 return 0;
1147 }
1148
1149 static const SlvFunctionsT slv9a_internals = {
1150 99
1151 ,"LRSlv"
1152 ,slv9a_create
1153 ,slv9a_destroy
1154 ,slv9a_eligible_solver
1155 ,slv9a_get_default_parameters
1156 ,slv9a_get_parameters
1157 ,slv9a_set_parameters
1158 ,slv9a_get_status
1159 ,slv9a_solve
1160 ,slv9a_presolve
1161 ,slv9a_iterate
1162 ,slv9a_resolve
1163 ,NULL
1164 ,slv9a_get_structural_matrix
1165 ,NULL
1166 };
1167
1168 int lrslv_register(void){
1169 return solver_register(&slv9a_internals);
1170 }
1171

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