/[ascend]/trunk/ascend4/solver/slv9a.c
ViewVC logotype

Contents of /trunk/ascend4/solver/slv9a.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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