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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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