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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 34631 byte(s)
Setting up web subdirectory in repository
1 aw0a 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