/[ascend]/trunk/base/generic/solver/slv9.c
ViewVC logotype

Contents of /trunk/base/generic/solver/slv9.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1175 - (show annotations) (download) (as text)
Thu Jan 18 00:34:22 2007 UTC (13 years, 10 months ago) by johnpye
File MIME type: text/x-csrc
File size: 165955 byte(s)
A few easily-fixed warnings dispelled
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Conditional Modeling Solver (CMSlv) module.
21 *//*
22 Conditional Modeling Solver
23 by Vicente Rico-Ramirez
24 Last in CVS: $Revision: 1.22 $ $Date: 2000/01/25 02:27:58 $ $Author: ballan $
25 */
26
27 #include <math.h>
28
29 #include <utilities/config.h>
30 #include <utilities/ascConfig.h>
31 #include <utilities/ascSignal.h>
32 #include <utilities/ascMalloc.h>
33 #include <general/tm_time.h>
34 #include <utilities/mem.h>
35 #include <compiler/compiler.h>
36 #include <compiler/instance_enum.h>
37 #include <general/list.h>
38 #include <general/mathmacros.h>
39 #include <compiler/fractions.h>
40 #include <compiler/dimen.h>
41 #include <compiler/functype.h>
42 #include <compiler/func.h>
43 #include "mtx.h"
44 #include "mtx_reorder.h"
45 #include "linsol.h"
46 #include "linsolqr.h"
47 #include "slv_types.h"
48 #include "var.h"
49 #include "rel.h"
50 #include "discrete.h"
51 #include "conditional.h"
52 #include "logrel.h"
53 #include "bnd.h"
54 #include "calc.h"
55 #include "relman.h"
56 #include "logrelman.h"
57 #include "bndman.h"
58 #include "cond_config.h"
59 #include "slv_common.h"
60 #include "slv_client.h"
61 #include "slv9.h"
62 #include "slv_stdcalls.h"
63 #include "slvDOF.h"
64
65 #ifdef ASC_WITH_CONOPT
66 # include <solver/conopt.h>
67 #else
68 # define MAX_INT MAXINT
69 # define MAX_REAL MAXDOUBLE
70 #endif
71
72 #if !defined(STATIC_CMSLV) && !defined(DYNAMIC_CMSLV)
73 int slv9_register(SlvFunctionsT *f)
74 {
75 UNUSED_PARAMETER(f);
76
77 ERROR_REPORTER_HERE(ASC_PROG_ERR,"CMSlv was not built in this copy of ASCEND");
78 return 1;
79 }
80 #else /* either STATIC_CMSLV or DYNAMIC_CMSLV is defined */
81
82 #ifdef DYNAMIC_CMSLV
83 /* do dynamic loading stuff. yeah, right */
84
85 #else /* following is used if STATIC_CMSLV is defined */
86
87
88 /*
89 * definitions to enable/disable the output of partial results in
90 * the solution of a problem
91 */
92 #define DEBUG FALSE
93 #define SHOW_LOGICAL_DETAILS FALSE
94 #define SHOW_BOUNDARY_ANALYSIS_DETAILS FALSE
95 #define SHOW_OPTIMIZATION_DETAILS FALSE
96 #define SHOW_BISECTION_DETAILS FALSE
97 #define SHOW_LINEAR_SEARCH_DETAILS FALSE
98 #define SHOW_LAGRANGE_DETAILS FALSE
99 #define DEBUG_CONSISTENCY FALSE
100 #define TEST_CONSISTENCY FALSE
101 #define USE_CONSISTENCY FALSE
102
103 /*
104 * system definitions
105 */
106 #define SLV9(s) ((slv9_system_t)(s))
107 #define SERVER (sys->slv)
108 #define slv9_PA_SIZE 26 /* MUST INCREMENT WHEN ADDING PARAMETERS */
109 #define LOGSOLVER_OPTION_PTR (sys->parm_array[0])
110 #define LOGSOLVER_OPTION ((*(char **)LOGSOLVER_OPTION_PTR))
111 #define NONLISOLVER_OPTION_PTR (sys->parm_array[1])
112 #define NONLISOLVER_OPTION ((*(char **)NONLISOLVER_OPTION_PTR))
113 #define OPTSOLVER_OPTION_PTR (sys->parm_array[2])
114 #define OPTSOLVER_OPTION ((*(char **)OPTSOLVER_OPTION_PTR))
115 #define TIME_LIMIT_PTR (sys->parm_array[3])
116 #define TIME_LIMIT ((*(int32 *)TIME_LIMIT_PTR))
117 #define ITER_LIMIT_PTR (sys->parm_array[4])
118 #define ITER_LIMIT ((*(int32 *)ITER_LIMIT_PTR))
119 #define ITER_BIS_LIMIT_PTR (sys->parm_array[5])
120 #define ITER_BIS_LIMIT ((*(int32 *)ITER_BIS_LIMIT_PTR))
121 #define TOO_SMALL_PTR (sys->parm_array[6])
122 #define TOO_SMALL ((*(real64 *)TOO_SMALL_PTR))
123 #define LINEAR_SEARCH_FACTOR_PTR (sys->parm_array[7])
124 #define LINEAR_SEARCH_FACTOR ((*(real64 *)LINEAR_SEARCH_FACTOR_PTR))
125 #define SHOW_MORE_IMPT_PTR (sys->parm_array[8])
126 #define SHOW_MORE_IMPT ((*(int32 *)SHOW_MORE_IMPT_PTR))
127 #define SHOW_LESS_IMPT_PTR (sys->parm_array[9])
128 #define SHOW_LESS_IMPT ((*(int32 *)SHOW_LESS_IMPT_PTR))
129 #define AUTO_RESOLVE_PTR (sys->parm_array[10])
130 #define AUTO_RESOLVE ((*(int32 *)AUTO_RESOLVE_PTR))
131 #define UNDEFINED_PTR (sys->parm_array[11])
132 #define UNDEFINED ((*(real64 *)UNDEFINED_PTR))
133 #define DOMLIM_PTR (sys->parm_array[12])
134 #define DOMLIM ((*(int32 *)DOMLIM_PTR))
135 #define OPT_ITER_LIMIT_PTR (sys->parm_array[13])
136 #define OPT_ITER_LIMIT ((*(int32 *)OPT_ITER_LIMIT_PTR))
137 #define INFINITY_PTR (sys->parm_array[14])
138 #define ASC_INFINITY ((*(real64 *)INFINITY_PTR))
139 #define OBJ_TOL_PTR (sys->parm_array[15])
140 #define OBJ_TOL ((*(real64 *)OBJ_TOL_PTR))
141 #define RTMAXJ_PTR (sys->parm_array[16])
142 #define RTMAXJ ((*(real64 *)RTMAXJ_PTR))
143 #define RHO_PTR (sys->parm_array[17])
144 #define RHO ((*(real64 *)RHO_PTR))
145
146
147 /*
148 * Client tokens of the different solvers: Conditional, Optimizer,
149 * Nonlinear, Logical. We will switch from one client token to
150 * another as the solution process occurs.
151 */
152 #define NUMBER_OF_CLIENTS 4
153 SlvClientToken token[NUMBER_OF_CLIENTS];
154 int32 solver_index[NUMBER_OF_CLIENTS];
155
156 /*
157 * indeces in arrays token and solver_index
158 */
159 #define CONDITIONAL_SOLVER 0
160 #define LOGICAL_SOLVER 1
161 #define NONLINEAR_SOLVER 2
162 #define OPTIMIZATION_SOLVER 3
163
164 /*
165 * Do we have an optimization problem ?. Global variable initialized
166 * to 0 (not optimizing)
167 */
168 static int32 g_optimizing = 0;
169
170 #if USE_CONSISTENCY
171 /*
172 * number of subregion visited during the solution of the conditional
173 * model.
174 */
175 static int32 g_subregions_visited;
176
177 #endif /* USE_CONSISTENCY */
178
179
180 /* auxiliar structures */
181 struct boolean_values {
182 int32 *pre_val; /* previous values of dis_discrete */
183 int32 *cur_val; /* current values of dis_discrete */
184 };
185
186 struct matching_cases {
187 int32 *case_list; /* list of cases */
188 int32 ncases; /* number of cases */
189 int32 diff_subregion; /* subregion ? */
190 };
191
192 struct real_values {
193 real64 *pre_values; /* previous values of var_variables */
194 real64 *cur_values; /* current values of var_variables */
195 };
196
197 struct opt_vector {
198 real64 *element; /* elements in colum of matrix */
199 };
200
201 struct opt_matrix {
202 struct opt_vector *cols; /* columns in matrix */
203 };
204
205 struct subregionID {
206 unsigned long ID_number;
207 int32 *bool_values;
208 };
209
210 struct ds_subregion_list {
211 int32 length,capacity;
212 struct subregionID *sub_stack;
213 };
214
215 struct ds_subregions_visited {
216 int32 length,capacity;
217 unsigned long *visited;
218 };
219
220 /*
221 * This solver's data structure (CMSlv)
222 */
223 struct slv9_system_structure {
224
225 /*
226 * Problem definition
227 */
228 slv_system_t slv; /* slv_system_t back-link */
229
230 struct rel_relation *obj; /* Objective function: NULL = none */
231 struct var_variable **vlist; /* Variable list (NULL terminated) */
232 struct rel_relation **rlist; /* Relation list (NULL terminated) */
233 struct dis_discrete **dvlist; /* Dis vars list (NULL terminated) */
234 struct logrel_relation **lrlist; /* Logrels list(NULL terminated)*/
235 struct bnd_boundary **blist; /* Variable list (NULL terminated) */
236 struct var_variable **mvlist;
237 struct dis_discrete **mdvlist; /* We will not touch the masters list,
238 * but they can provide very useful
239 * information to the conditional
240 * solver since the master index does
241 * not change.
242 */
243
244 /*
245 * for optimization at boundaries
246 */
247 struct opt_matrix *coeff_matrix; /* Matrix for optimization problem */
248 struct opt_vector *opt_var_values; /* Values of vars in opt problem */
249 int32 subregions; /* number of subregions at cur bnd */
250 mtx_matrix_t lin_mtx; /* Matrix to define the linear system
251 * for calculation of the lagrange
252 * multipliers
253 */
254 /*
255 * For search consistency analysis
256 */
257 struct ds_subregion_list subregion_list;
258 /*
259 * Information about the subregions
260 * visited during the solution of the
261 * conditional model
262 */
263 struct ds_subregions_visited subregions_visited;
264 /*
265 * ID number of the subregions
266 * visited
267 */
268 int32 *bool_mindex; /* master indices of boolean vars in
269 * the problem associated with WHENs
270 */
271 int32 need_consistency_analysis; /* Is the consistency analysis needed */
272
273
274 /*
275 * Solver information
276 */
277 int32 integrity; /* Has the system been created ? */
278 int32 presolved; /* Has the system been presolved ? */
279 slv_parameters_t p; /* Parameters */
280 slv_status_t s; /* Status (as of iteration end) */
281 int32 cap; /* Order of matrix/vectors */
282 int32 rank; /* Symbolic rank of problem */
283 int32 vused; /* Free and incident variables */
284 int32 vtot; /* length of varlist */
285 int32 mvtot; /* length of master varlist */
286 int32 rused; /* Included relations */
287 int32 rtot; /* length of rellist */
288 real64 clock; /* CPU time */
289 int32 nliter; /* iterations in nonlinear solver */
290
291 void *parm_array[slv9_PA_SIZE]; /* array of pointers to param values */
292 struct slv_parameter pa[slv9_PA_SIZE]; /* &pa[0] => sys->p.parms */
293
294 #ifdef ASC_WITH_CONOPT
295 /*
296 * Data for optimizer at boundaries (CONOPT)
297 */
298 struct conopt_data con;
299 #endif
300 };
301
302
303 /*
304 * Integrity checks
305 * ----------------
306 * check_system(sys)
307 */
308
309 #define OK ((int32)813029392)
310 #define DESTROYED ((int32)103289182)
311
312 /*
313 * Checks sys for NULL and for integrity.
314 */
315 static int check_system(slv9_system_t sys)
316 {
317 if( sys == NULL ) {
318 FPRINTF(ASCERR,"ERROR: (slv9) check_system\n");
319 FPRINTF(ASCERR," NULL system handle.\n");
320 return 1;
321 }
322
323 switch( sys->integrity ) {
324 case OK:
325 return 0;
326 case DESTROYED:
327 FPRINTF(ASCERR,"ERROR: (slv9) check_system\n");
328 FPRINTF(ASCERR," System was recently destroyed.\n");
329 return 1;
330 default:
331 FPRINTF(ASCERR,"ERROR: (slv9) check_system\n");
332 FPRINTF(ASCERR," System reused or never allocated.\n");
333 return 1;
334 }
335 }
336
337
338 /*
339 * General input/output routines
340 * -----------------------------
341 * print_var_name(out,sys,var)
342 */
343 #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
344
345
346 /*
347 * Array operations
348 * ----------------
349 * destroy_array(p)
350 * create_array(len,type)
351 * create_zero_array(len,type)
352 */
353 #define destroy_array(p) \
354 if( (p) != NULL ) ascfree((p))
355 #define create_array(len,type) \
356 ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
357 #define create_zero_array(len,type) \
358 ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
359
360
361 /*
362 * Search Consistency Analysis during iterative process
363 * ---------------------------------------------------------
364 * The caller of this functions is in charge of
365 * defining the extension of the analysis by passing an integer
366 * which will tell us if 1) the analysis consider only the current
367 * and the previous alternatives or 2) the analysis consider all the
368 * alternatives visited at current the state of the solution
369 * algorithm.
370 */
371
372 /*
373 * Handling dynamic allocation of the structural information
374 */
375
376 #define alloc_array(nelts,type) \
377 ((nelts) > 0 ? (type *)ascmalloc((nelts)*sizeof(type)) : NULL)
378 #define copy_nums(from,too,nnums) \
379 asc_memcpy((from),(too),(nnums)*sizeof(int32))
380 #define copy_subregions(from,too,nsubs) \
381 asc_memcpy((from),(too),(nsubs)*sizeof(struct subregionID))
382
383 #if TEST_CONSISTENCY
384 /*
385 * Appends the subregion_visited into the list
386 */
387 static void append_subregion(struct ds_subregion_list *sl,
388 struct subregionID sub)
389 {
390 if( sl->length == sl->capacity ) {
391 int32 newcap;
392 struct subregionID *newlist;
393
394 newcap = sl->capacity + 10;
395 newlist = alloc_array(newcap,struct subregionID);
396 copy_subregions((char *)sl->sub_stack,(char *)newlist,sl->length);
397 if( sl->sub_stack != NULL ) {
398 ascfree(sl->sub_stack);
399 }
400 sl->sub_stack = newlist;
401 sl->capacity = newcap;
402 }
403
404 sl->sub_stack[sl->length++] = sub;
405 }
406
407 /*
408 * Appends the subregion_visited into the list
409 */
410 static void append_sub_visited(struct ds_subregions_visited *sv,
411 unsigned long sub_visited)
412 {
413 if( sv->length == sv->capacity ) {
414 int32 newcap;
415 unsigned long *newlist;
416
417 newcap = sv->capacity + 10;
418 newlist = alloc_array(newcap,unsigned long);
419 copy_nums((char *)sv->visited,(char *)newlist,sv->length);
420 if( sv->visited != NULL ) {
421 ascfree(sv->visited);
422 }
423 sv->visited = newlist;
424 sv->capacity = newcap;
425 }
426
427 sv->visited[sv->length++] = sub_visited;
428 }
429
430 static unsigned long powoftwo (int32 expo)
431 {
432 unsigned long val;
433 int32 c;
434
435 val = 1;
436 for (c=1; c<= expo; c++) {
437 val = val * 2;
438 }
439
440 return val;
441 }
442
443
444 /*
445 * Storage information (boolean values) about a subregion so that
446 * we can visit it later for interactive strucutral analysis
447 */
448 static void ID_and_storage_subregion_information(slv_system_t server,
449 SlvClientToken asys)
450 {
451 slv9_system_t sys;
452 struct dis_discrete **bvlist;
453 struct dis_discrete *cur_dis;
454 struct ds_subregion_list *sl;
455 struct ds_subregions_visited *sv;
456 struct subregionID *sub;
457 dis_filter_t dfilter;
458 unsigned long val, visited, sID;
459 int32 d, numdvs, numdvf, dcount;
460 int32 len, s, found;
461
462 sys = SLV9(asys);
463 check_system(sys);
464
465 bvlist = sys->mdvlist;
466 if(bvlist == NULL ) {
467 FPRINTF(ASCERR,"ERROR: ID_and_storage_subregion_information.\n");
468 FPRINTF(ASCERR," Master discrete var list was never set.\n");
469 return;
470 }
471 numdvs = slv_get_num_master_dvars(server);
472
473 dfilter.matchbits = (DIS_INWHEN | DIS_BOOLEAN);
474 dfilter.matchvalue = (DIS_INWHEN | DIS_BOOLEAN);
475 numdvf = slv_count_master_dvars(server,&dfilter);
476
477 if (numdvf > 0) {
478 sub = (struct subregionID *)(ascmalloc(sizeof(struct subregionID)));
479 sub->bool_values = (int32 *)(ascmalloc(numdvf*sizeof(int32)));
480 } else {
481 FPRINTF(ASCERR,"ERROR: ID_and_storage_subregion_information.\n");
482 FPRINTF(ASCERR," No boolean variables in the problem \n");
483 return;
484 }
485
486 dcount = 0;
487 val = 0;
488 for (d=0; d<numdvs; d++) {
489 cur_dis = bvlist[d];
490 if (dis_apply_filter(cur_dis,&dfilter)) {
491 sub->bool_values[dcount] = dis_value(cur_dis);
492 dcount++;
493 if (sub->bool_values[dcount - 1] == 1) {
494 val = val + powoftwo(numdvf - dcount);
495 }
496 }
497 }
498 if ( (val == 0 ) && (numdvf > 0) ) {
499 val = powoftwo(numdvf);
500 }
501 sub->ID_number = val;
502 #if DEBUG_CONSISTENCY
503 FPRINTF(ASCERR,"ID of alternative is %ul \n", val);
504 #endif /* DEBUG_CONSISTENCY */
505 visited = val;
506 found = 0;
507 len = sys->subregions_visited.length;
508 if (len > 0) {
509 for (s=0; s<len; s++) {
510 sID = sys->subregions_visited.visited[s];
511 if (visited == sID) {
512 found = 1;
513 break;
514 }
515 }
516 }
517
518 sv = &(sys->subregions_visited);
519 append_sub_visited(sv,visited);
520
521 if (found == 0) {
522 #if DEBUG_CONSISTENCY
523 FPRINTF(ASCERR,"Saving alternative\n");
524 #endif /* DEBUG_CONSISTENCY */
525 sl = &(sys->subregion_list);
526 append_subregion(sl,(*sub));
527 } else {
528 destroy_array(sub->bool_values);
529 ascfree(sub);
530 }
531
532 return;
533 }
534 #endif
535
536 /*
537 * Destroys subregion information
538 */
539 static void destroy_subregion_information(SlvClientToken asys)
540 {
541 slv9_system_t sys;
542 struct subregionID *sub;
543 int32 lens, s;
544
545 sys = SLV9(asys);
546 check_system(sys);
547
548 if (sys->subregions_visited.visited != NULL) {
549 destroy_array(sys->subregions_visited.visited);
550 }
551
552 lens = sys->subregion_list.length;
553 if (lens != 0) {
554 for (s=0; s<lens; s++) {
555 sub = &(sys->subregion_list.sub_stack[s]);
556 if (sub->bool_values != NULL) {
557 destroy_array(sub->bool_values);
558 }
559 }
560 }
561
562 if (sys->subregion_list.sub_stack != NULL) {
563 destroy_array(sys->subregion_list.sub_stack);
564 }
565
566 if (sys->bool_mindex != NULL) {
567 destroy_array(sys->bool_mindex);
568 }
569 }
570
571
572 #if 0 /** unused function eligible_set_for_neighboring_subregions */
573 /* might be used if DEBUG_CONSISTENCY on */
574
575 /*
576 * Storing original values of boolean variables
577 */
578 static void store_original_bool_values(struct gl_list_t *bollist,
579 struct boolean_values *bval)
580 {
581 struct dis_discrete *dvar;
582 int32 d, dlen;
583
584 dlen = gl_length(bollist);
585 bval->pre_val = create_array(dlen,int32);
586 bval->cur_val = create_array(dlen,int32);
587 for (d=1; d<=dlen; d++){
588 dvar = (struct dis_discrete *)gl_fetch(bollist,d);
589 bval->cur_val[d-1] = dis_value(dvar);
590 bval->pre_val[d-1] = dis_previous_value(dvar);
591 }
592 }
593
594 /*
595 * Restoring original values of boolean variables
596 */
597 static void restore_original_bool_values(struct gl_list_t *bollist,
598 struct boolean_values *bval)
599 {
600 struct dis_discrete *dvar;
601 int32 d, dlen;
602
603 dlen = gl_length(bollist);
604 for (d=1; d<=dlen; d++){
605 dvar = (struct dis_discrete *)gl_fetch(bollist,d);
606 dis_set_boolean_value(dvar,bval->cur_val[d-1]);
607 dis_set_value(dvar,bval->cur_val[d-1]);
608 dis_set_previous_value(dvar,bval->pre_val[d-1]);
609 }
610 destroy_array(bval->cur_val);
611 destroy_array(bval->pre_val);
612 }
613
614 #endif /* if 0 */
615
616 /*
617 * the first element of cur_cases is in position one. The result is
618 * the same array, but ordered and starting in position zero
619 */
620 static void cases_reorder(int32 *cur_cases, int32 *correct_cases, int32 ncases)
621 {
622 int32 cur_case,pos=0,tmp_num,c,ind;
623
624 for (c=1; c<=ncases; c++) {
625 tmp_num = 0;
626 for (ind=1; ind<=ncases; ind++) {
627 cur_case = cur_cases[ind];
628 if (tmp_num < cur_case) {
629 pos = ind;
630 tmp_num = cur_case;
631 }
632 }
633 cur_cases[pos] = 0;
634 correct_cases[ncases-c] = tmp_num;
635 }
636
637 return;
638 }
639
640 #if 0 /** unused function eligible_set_for_neighboring_subregions */
641 /* might appear if debug_consistency is true. */
642 /*
643 * Restoring orignal configuration of the system
644 */
645 static void restore_configuration(slv_system_t server,
646 struct gl_list_t *bollist)
647
648 {
649 int32 *cur_cases, *correct_cases;
650 int32 ncases;
651
652 cur_cases = cases_matching(bollist,&ncases);
653 correct_cases = create_array(ncases,int32);
654 cases_reorder(cur_cases,correct_cases,ncases);
655 set_active_rels_in_subregion(server,correct_cases,ncases,bollist);
656 set_active_vars_in_subregion(server);
657 destroy_array(cur_cases);
658 destroy_array(correct_cases);
659 }
660
661
662 /*
663 * Get the list of boolean variables in the problem which are
664 * associated with a WHEN
665 */
666 static struct gl_list_t *get_list_of_booleans(slv_system_t server,
667 SlvClientToken asys)
668 {
669 slv9_system_t sys;
670 struct dis_discrete **bvlist;
671 struct dis_discrete *cur_dis;
672 struct gl_list_t *boolvars;
673 dis_filter_t dfilter;
674 int32 numdvf, numdvs, d, dcount;
675
676 sys = SLV9(asys);
677 check_system(sys);
678
679 bvlist = sys->mdvlist;
680 if(bvlist == NULL ) {
681 FPRINTF(ASCERR,"ERROR: (slv9) get_list_of_booleans.\n");
682 FPRINTF(ASCERR," Master discrete var list was never set.\n");
683 return NULL;
684 }
685 numdvs = slv_get_num_master_dvars(server);
686
687 dfilter.matchbits = (DIS_INWHEN | DIS_BOOLEAN);
688 dfilter.matchvalue = (DIS_INWHEN | DIS_BOOLEAN);
689 numdvf = slv_count_master_dvars(server,&dfilter);
690
691 if (numdvf == 0) {
692 FPRINTF(ASCERR,"ERROR: (slv9) get_list_of_booleans.\n");
693 FPRINTF(ASCERR," No boolean variables in the problem \n");
694 return NULL;
695 }
696
697 sys->bool_mindex = (int32 *)(ascmalloc(numdvf*sizeof(int32)));
698 boolvars = gl_create(numdvf);
699
700 dcount = 0;
701 for (d=0; d<numdvs; d++) {
702 cur_dis = bvlist[d];
703 if (dis_apply_filter(cur_dis,&(dfilter))) {
704 gl_append_ptr(boolvars,cur_dis);
705 sys->bool_mindex[dcount] = d;
706 dcount++;
707 }
708 }
709
710 return boolvars;
711 }
712
713 #endif /* 0*/
714
715 /*
716 * Get the eligible var list for each alternative
717 * Return:
718 * 1 means everything went right
719 * 0 means the analysis has failed with the current parititioning
720 * -1 means a memory problem has occurred
721 */
722 static int32 get_eligible_set(slv_system_t server,struct gl_list_t *disvars,
723 int32 *terminate)
724 {
725 struct var_variable **vslist;
726 struct var_variable **vmlist;
727 struct var_variable *mvar, *svar;
728 var_filter_t vfilter;
729 int32 *cur_cases;
730 int32 *correct_cases;
731 int32 *vars;
732 int32 v, count=0, ind;
733 int32 ncases;
734 int32 mnum;
735 int32 status,dof;
736
737 vslist = slv_get_solvers_var_list(server);
738 vmlist = slv_get_master_var_list(server);
739 mnum = slv_get_num_master_vars(server);
740 for (v=0; v<mnum; v++) {
741 mvar = vmlist[v];
742 var_set_eligible_in_subregion(mvar,FALSE);
743 }
744
745 cur_cases = cases_matching(disvars,&ncases);
746 correct_cases = create_array(ncases,int32);
747 cases_reorder(cur_cases,correct_cases,ncases);
748 set_active_rels_in_subregion(server,correct_cases,ncases,disvars);
749 set_active_vars_in_subregion(server);
750 destroy_array(cur_cases);
751 destroy_array(correct_cases);
752
753 #if DEBUG_CONSISTENCY
754 FPRINTF(ASCERR,"Analyzing alternative:\n");
755 #endif /* DEBUG_CONSISTENCY */
756
757 if (!slvDOF_status(server,(&status),(&dof))) {
758 FPRINTF(ASCERR,"ERROR in combinatorial search\n");
759 FPRINTF(ASCERR,"Combinatorial search aborted\n");
760 return -1;
761 } else {
762 if (status == 3) {
763 #if DEBUG_CONSISTENCY
764 FPRINTF(ASCERR,"Alternative is structurally singular\n");
765 #endif /* DEBUG_CONSISTENCY */
766 (*terminate) = 0;
767 return 0;
768 } else {
769 if (status == 4) {
770 #if DEBUG_CONSISTENCY
771 FPRINTF(ASCERR,"Alternative is overspecified\n");
772 #endif /* DEBUG_CONSISTENCY */
773 (*terminate) = 0;
774 return 0;
775 }
776 }
777 }
778
779 if (status == 1) {
780 (*terminate) = 0;
781 #if DEBUG_CONSISTENCY
782 FPRINTF(ASCERR,"Alternative has % d degrees of freedom.\n", dof);
783
784 #endif /* DEBUG_CONSISTENCY */
785 if (slvDOF_eligible(server,&(vars))) {
786 count = 0;
787 while (vars[count] != -1) {
788 ind = vars[count];
789 svar = vslist[ind];
790 v = var_mindex(svar);
791 mvar = vmlist[v];
792 var_set_eligible_in_subregion(mvar,TRUE);
793 count++;
794 }
795 destroy_array(vars);
796 }
797 if (dof > count) {
798 #if DEBUG_CONSISTENCY
799 FPRINTF(ASCERR,
800 "Alternative does not have enough number of eligible vars\n");
801 #endif /* DEBUG_CONSISTENCY */
802 return 0;
803 }
804 }
805
806 if (status == 2) {
807 #if DEBUG_CONSISTENCY
808 FPRINTF(ASCERR,"Alternative is square.\n");
809 #endif /* DEBUG_CONSISTENCY */
810 }
811
812 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR
813 | VAR_ELIGIBLE_IN_SUBREGION);
814 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
815
816 for (v=0; v<mnum; v++) {
817 mvar = vmlist[v];
818 if(var_apply_filter(mvar,&vfilter)) {
819 var_set_eligible(mvar,FALSE);
820 }
821 var_set_eligible_in_subregion(mvar,FALSE);
822 }
823
824 return 1;
825 }
826
827 /*
828 * Get the eligible set of variables for each of the alternatives generated
829 * by modifying the values of the boolean variables with the values stored
830 * during the solution process
831 * Return:
832 * 1 means everything went right
833 * 0 means the analysis has failed with the current partitioning
834 * -1 means a memory problem or wierdness has occurred
835 */
836 static int32 do_search_alternatives(slv_system_t server, SlvClientToken asys,
837 struct gl_list_t *disvars,
838 int32 *terminate, int32 all_sub)
839 {
840 slv9_system_t sys;
841 struct dis_discrete *cur_dis;
842 struct subregionID *sub;
843 int32 *values = NULL;
844 int32 dlen, test;
845 int32 lens, lenv, v, s, d;
846 int32 result;
847 unsigned long visited, vID;
848
849
850 sys = SLV9(asys);
851
852 dlen = gl_length(disvars);
853 lenv = sys->subregions_visited.length;
854 lens = sys->subregion_list.length;
855
856 if (all_sub == 0) { /* current and previous subregion */
857 for (v=lenv-2; v<lenv; v++) {
858 test = 0;
859 vID = sys->subregions_visited.visited[v];
860 for (s=lens-1; s>=0; s--) {
861 sub = &(sys->subregion_list.sub_stack[s]);
862 visited = sub->ID_number;
863 if (vID == visited) {
864 values = sub->bool_values;
865 test = 1;
866 FPRINTF(ASCERR,"s = %d \n",s);
867 break;
868 }
869 }
870
871 if (test == 0) {
872 FPRINTF(ASCERR,"ERROR: (slv9) do_search_alternatives \n");
873 FPRINTF(ASCERR," subregion not found \n");
874 return -1;
875 }
876 #if DEBUG_CONSISTENCY
877 FPRINTF(ASCERR,"Alternative = %ul \n", vID);
878 #endif /* DEBUG_CONSISTENCY */
879 for (d=0; d<dlen; d++) {
880 assert(values != NULL); /* if null, test was 0 above and we returned, in theory */
881 cur_dis = (struct dis_discrete *)(gl_fetch(disvars,d+1));
882 if (values[d] == 1) {
883 dis_set_boolean_value(cur_dis,TRUE);
884 } else {
885 dis_set_boolean_value(cur_dis,FALSE);
886 }
887 }
888 result = get_eligible_set(server,disvars,terminate);
889 if (result != 1) {
890 return result;
891 }
892 }
893
894 } else { /* all visited subregions */
895
896 for (s=lens-1; s>=0; s--) {
897 sub = &(sys->subregion_list.sub_stack[s]);
898 values = sub->bool_values;
899 vID = sub->ID_number;
900 #if DEBUG_CONSISTENCY
901 FPRINTF(ASCERR,"Alternative = %ul \n", vID);
902 #endif /* DEBUG_CONSISTENCY */
903 for (d=0; d<dlen; d++) {
904 cur_dis = (struct dis_discrete *)(gl_fetch(disvars,d+1));
905 if (values[d] == 1) {
906 dis_set_boolean_value(cur_dis,TRUE);
907 } else {
908 dis_set_boolean_value(cur_dis,FALSE);
909 }
910 }
911 result = get_eligible_set(server,disvars,terminate);
912 if (result != 1) {
913 return result;
914 }
915 }
916 }
917
918 return 1;
919 }
920
921
922 /*
923 * Perform consistency analysis for the visited/current-previous subregions.
924 * If all_subs is 1, the analysis takes in account all of the subregions
925 * visited by the solution algorithm at the current point if the solution
926 * procedure. If all_subs is 0, the analysis is only for the current
927 * and previous subregion.
928 */
929 static int32 consistency(slv_system_t server, SlvClientToken asys,
930 struct gl_list_t *bollist,
931 int32 all_subs, int32 *terminate)
932 {
933 slv9_system_t sys;
934 struct var_variable **vmlist;
935 struct var_variable *mvar;
936 var_filter_t vfilter;
937 int32 *globeli = NULL;
938 int32 dlen;
939 int32 mnum, v, elnum;
940 int32 result;
941 int32 iter;
942
943 sys = SLV9(asys);
944 check_system(sys);
945
946 /*
947 * Initializing eligible bit for variables
948 */
949 vmlist = slv_get_master_var_list(server);
950 mnum = slv_get_num_master_vars(server);
951 for (v=0; v<mnum; v++) {
952 mvar = vmlist[v];
953 var_set_eligible(mvar,TRUE);
954 }
955
956 dlen = gl_length(bollist);
957
958 #if DEBUG_CONSISTENCY
959 FPRINTF(ASCERR,"S e a r c h i n g \n");
960 #endif /* DEBUG_CONSISTENCY */
961 result = do_search_alternatives(server,asys,bollist,terminate,all_subs);
962
963 if (result != 1) {
964 #if DEBUG_CONSISTENCY
965 FPRINTF(ASCERR,"returning failed search after S e a r c h \n");
966 #endif /* DEBUG_CONSISTENCY */
967 return result;
968 }
969
970 /*
971 * Getting the "globally" eligible variables
972 */
973 vfilter.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_ELIGIBLE | VAR_FIXED);
974 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ELIGIBLE);
975 elnum = slv_count_master_vars(server,&vfilter);
976
977 if (elnum > 0) {
978 globeli = ASC_NEW_ARRAY(int32,elnum);
979 elnum = 0;
980 for (v=0; v<mnum; v++) {
981 mvar = vmlist[v];
982 if(var_apply_filter(mvar,&vfilter)) {
983 #if DEBUG_CONSISTENCY
984 FPRINTF(ASCERR,"Eligible index = %d \n",v);
985 #endif /* DEBUG_CONSISTENCY */
986 globeli[elnum] = v;
987 elnum++;
988 }
989 }
990 }
991
992 /*
993 * Recursively analysis
994 */
995
996 if ((*terminate) == 1) {
997 if (elnum != 0) {
998 #if DEBUG_CONSISTENCY
999 FPRINTF(ASCERR,"All alternatives are square but the \n");
1000 FPRINTF(ASCERR,"Eligible set is not null\n");
1001 #endif /* DEBUG_CONSISTENCY */
1002 destroy_array(globeli);
1003 }
1004 return 1;
1005 } else {
1006 if (elnum == 0) {
1007 #if DEBUG_CONSISTENCY
1008 FPRINTF(ASCERR,"No globally eligible variables to be fixed.\n");
1009 #endif /* DEBUG_CONSISTENCY */
1010 return 0;
1011 }
1012
1013 for (v=0; v<elnum; v++) {
1014 iter = 1;
1015 mvar = vmlist[globeli[v]];
1016 var_set_fixed(mvar,TRUE);
1017 var_set_potentially_fixed(mvar,TRUE);
1018 #if DEBUG_CONSISTENCY
1019 FPRINTF(ASCERR,"Fixing index = %d \n",globeli[v]);
1020 FPRINTF(ASCERR,"N e s t e d S e a r c h \n");
1021 #endif /* DEBUG_CONSISTENCY */
1022 result = consistency(server,asys,bollist,all_subs,&iter);
1023
1024 if (result != 1) {
1025 #if DEBUG_CONSISTENCY
1026 FPRINTF(ASCERR,"%d eliminated\n",globeli[v]);
1027 #endif /* DEBUG_CONSISTENCY */
1028 var_set_fixed(mvar,FALSE);
1029 var_set_potentially_fixed(mvar,FALSE);
1030 continue;
1031 } else {
1032 if (iter == 1) {
1033 (*terminate) = 1;
1034 #if DEBUG_CONSISTENCY
1035 FPRINTF(ASCERR,"%d Acepted \n",globeli[v]);
1036 #endif /* DEBUG_CONSISTENCY */
1037 destroy_array(globeli);
1038 return 1;
1039 } else {
1040 var_set_fixed(mvar,FALSE);
1041 var_set_potentially_fixed(mvar,FALSE);
1042 continue;
1043 }
1044 }
1045 }
1046 destroy_array(globeli);
1047 #if DEBUG_CONSISTENCY
1048 FPRINTF(ASCERR,"returning 0 after nested search\n");
1049 #endif /* DEBUG_CONSISTENCY */
1050 return 0;
1051 }
1052 }
1053
1054 #if 0 /** unused function eligible_set_for_neighboring_subregions */
1055 /* might appear if debug_consistency is true. */
1056
1057 /*
1058 * Get a set of globally eligible variables. Eligible for all the subregions
1059 * visited, or for the previous and current subregions.
1060 */
1061 static int32 get_globally_eligible(slv_system_t server, SlvClientToken asys,
1062 struct gl_list_t *bollist,
1063 int32 all_subs, int32 **eliset)
1064 {
1065 slv9_system_t sys;
1066 struct var_variable **vmlist;
1067 struct var_variable *mvar;
1068 var_filter_t vfilter;
1069 int32 dlen;
1070 int32 mnum, v, elnum;
1071 int32 terminate;
1072 int32 result;
1073
1074 sys = SLV9(asys);
1075 check_system(sys);
1076 /*
1077 * Initializing eligible bit for variables
1078 */
1079 vmlist = slv_get_master_var_list(server);
1080 mnum = slv_get_num_master_vars(server);
1081 for (v=0; v<mnum; v++) {
1082 mvar = vmlist[v];
1083 var_set_eligible(mvar,TRUE);
1084 }
1085
1086 dlen = gl_length(bollist);
1087
1088 /*
1089 * initializing
1090 */
1091 *eliset = NULL;
1092 terminate = 1;
1093
1094 #if DEBUG_CONSISTENCY
1095 FPRINTF(ASCERR,"S e a r c h i n g \n");
1096 #endif /* DEBUG_CONSISTENCY */
1097 result = do_search_alternatives(server,asys,bollist,&terminate,all_subs);
1098
1099 if (result != 1) {
1100 if (terminate == 0) {
1101 #if DEBUG_CONSISTENCY
1102 FPRINTF(ASCERR,"ERROR: some alternatives are either singular or\n");
1103 FPRINTF(ASCERR,"overspecified. All the alternatives have to be\n");
1104 FPRINTF(ASCERR,
1105 "either square or underspecified to complete the analysis\n");
1106 #endif /* DEBUG_CONSISTENCY */
1107 }
1108 return 0;
1109 }
1110
1111 /*
1112 * Getting the "globally" eligible variables
1113 */
1114 vfilter.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_ELIGIBLE | VAR_FIXED);
1115 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ELIGIBLE);
1116 elnum = slv_count_master_vars(server,&vfilter);
1117
1118 *eliset = (int32 *)ascmalloc((elnum+1)*sizeof(int32));
1119 elnum = 0;
1120 for (v=0; v<mnum; v++) {
1121 mvar = vmlist[v];
1122 if(var_apply_filter(mvar,&vfilter)) {
1123 #if DEBUG_CONSISTENCY
1124 FPRINTF(ASCERR,"Eligible index = %d \n",v);
1125 FPRINTF(ASCERR,"Variable : \n");
1126 print_var_name(ASCERR,sys,mvar);
1127 #endif /* DEBUG_CONSISTENCY */
1128 (*eliset)[elnum] = v;
1129 elnum++;
1130 }
1131 }
1132 (*eliset)[elnum] = -1;
1133
1134 if (elnum == 0) {
1135 if (terminate == 0) {
1136 #if DEBUG_CONSISTENCY
1137 FPRINTF(ASCERR,
1138 "Some alternatives are underspecified, but there does\n");
1139 FPRINTF(ASCERR,"not exist a set of eligible variables consistent \n");
1140 FPRINTF(ASCERR,"with all the alternatives\n");
1141 #endif /* DEBUG_CONSISTENCY */
1142 } else {
1143 #if DEBUG_CONSISTENCY
1144 FPRINTF(ASCERR,"All alternatives are already square\n");
1145 #endif /* DEBUG_CONSISTENCY */
1146 }
1147 return 0;
1148 } else {
1149 if (terminate == 1) {
1150 #if DEBUG_CONSISTENCY
1151 FPRINTF(ASCERR,"All alternatives are square but the \n");
1152 FPRINTF(ASCERR,"Eligible set is not null\n");
1153 #endif /* DEBUG_CONSISTENCY */
1154 }
1155 }
1156 return 1;
1157 }
1158
1159
1160
1161 /*
1162 * Store and restore values of the boolean variables of the problem
1163 * and calls for the the set of globally eligible variables.If all_subs
1164 * is 1, the analysis takes in account all of the subregions visited
1165 * by the solution algorithm at the current point if the solution
1166 * procedure. If all_subs is 0, the analysis is only for the current
1167 * and previous subregion.
1168 */
1169 static int32 consistent_eligible_set_for_subregions(slv_system_t server,
1170 SlvClientToken asys,
1171 int32 **vlist,
1172 int32 all_subs)
1173 {
1174 struct gl_list_t *blist;
1175 struct boolean_values bval;
1176 int32 result;
1177
1178 if (server==NULL || vlist == NULL) {
1179 FPRINTF(ASCERR,
1180 "consistent_eligible_set_for_subregions called with NULL.\n");
1181 return 0;
1182 }
1183
1184 blist = get_list_of_booleans(server,asys);
1185
1186 if ( (blist == NULL) || (gl_length(blist) == 0) ) {
1187 FPRINTF(ASCERR,"ERROR: (slv9) consistent_eligible_set_for_subregions \n");
1188 FPRINTF(ASCERR," List of boolean vars could not be found\n");
1189 return 0;
1190 }
1191
1192 store_original_bool_values(blist,&(bval));
1193 result = get_globally_eligible(server,asys,blist,all_subs,vlist);
1194
1195 restore_original_bool_values(blist,&(bval));
1196 restore_configuration(server,blist);
1197 gl_destroy(blist);
1198
1199 if (result == 1) {
1200 return 1;
1201 } else {
1202 return 0;
1203 }
1204
1205 }
1206
1207 /*
1208 * Store and restore values of the boolean variables of the problem
1209 * and calls for the consistency analysis of the subregions.If all_subs
1210 * is 1, the analysis takes in account all of the subregions visited
1211 * by the solution algorithm at the current point if the solution
1212 * procedure. If all_subs is 0, the analysis is only for the current
1213 * and previous subregion.
1214 */
1215 static int32 analyze_subregions(slv_system_t server,SlvClientToken asys,
1216 int32 **vlist, int32 all_subs)
1217 {
1218 slv9_system_t sys;
1219 struct var_variable ** vmlist;
1220 struct var_variable *mvar;
1221 struct gl_list_t *blist;
1222 struct boolean_values bval;
1223 var_filter_t vfilter;
1224 int32 mnum, elnum, v;
1225 int32 result;
1226 int32 terminate;
1227
1228 sys = SLV9(asys);
1229 check_system(sys);
1230
1231 if (server==NULL || vlist == NULL) {
1232 FPRINTF(ASCERR,"(slv9) analyze_subregions called with NULL.\n");
1233 return 0;
1234 }
1235
1236 blist = get_list_of_booleans(server,asys);
1237 if ( (blist == NULL) || (gl_length(blist) == 0) ) {
1238 FPRINTF(ASCERR,"ERROR: (slv9) analyze_subregions \n");
1239 FPRINTF(ASCERR," List of boolean vars could not be found\n");
1240 return 0;
1241 }
1242
1243 store_original_bool_values(blist,&(bval));
1244 /*
1245 * initializing
1246 */
1247 terminate = 1;
1248 (*vlist) = NULL;
1249
1250 vmlist = slv_get_master_var_list(server);
1251 mnum = slv_get_num_master_vars(server);
1252
1253 vfilter.matchbits = (VAR_POTENTIALLY_FIXED);
1254 vfilter.matchvalue = (VAR_POTENTIALLY_FIXED);
1255
1256 result = consistency(server,asys,blist,all_subs,&terminate);
1257
1258 if (result == 1) {
1259 /*
1260 * Getting the set of eligible variables
1261 */
1262 elnum = slv_count_master_vars(server,&vfilter);
1263 *vlist = (int32 *)ascmalloc((elnum+1)*sizeof(int32));
1264 elnum = 0;
1265 for (v=0; v<mnum; v++) {
1266 mvar = vmlist[v];
1267 if(var_apply_filter(mvar,&vfilter)) {
1268 var_set_fixed(mvar,FALSE);
1269 var_set_potentially_fixed(mvar,FALSE);
1270 #if DEBUG_CONSISTENCY
1271 FPRINTF(ASCERR,"Variable in consistent set: \n");
1272 print_var_name(ASCERR,sys,mvar);
1273 #endif /* DEBUG_CONSISTENCY */
1274 (*vlist)[elnum] = v;
1275 elnum++;
1276 }
1277 }
1278 (*vlist)[elnum] = -1;
1279
1280 restore_original_bool_values(blist,&(bval));
1281 restore_configuration(server,blist);
1282 gl_destroy(blist);
1283 return 1;
1284 } else {
1285 for (v=0; v<mnum; v++) {
1286 mvar = vmlist[v];
1287 if(var_apply_filter(mvar,&vfilter)) {
1288 var_set_fixed(mvar,FALSE);
1289 var_set_potentially_fixed(mvar,FALSE);
1290 }
1291 }
1292 restore_original_bool_values(blist,&(bval));
1293 restore_configuration(server,blist);
1294 gl_destroy(blist);
1295 return 0;
1296 }
1297 }
1298
1299
1300 /*
1301 * Finds the globally eligible set of variables only for the current and
1302 * previous subregions
1303 */
1304 static int32 eligible_set_for_neighboring_subregions(slv_system_t server,
1305 SlvClientToken asys,
1306 int32 **vlist)
1307 {
1308 slv9_system_t sys;
1309
1310 sys = SLV9(asys);
1311 check_system(sys);
1312
1313 if( sys->mdvlist == NULL ) {
1314 FPRINTF(ASCERR,"ERROR: (slv9) eligible_set_for_neighboring_subregions\n");
1315 FPRINTF(ASCERR," Discrete Variable list was never set.\n");
1316 return 0;
1317 }
1318
1319 if (!(sys->need_consistency_analysis)) {
1320 FPRINTF(ASCERR,"Globally eligible set not necessary\n");
1321 FPRINTF(ASCERR,"All the subregions have the same structure \n");
1322 return 0;
1323 }
1324
1325 if (consistent_eligible_set_for_subregions(server,asys,vlist,0)) {
1326 return 1;
1327 }
1328
1329 return 0;
1330 }
1331
1332
1333 /*
1334 * Perform the consistency analysis algorithm only for the current and
1335 * previous subregions
1336 */
1337 static int32 consistency_for_neighboring_subregions(slv_system_t server,
1338 SlvClientToken asys,
1339 int32 **vlist)
1340 {
1341 slv9_system_t sys;
1342
1343 sys = SLV9(asys);
1344 check_system(sys);
1345
1346 if( sys->mdvlist == NULL ) {
1347 FPRINTF(ASCERR,"ERROR: (slv9) consistency_for_neighboring_subregions\n");
1348 FPRINTF(ASCERR," Discrete Variable list was never set.\n");
1349 return 0;
1350 }
1351
1352 if (!(sys->need_consistency_analysis)) {
1353 FPRINTF(ASCERR,"consistency_analysis is not required\n");
1354 FPRINTF(ASCERR,"All the subregions have the same structure \n");
1355 return 0;
1356 }
1357
1358 if (analyze_subregions(server,asys,vlist,0)) {
1359 return 1;
1360 }
1361
1362 return 0;
1363 }
1364
1365
1366
1367 /*
1368 * Consistency analysis for visisted subregions. This function
1369 * gets the subregions that the solution algorithm has visited and
1370 * analyzes them.
1371 */
1372 static int32 eligible_set_for_subregions(slv_system_t server,
1373 SlvClientToken asys,
1374 int32 **vlist)
1375 {
1376 slv9_system_t sys;
1377
1378 sys = SLV9(asys);
1379 check_system(sys);
1380
1381 if( sys->mdvlist == NULL ) {
1382 FPRINTF(ASCERR,"ERROR: (slv9) eligible_set_for_subregions \n");
1383 FPRINTF(ASCERR," Discrete Variable list was never set.\n");
1384 return 0;
1385 }
1386
1387 if (!(sys->need_consistency_analysis)) {
1388 FPRINTF(ASCERR,"Globally eligible set not necessary \n");
1389 FPRINTF(ASCERR,"All the subregions have the same structure \n");
1390 return 0;
1391 }
1392
1393 if (consistent_eligible_set_for_subregions(server,asys,vlist,1)) {
1394 return 1;
1395 }
1396
1397 return 0;
1398 }
1399
1400
1401 /*
1402 * Consistency analysis for visisted subregions. This function
1403 * gets the subregions that the solution algorithm has visited and
1404 * analyzes them.
1405 */
1406 static int32 consistency_analysis_for_subregions(slv_system_t server,
1407 SlvClientToken asys,
1408 int32 **vlist)
1409 {
1410 slv9_system_t sys;
1411
1412 sys = SLV9(asys);
1413 check_system(sys);
1414
1415 if( sys->mdvlist == NULL ) {
1416 FPRINTF(ASCERR,"ERROR: (slv9) consistency_analysis_for_subregions\n");
1417 FPRINTF(ASCERR," Discrete Variable list was never set.\n");
1418 return 0;
1419 }
1420
1421 if (!(sys->need_consistency_analysis)) {
1422 FPRINTF(ASCERR,"consistency_analysis is not required\n");
1423 FPRINTF(ASCERR,"All the subregions have the same structure \n");
1424 return 0;
1425 }
1426
1427 if (analyze_subregions(server,asys,vlist,1)) {
1428 return 1;
1429 }
1430
1431 return 0;
1432 }
1433
1434 #endif /*#if 0 unused functions */
1435
1436
1437 /*
1438 * Handling of solution of the Logical Equations
1439 * ---------------------------------------------------------
1440 * This is made this way because it is a process which will be
1441 * required very often.
1442 */
1443
1444 /*
1445 * Solution of the logical relations encountered in the system based on
1446 * the current values of the discrete variables.
1447 */
1448 static void solve_logical_relations(slv_system_t server)
1449 {
1450 slv_set_client_token(server,token[LOGICAL_SOLVER]);
1451 slv_set_solver_index(server,solver_index[LOGICAL_SOLVER]);
1452 slv_presolve(server);
1453 #if SHOW_LOGICAL_DETAILS
1454 FPRINTF(ASCERR,"Solving Logical Relations\n");
1455 #endif /* SHOW_LOGICAL_DETAILS */
1456 slv_solve(server);
1457 }
1458
1459
1460
1461 /*
1462 * Handling the modification of parameters in external solvers
1463 * ---------------------------------------------------------
1464 */
1465
1466 /*
1467 * different types of parameter values
1468 */
1469 union param_value {
1470 int32 i;
1471 real64 r;
1472 int32 b;
1473 char *c;
1474 };
1475
1476
1477 /*
1478 * Setting the value of a parameter in a subsidiary solver
1479 */
1480 static void set_param_in_solver(slv_system_t server, int32 solver,
1481 enum parm_type types, char *param,
1482 union param_value *value)
1483 {
1484 slv_parameters_t p;
1485 int32 len,length;
1486
1487 slv_set_client_token(server,token[solver]);
1488 slv_set_solver_index(server,solver_index[solver]);
1489 slv_get_parameters(server,&p);
1490 length = p.num_parms;
1491 for (len = 0; len < length; len++) {
1492 if (p.parms[len].type == types) {
1493 switch (p.parms[len].type) {
1494 case bool_parm:
1495 if (strcmp(param,p.parms[len].name) == 0) {
1496 p.parms[len].info.b.value = value->b;
1497 }
1498 break;
1499 case real_parm:
1500 if (strcmp(param,p.parms[len].name) == 0) {
1501 p.parms[len].info.r.value = value->r;
1502 }
1503 break;
1504 case char_parm:
1505 if (strcmp(param,p.parms[len].name) == 0) {
1506 p.parms[len].info.c.value = value->c;
1507 }
1508 break;
1509 case int_parm:
1510 if (strcmp(param,p.parms[len].name) == 0) {
1511 p.parms[len].info.i.value = value->i;
1512 }
1513 break;
1514 default:
1515 break;
1516 }
1517 }
1518 }
1519 return;
1520 }
1521
1522
1523 /*
1524 * Analysis of Discrete Variables
1525 * -------------------------------
1526 */
1527
1528 /*
1529 * Compare current values of the discrete variables with their previous values
1530 * in order to know if some of them have changed.
1531 */
1532 static int32 some_dis_vars_changed(slv_system_t server, SlvClientToken asys)
1533 {
1534 struct dis_discrete **dv, *cur_dis;
1535 int32 numdvs, ind;
1536 slv9_system_t sys;
1537
1538 sys = SLV9(asys);
1539 check_system(sys);
1540
1541 if( sys->dvlist == NULL ) {
1542 FPRINTF(ASCERR,"ERROR: (slv9) some_dis_vars_changed\n");
1543 FPRINTF(ASCERR," Discrete variable list was never set.\n");
1544 return 0;
1545 }
1546
1547 dv = sys->dvlist;
1548 numdvs = slv_get_num_solvers_dvars(server);
1549 for( ind = 0; ind < numdvs; ind++ ) {
1550 cur_dis = dv[ind];
1551 #if SHOW_LOGICAL_DETAILS
1552 FPRINTF(ASCERR,"Boundary index = %d \n",ind);
1553 FPRINTF(ASCERR,"Current Value = %d\n",dis_value(cur_dis));
1554 FPRINTF(ASCERR,"Previous Value = %d\n",dis_previous_value(cur_dis));
1555 #endif /* SHOW_LOGICAL_DETAILS */
1556 if ((dis_kind(cur_dis)==e_dis_boolean_t ) && dis_inwhen(cur_dis) ) {
1557 if (dis_value(cur_dis) != dis_previous_value(cur_dis)) {
1558 return 1;
1559 }
1560 }
1561 }
1562 return 0;
1563 }
1564
1565 /*
1566 * Compare the original value of a discrete boolean variable (before
1567 * perturbation of boundaries) with its values after a solution
1568 * of the logical relations with some perturbed values for boundaries.
1569 * If those values are different, the bit VAL_MODIFIED is set to
1570 * TRUE. This will give us the boolean variable which will change as a
1571 * consequence of a boundary crossing.
1572 */
1573 static void search_for_modified_dvars(struct dis_discrete **dv,
1574 int32 numdvs,
1575 struct boolean_values *bval)
1576 {
1577 struct dis_discrete *cur_dis;
1578 int32 d;
1579 int32 orig_value;
1580
1581 for (d=0; d<numdvs; d++) {
1582 cur_dis = dv[d];
1583 if(dis_inwhen(cur_dis) && dis_boolean(cur_dis)) {
1584 orig_value = bval->cur_val[d];
1585 if (orig_value != dis_value(cur_dis)) {
1586 dis_set_val_modified(cur_dis,TRUE);
1587 }
1588 }
1589 }
1590 }
1591
1592
1593 /*
1594 * Analysis of Boundaries
1595 * ----------------------------
1596 */
1597
1598 /*
1599 * Evaluates the current status (satisfied? , at zero?) of a boundary
1600 * and sets its flags accordingly. At the same time, it updates the
1601 * residual of the relation included in the boundary (see
1602 * bndman_calc_satisfied).
1603 */
1604
1605 static void update_boundaries(slv_system_t server, SlvClientToken asys)
1606 {
1607 struct bnd_boundary **bp;
1608 int32 numbnds, ind, value;
1609 slv9_system_t sys;
1610
1611 sys = SLV9(asys);
1612 check_system(sys);
1613
1614 if( sys->blist == NULL ) {
1615 FPRINTF(ASCERR,"ERROR: (slv9) update_boundaries.\n");
1616 FPRINTF(ASCERR," Boundary list was never set.\n");
1617 return;
1618 }
1619
1620 bp = sys->blist;
1621 numbnds = slv_get_num_solvers_bnds(server);
1622
1623 for( ind = 0; ind < numbnds; ++ind ) {
1624 value = bnd_status_cur(bp[ind]);
1625 bnd_set_pre_status(bp[ind],value);
1626 value = bndman_calc_satisfied(bp[ind]);
1627 bnd_set_cur_status(bp[ind],value);
1628 if ( (bnd_status_cur(bp[ind]) != bnd_status_pre(bp[ind])) &&
1629 bnd_kind(bp[ind]) == e_bnd_rel && !bnd_at_zero(bp[ind])) {
1630 bnd_set_crossed(bp[ind],TRUE);
1631 } else {
1632 bnd_set_crossed(bp[ind],FALSE);
1633 }
1634 if (bnd_kind(bp[ind]) == e_bnd_rel) {
1635 value = bndman_calc_at_zero(bp[ind]);
1636 bnd_set_at_zero(bp[ind],value);
1637 } else {
1638 bnd_set_at_zero(bp[ind],FALSE);
1639 }
1640 }
1641 }
1642
1643
1644 /*
1645 * Look for some boundary with the CROSSED bit active. If this boundary
1646 * is used in some logical relation, the function returns 1, else returns 0
1647 */
1648 static int32 some_boundaries_crossed(slv_system_t server, SlvClientToken asys)
1649 {
1650 struct bnd_boundary **bp, *cur_bnd;
1651 int32 numbnds, ind;
1652 slv9_system_t sys;
1653
1654 sys = SLV9(asys);
1655 check_system(sys);
1656
1657 if( sys->blist == NULL ) {
1658 FPRINTF(ASCERR,"ERROR: (slv9) some_boundaries_crossed\n");
1659 FPRINTF(ASCERR," Boundary list was never set.\n");
1660 return 0;
1661 }
1662
1663 bp = sys->blist;
1664 numbnds = slv_get_num_solvers_bnds(server);
1665 for( ind = 0; ind < numbnds; ++ind ) {
1666 cur_bnd = bp[ind];
1667 if (bnd_crossed(cur_bnd) && bnd_in_logrel(cur_bnd)) {
1668 return 1;
1669 }
1670 }
1671 return 0;
1672 }
1673
1674 /*
1675 * Look for some boundary with the AT_ZERO bit active.If this boundary
1676 * is used in some logical relation, the function returns 1, else returns 0
1677 */
1678 static int32 some_boundaries_at_zero(slv_system_t server, SlvClientToken asys)
1679 {
1680 struct bnd_boundary **bp, *cur_bnd;
1681 int32 numbnds, ind;
1682 slv9_system_t sys;
1683
1684 sys = SLV9(asys);
1685 check_system(sys);
1686
1687 if( sys->blist == NULL ) {
1688 FPRINTF(ASCERR,"ERROR: (slv9) some_boundaries_at_zero\n");
1689 FPRINTF(ASCERR," Boundary list was never set.\n");
1690 return 0;
1691 }
1692
1693 bp = sys->blist;
1694 numbnds = slv_get_num_solvers_bnds(server);
1695 for( ind = 0; ind < numbnds; ++ind ) {
1696 cur_bnd = bp[ind];
1697 if (bnd_at_zero(cur_bnd) && bnd_in_logrel(cur_bnd)) {
1698 return 1;
1699 }
1700 }
1701 return 0;
1702 }
1703
1704 /*
1705 * Perform the combinatorial perturbation of the boundaries which are
1706 * at their zero. That means: We are going to perform a combinatorial
1707 * search, changing the truth value of a SATISFIED term (for the
1708 * specified boundaries) ON and OFF, and finding the boolean variables
1709 * affected for those changes in value of the SATISFIED terms.
1710 */
1711 static void do_perturbation_combinations(slv_system_t server,
1712 struct boolean_values *bval,
1713 struct bnd_boundary **bp,
1714 struct dis_discrete **dv,
1715 int32 numdvs,int32 *bndatzero,
1716 int32 ind, int32 numbndf)
1717 {
1718 slv_status_t status;
1719 int32 indpo;
1720
1721 if(ind<(numbndf-1)) {
1722 indpo = ind + 1;
1723 bnd_set_perturb(bp[bndatzero[ind]],TRUE);
1724 do_perturbation_combinations(server,bval,bp,dv,numdvs,
1725 bndatzero,indpo,numbndf);
1726 bnd_set_perturb(bp[bndatzero[ind]],FALSE);
1727 do_perturbation_combinations(server,bval,bp,dv,numdvs,
1728 bndatzero,indpo,numbndf);
1729 } else {
1730 if (ind < numbndf) {
1731 bnd_set_perturb(bp[bndatzero[ind]],TRUE);
1732 solve_logical_relations(server);
1733 slv_get_status(server,&status);
1734 if(!status.converged) {
1735 FPRINTF(ASCERR,"WARNING: \n");
1736 FPRINTF(ASCERR,"(slv9) do_perturbation_combinations\n");
1737 FPRINTF(ASCERR," Not convergence in logical solver \n");
1738 } else {
1739 search_for_modified_dvars(dv,numdvs,bval);
1740 }
1741 bnd_set_perturb(bp[bndatzero[ind]],FALSE);
1742 solve_logical_relations(server);
1743 slv_get_status(server,&status);
1744 if(!status.converged) {
1745 FPRINTF(ASCERR,"WARNING: \n");
1746 FPRINTF(ASCERR,"(slv9) do_perturbation_combinations\n");
1747 FPRINTF(ASCERR," Not convergence in logical solver \n");
1748 } else {
1749 search_for_modified_dvars(dv,numdvs,bval);
1750 }
1751 } else {
1752 FPRINTF(ASCERR,"ERROR: (slv9) do_perturbation_combinations\n");
1753 FPRINTF(ASCERR," Wrong boundary index as argument\n");
1754 }
1755 }
1756 return;
1757 }
1758
1759
1760 /*
1761 * Perform the combinatorial search of the subregions. That means:
1762 * We perform a combinatorial search, changing the value of the
1763 * discrete variables (given in disvars) TRUE and FALSE, and
1764 * finding which cases (in the WHENs) applies for each of the
1765 * combinations.
1766 */
1767 static void do_dvar_values_combinations(struct gl_list_t *disvars,
1768 struct matching_cases *cases,
1769 int numdvf, int d,
1770 int *pos_cases)
1771 {
1772 struct dis_discrete *cur_dis;
1773 int32 *cur_cases;
1774 int32 ncases, dpo;
1775
1776 if(d < numdvf) {
1777 dpo = d + 1;
1778 cur_dis = (struct dis_discrete *)(gl_fetch(disvars,d));
1779 dis_set_boolean_value(cur_dis,TRUE);
1780 do_dvar_values_combinations(disvars,cases,numdvf,dpo,pos_cases);
1781 dis_set_boolean_value(cur_dis,FALSE);
1782 do_dvar_values_combinations(disvars,cases,numdvf,dpo,pos_cases);
1783 } else {
1784 if (d == numdvf) {
1785 cur_dis = (struct dis_discrete *)(gl_fetch(disvars,d));
1786 dis_set_boolean_value(cur_dis,TRUE);
1787 cur_cases = cases_matching(disvars,&ncases);
1788 cases[(*pos_cases)].case_list = cur_cases;
1789 cases[(*pos_cases)].ncases = ncases;
1790 cases[(*pos_cases)].diff_subregion = 1;
1791 (*pos_cases)++;
1792 dis_set_boolean_value(cur_dis,FALSE);
1793 cur_cases = cases_matching(disvars,&ncases);
1794 cases[(*pos_cases)].case_list = cur_cases;
1795 cases[(*pos_cases)].ncases = ncases;
1796 cases[(*pos_cases)].diff_subregion = 1;
1797 (*pos_cases)++;
1798 } else {
1799 FPRINTF(ASCERR,"ERROR: (slv9) do_dvar_values_combinations\n");
1800 FPRINTF(ASCERR," Wrong discrete var index as argument\n");
1801 }
1802 }
1803 return;
1804 }
1805
1806
1807 /*
1808 * Orders of the elements of each array of cases,
1809 * so that we can compare them easier.
1810 */
1811 static void order_case(int32 *case_list, int32 *newcaselist, int ncases)
1812 {
1813 int32 cur_case,pos=0,tmp_num,c,ind;
1814
1815 for (c=1; c<=ncases; c++) {
1816 tmp_num = 0;
1817 for (ind=1; ind<=ncases; ind++) {
1818 cur_case = case_list[ind];
1819 if (tmp_num < cur_case) {
1820 pos = ind;
1821 tmp_num = cur_case;
1822 }
1823 }
1824 case_list[pos] = 0;
1825 newcaselist[ncases-c] = tmp_num;
1826 }
1827 }
1828
1829
1830
1831 /*
1832 * Calls for the ordering of the elements of each array of cases,
1833 * so that we can compare them easier.
1834 */
1835 static void order_cases(struct matching_cases *cases,int pos_cases)
1836 {
1837 int32 *caselist;
1838 int32 cur_ncase,c;
1839 int32 *newcaselist;
1840
1841 for (c=0; c<pos_cases;c++) {
1842 caselist = cases[c].case_list;
1843 cur_ncase = cases[c].ncases;
1844 if (cur_ncase > 1) {
1845 newcaselist = create_array(cur_ncase,int32);
1846 order_case(caselist,newcaselist,cur_ncase);
1847 cases[c].case_list = newcaselist;
1848 destroy_array(caselist);
1849 } else {
1850 if (cur_ncase == 1) {
1851 newcaselist = create_array(1,int32);
1852 newcaselist[0] = caselist[1];
1853 cases[c].case_list = newcaselist;
1854 destroy_array(caselist);
1855 }
1856 }
1857 }
1858
1859 }
1860
1861
1862
1863 /*
1864 * Compare two arrays of cases (integer numbers). It returns 1 if they are
1865 * equal, else it returns 0.
1866 */
1867 static int32 compare_case(int32 *cur_set, int32 *comp_set, int cur_ncases)
1868 {
1869 int32 cur_case, comp_case, ind;
1870
1871 for (ind=0; ind<cur_ncases; ind++) {
1872 cur_case = cur_set[ind];
1873 comp_case = comp_set[ind];
1874 if (cur_case != comp_case) {
1875 return 0;
1876 }
1877 }
1878 return 1;
1879 }
1880
1881
1882 /*
1883 * Compare the arrays of cases so that we can find the number of
1884 * different alternatives (subregions)
1885 */
1886 static void compare_cases(struct matching_cases *cases,int pos_cases)
1887 {
1888 int32 *cur_set, *comp_set, cur_ncases, comp_ncases;
1889 int32 c,d;
1890
1891 for (c=0; c<pos_cases; c++) {
1892 cur_set = cases[c].case_list;
1893 cur_ncases = cases[c].ncases;
1894 if (cur_ncases == 0) {
1895 cases[c].diff_subregion = 0;
1896 continue;
1897 }
1898 for(d=0; d<c; d++) {
1899 comp_set = cases[d].case_list;
1900 comp_ncases = cases[d].ncases;
1901 if (cur_ncases != comp_ncases) {
1902 continue;
1903 } else {
1904 if (compare_case(cur_set,comp_set,cur_ncases)) {
1905 cases[c].diff_subregion = 0;
1906 break;
1907 }
1908 }
1909 }
1910 }
1911 }
1912
1913
1914 /*
1915 * Finds if my current point lies at a "real" boundary. By "real"
1916 * I mean a boundary which really causes a change in the
1917 * configuration. It returns 0 if the boundary at zero does not
1918 * affect the configuration. If the configuration is affected,
1919 * this function will find the number of subregions existing
1920 * for the current point as well as the cases (in WHENs) corresponding
1921 * to each of the subregions. At the end, the number of subregions is
1922 * n_subregions and the cases applying for each of them is stored
1923 * in the structure subregions.
1924 */
1925 static int32 at_a_boundary(slv_system_t server, SlvClientToken asys,
1926 int32 *n_subregions,
1927 struct matching_cases **subregions,
1928 int32 *cur_subregion,
1929 struct gl_list_t *disvars)
1930 {
1931 slv9_system_t sys;
1932 struct bnd_boundary **bp, *cur_bnd;
1933 struct dis_discrete **dv, *cur_dis;
1934 struct boolean_values bval;
1935 dis_filter_t dfilter;
1936 bnd_filter_t bfilter;
1937 struct matching_cases *cases;
1938 int32 *bndatzero;
1939 int32 *dvarmodified;
1940 int32 *cur_cases;
1941 int32 *caselist, *newcaselist;
1942 int32 numbnds, numbndf, b, ind;
1943 int32 numdvs, numdvf, d;
1944 int32 cur_ncases, assign_cur_sub;
1945 int32 pos_cases, comb;
1946 char *param;
1947 union param_value u;
1948
1949 sys = SLV9(asys);
1950 check_system(sys);
1951
1952 if( sys->blist == NULL ) {
1953 FPRINTF(ASCERR,"ERROR: (slv9) at_a_boundary\n");
1954 FPRINTF(ASCERR," Boundary list was never set.\n");
1955 return 0;
1956 }
1957
1958 if( sys->dvlist == NULL ) {
1959 FPRINTF(ASCERR,"ERROR: (slv9) at_a_boundary\n");
1960 FPRINTF(ASCERR," Discrete Variable list was never set.\n");
1961 return 0;
1962 }
1963
1964 if (!some_boundaries_at_zero(server,asys)) {
1965 return 0;
1966 }
1967
1968 bp = sys->blist;
1969 numbnds = slv_get_num_solvers_bnds(server);
1970 bfilter.matchbits = (BND_AT_ZERO);
1971 bfilter.matchvalue = (BND_AT_ZERO);
1972 numbndf = slv_count_solvers_bnds(server,&bfilter);
1973 bndatzero = create_array(numbndf,int32);
1974 ind = 0;
1975 for (b=0; b<numbnds; b++) {
1976 cur_bnd = bp[b];
1977 bnd_set_perturb(cur_bnd,FALSE);
1978 if (bnd_at_zero(cur_bnd)) {
1979 #if SHOW_BOUNDARY_ANALYSIS_DETAILS
1980 FPRINTF(ASCERR,"boundary at zero = %d\n",b);
1981 #endif /* SHOW_BOUNDARY_ANALYSIS_DETAILS */
1982 bndatzero[ind] = b;
1983 ind++;
1984 }
1985 }
1986
1987 dv = sys->dvlist;
1988 numdvs = slv_get_num_solvers_dvars(server);
1989 bval.cur_val = create_array(numdvs,int32);
1990 bval.pre_val = create_array(numdvs,int32);
1991
1992 for (d=0; d<numdvs; d++) {
1993 cur_dis = dv[d];
1994 dis_set_val_modified(cur_dis,FALSE);
1995 bval.cur_val[d] = dis_value(cur_dis);
1996 bval.pre_val[d] = dis_previous_value(cur_dis);
1997 }
1998
1999 #if SHOW_BOUNDARY_ANALYSIS_DETAILS
2000 FPRINTF(ASCERR,"Executing combinatorial perturbation of boundaries\n");
2001 #endif /* SHOW_BOUNDARY_ANALYSIS_DETAILS */
2002
2003 /*
2004 * Setting the value of the perturbation mode flag in the logical solver
2005 * to 1.
2006 * PERTURB_BOUNDARY is a boolean parameter of the logical solver
2007 * LRSlv. This parameter tells the solver whether it should change
2008 * the truth value of the SATISFIED terms or not (only for the
2009 * boundaries specified). This trick is important while finding
2010 * the number of subregions around a/several boundary(ies).
2011 */
2012 param = "perturbboundaries";
2013 u.b = 1;
2014 set_param_in_solver(server,LOGICAL_SOLVER,bool_parm,param,&u);
2015
2016 ind = 0;
2017 do_perturbation_combinations(server,&(bval),bp,dv,numdvs,
2018 bndatzero,ind,numbndf);
2019 /*
2020 * Setting the value of the perturbation mode flag in the logical solver
2021 * to 0.
2022 */
2023 u.b = 0;
2024 set_param_in_solver(server,LOGICAL_SOLVER,bool_parm,param,&u);
2025
2026 destroy_array(bndatzero);
2027
2028 dfilter.matchbits = (DIS_VAL_MODIFIED);
2029 dfilter.matchvalue = (DIS_VAL_MODIFIED);
2030 numdvf = slv_count_solvers_dvars(server,&dfilter);
2031
2032 if (numdvf == 0) {
2033 FPRINTF(ASCERR,"Not really at a boundary\n");
2034 for (d=0; d<numdvs; d++) {
2035 cur_dis = dv[d];
2036 dis_set_boolean_value(cur_dis,bval.cur_val[d]);
2037 dis_set_value(cur_dis,bval.cur_val[d]);
2038 dis_set_previous_value(cur_dis,bval.pre_val[d]);
2039 }
2040 destroy_array(bval.cur_val);
2041 destroy_array(bval.pre_val);
2042 return 0;
2043 }
2044
2045 dvarmodified = create_array(numdvf,int32);
2046 ind = 0;
2047 for (d=0; d<numdvs; d++) {
2048 cur_dis = dv[d];
2049 if(dis_val_modified(cur_dis)) {
2050 dvarmodified[ind] = d;
2051 gl_append_ptr(disvars,cur_dis);
2052 dis_set_val_modified(cur_dis,FALSE);
2053 ind++;
2054 }
2055 }
2056
2057 for (d=0; d<numdvs; d++) {
2058 cur_dis = dv[d];
2059 dis_set_boolean_value(cur_dis,bval.cur_val[d]);
2060 dis_set_value(cur_dis,bval.cur_val[d]);
2061 dis_set_previous_value(cur_dis,bval.pre_val[d]);
2062 }
2063
2064 pos_cases = 1;
2065 for (d = 1; d<=numdvf; d++) {
2066 pos_cases = pos_cases * 2;
2067 }
2068
2069 cases = (struct matching_cases *)
2070 (ascmalloc((pos_cases)*sizeof(struct matching_cases)));
2071
2072 #if SHOW_BOUNDARY_ANALYSIS_DETAILS
2073 FPRINTF(ASCERR,"Executing combinatorial search for subregions \n");
2074 #endif /* SHOW_BOUNDARY_ANALYSIS_DETAILS */
2075
2076 d = 1;
2077 comb = 0;
2078 do_dvar_values_combinations(disvars,cases,numdvf,d,&(comb));
2079
2080 order_cases(cases,pos_cases);
2081 compare_cases(cases,pos_cases);
2082
2083 (*n_subregions) = 0;
2084 for(comb=0; comb<pos_cases;comb++) {
2085 if(cases[comb].diff_subregion) {
2086 (*n_subregions)++;
2087 }
2088 }
2089
2090 if ((*n_subregions)==0) {
2091 FPRINTF(ASCERR,"ERROR: at least one subregion must be found\n");
2092 for (d=0; d<numdvs; d++) {
2093 cur_dis = dv[d];
2094 dis_set_boolean_value(cur_dis,bval.cur_val[d]);
2095 dis_set_value(cur_dis,bval.cur_val[d]);
2096 dis_set_previous_value(cur_dis,bval.pre_val[d]);
2097 }
2098 destroy_array(bval.cur_val);
2099 destroy_array(bval.pre_val);
2100 for(comb=0; comb<pos_cases;comb++) {
2101 destroy_array(cases[comb].case_list);
2102 }
2103 destroy_array(cases);
2104 return 0;
2105 }
2106
2107 if ((*n_subregions)==1) {
2108 FPRINTF(ASCERR,"Not really at a boundary\n");
2109 for (d=0; d<numdvs; d++) {
2110 cur_dis = dv[d];
2111 dis_set_boolean_value(cur_dis,bval.cur_val[d]);
2112 dis_set_value(cur_dis,bval.cur_val[d]);
2113 dis_set_previous_value(cur_dis,bval.pre_val[d]);
2114 }
2115 destroy_array(bval.cur_val);
2116 destroy_array(bval.pre_val);
2117 for(comb=0; comb<pos_cases;comb++) {
2118 destroy_array(cases[comb].case_list);
2119 }
2120 destroy_array(cases);
2121 return 0;
2122 }
2123
2124 if ((*n_subregions) > 0) {
2125 (*subregions) = (struct matching_cases *)
2126 (ascmalloc(((*n_subregions))*sizeof(struct matching_cases)));
2127 (*n_subregions) = 0;
2128 for(comb=0; comb<pos_cases;comb++) {
2129 if(cases[comb].diff_subregion) {
2130 (*subregions)[(*n_subregions)].case_list = cases[comb].case_list;
2131 cases[comb].case_list = NULL;
2132 (*subregions)[(*n_subregions)].ncases = cases[comb].ncases;
2133 cases[comb].ncases = 0;
2134 (*subregions)[(*n_subregions)].diff_subregion = 1;
2135 (*n_subregions)++;
2136 }
2137 }
2138 }
2139
2140 for(comb=0; comb<pos_cases;comb++) {
2141 destroy_array(cases[comb].case_list);
2142 }
2143 destroy_array(cases);
2144
2145
2146 assign_cur_sub = 0;
2147 /*
2148 * Finding the subregion corresponding to the "original" configuration
2149 */
2150 for (d=0; d<numdvs; d++) {
2151 cur_dis = dv[d];
2152 dis_set_boolean_value(cur_dis,bval.cur_val[d]);
2153 dis_set_value(cur_dis,bval.cur_val[d]);
2154 dis_set_previous_value(cur_dis,bval.pre_val[d]);
2155 }
2156 cur_cases = cases_matching(disvars,&cur_ncases);
2157 caselist = cur_cases;
2158 if (cur_ncases > 1) {
2159 newcaselist = create_array(cur_ncases,int32);
2160 order_case(caselist,newcaselist,cur_ncases);
2161 cur_cases = newcaselist;
2162 destroy_array(caselist);
2163 } else {
2164 if (cur_ncases == 1) {
2165 newcaselist = create_array(1,int32);
2166 newcaselist[0] = caselist[1];
2167 cur_cases = newcaselist;
2168 destroy_array(caselist);
2169 }
2170 }
2171 for(comb=0; comb<(*n_subregions);comb++) {
2172 if ((*subregions)[comb].ncases == cur_ncases) {
2173 if(compare_case((*subregions)[comb].case_list,cur_cases,cur_ncases)) {
2174 (*cur_subregion) = comb;
2175 assign_cur_sub = 1;
2176 break;
2177 }
2178 }
2179 }
2180
2181 if (!assign_cur_sub) {
2182 FPRINTF(ASCERR,"PANIC: original configuration not found\n");
2183 }
2184
2185 destroy_array(cur_cases);
2186 destroy_array(dvarmodified);
2187 destroy_array(bval.cur_val);
2188 destroy_array(bval.pre_val);
2189 return 1;
2190 }
2191
2192
2193 /*
2194 * If some boundary(ies) has been crossed in the iterative scheme,
2195 * this function finds the boundary crossed (the first one, if many).
2196 * It returns the factor (less than 1) by which the step length has
2197 * to be multiplied son that the new point will lie precisely
2198 * at that boundary. This factor is found by using the method of
2199 * bisection
2200 */
2201 static real64 return_to_first_boundary(slv_system_t server,
2202 SlvClientToken asys,
2203 struct real_values *rvalues,
2204 var_filter_t *vfilter)
2205 {
2206 slv9_system_t sys;
2207 struct bnd_boundary **bp, *cur_bnd;
2208 struct var_variable **incidences, **bnd_incidences;
2209 struct var_variable *cur_var;
2210 bnd_filter_t bfilter;
2211 struct boolean_values bval;
2212 real64 factor=0.0,fup,flo,newvalue;
2213 int32 *bndcrossed;
2214 int32 *inc_vars;
2215 int32 count,n_incidences,inc,conv_flag,still_crossed;
2216 int32 numbnds,numbndf,b,ind;
2217 int32 iter,n_iterations;
2218 FILE *lif;
2219
2220 sys = SLV9(asys);
2221 check_system(sys);
2222 lif = LIF(sys);
2223
2224 if( sys->blist == NULL ) {
2225 FPRINTF(ASCERR,"ERROR: (slv9) return_to_first_boundary\n");
2226 FPRINTF(ASCERR," Boundary list was never set.\n");
2227 return 1.0;
2228 }
2229
2230 if (!some_boundaries_crossed(server,asys)) {
2231 return 1.0;
2232 }
2233
2234 bp = sys->blist;
2235 numbnds = slv_get_num_solvers_bnds(server);
2236 bfilter.matchbits = (BND_CROSSED);
2237 bfilter.matchvalue = (BND_CROSSED);
2238 numbndf = slv_count_solvers_bnds(server,&bfilter);
2239 bndcrossed = create_array(numbndf,int32);
2240 bval.cur_val = create_array(numbndf,int32);
2241 bval.pre_val = create_array(numbndf,int32);
2242 ind = 0;
2243 for (b=0; b<numbnds; b++) {
2244 cur_bnd = bp[b];
2245 if (bnd_crossed(cur_bnd)) {
2246 bndcrossed[ind] = b;
2247 bval.cur_val[ind] = bnd_status_cur(cur_bnd);
2248 bval.pre_val[ind] = bnd_status_pre(cur_bnd);
2249 ind++;
2250 }
2251 }
2252
2253 count = 0;
2254 for (b=0; b<numbndf; b++) {
2255 cur_bnd = bp[bndcrossed[b]];
2256 n_incidences = bnd_n_real_incidences(cur_bnd);
2257 count = count + n_incidences;
2258 }
2259
2260 incidences = (struct var_variable **)
2261 ( ascmalloc((count)*sizeof(struct var_variable *)));
2262 inc_vars = create_array(count,int32);
2263 count = 0;
2264 for (b=0; b<numbndf; b++) {
2265 cur_bnd = bp[bndcrossed[b]];
2266 bnd_incidences = bnd_real_incidence(cur_bnd);
2267 n_incidences = bnd_n_real_incidences(cur_bnd);
2268 #if SHOW_BOUNDARY_ANALYSIS_DETAILS
2269 FPRINTF(lif,"boundary crossed = %d\n",bndcrossed[b]);
2270 FPRINTF(lif,"previous boundary status = %d\n",bval.pre_val[b]);
2271 FPRINTF(lif,"current boundary status = %d\n",bval.cur_val[b]);
2272 #endif /* SHOW_BOUNDARY_ANALYSIS_DETAILS */
2273 for (inc=0; inc<n_incidences; inc++) {
2274 incidences[count] = bnd_incidences[inc];
2275 inc_vars[count] = var_mindex(incidences[count]);
2276 count++;
2277 }
2278 }
2279
2280 /* bisection to find first boundary crossed */
2281 fup = 1.0;
2282 flo = 0.0;
2283 conv_flag = 0;
2284 iter = 0;
2285
2286 /*
2287 * Maximum number of bisection iterations. This must be modified
2288 * so that it becomes a parameter to be defined by the user
2289 */
2290 n_iterations = ITER_BIS_LIMIT;
2291
2292 #if SHOW_BOUNDARY_ANALYSIS_DETAILS
2293 for (inc=0; inc<count; inc++) {
2294 cur_var = incidences[inc];
2295 if(var_apply_filter(cur_var,vfilter)) {
2296 FPRINTF(lif,"Variable ");
2297 print_var_name(lif,sys,cur_var); PUTC('\n',lif);
2298 FPRINTF(lif,
2299 "previous value = %f\n",rvalues->pre_values[inc_vars[inc]]);
2300 FPRINTF(lif,"current value = %f\n",rvalues->cur_values[inc_vars[inc]]);
2301 }
2302 }
2303 #endif /* SHOW_BOUNDARY_ANALYSIS_DETAILS */
2304
2305 while (conv_flag == 0) {
2306 iter++;
2307 if (iter>n_iterations) {
2308 FPRINTF(ASCERR,"ERROR: (slv9) return_to_first_boundary\n");
2309 FPRINTF(ASCERR,"Could not find the first boundary crossed \n");
2310 FPRINTF(ASCERR,"Returning the last factor calculated\n");
2311 break;
2312 }
2313 still_crossed = 0;
2314 factor = ( fup + flo ) / 2.0;
2315 #if SHOW_BISECTION_DETAILS
2316 FPRINTF(lif,"fup = %f\n",fup);
2317 FPRINTF(lif,"flo = %f\n",flo);
2318 FPRINTF(lif,"factor = %f\n",factor);
2319 #endif /* SHOW_BISECTION_DETAILS */
2320 for (inc=0; inc<count; inc++) {
2321 cur_var = incidences[inc];
2322 if(var_apply_filter(cur_var,vfilter)) {
2323 newvalue = rvalues->pre_values[inc_vars[inc]] + factor *
2324 ( rvalues->cur_values[inc_vars[inc]] -
2325 rvalues->pre_values[inc_vars[inc]] );
2326 var_set_value(cur_var,newvalue);
2327 #if SHOW_BISECTION_DETAILS
2328 FPRINTF(lif,"Variable ");
2329 print_var_name(lif,sys,cur_var); PUTC('\n',lif);
2330 FPRINTF(lif,"value after factor = %f\n",newvalue);
2331 #endif /* SHOW_BISECTION_DETAILS */
2332 }
2333 }
2334
2335 update_boundaries(server,asys);
2336 for (b=0; b<numbndf; b++) {
2337 cur_bnd = bp[bndcrossed[b]];
2338 #if SHOW_BISECTION_DETAILS
2339 FPRINTF(lif,"previous status = %d\n", bval.pre_val[b]);
2340 FPRINTF(lif,"status aftert factor = %d\n",bnd_status_cur(cur_bnd));
2341 #endif /* SHOW_BISECTION_DETAILS */
2342 if ( bnd_status_cur(cur_bnd) != bval.pre_val[b] ) {
2343 still_crossed = 1;
2344 }
2345 }
2346 #if SHOW_BISECTION_DETAILS
2347 FPRINTF(lif,"still_crossed = %d\n",still_crossed);
2348 #endif /* SHOW_BISECTION_DETAILS */
2349 if (still_crossed) {
2350 fup = factor;
2351 } else {
2352 flo = factor;
2353 for (b=0; b<numbndf; b++) {
2354 cur_bnd = bp[bndcrossed[b]];
2355 bnd_set_pre_status(cur_bnd,bval.pre_val[b]);
2356 bnd_set_cur_status(cur_bnd,bval.pre_val[b]);
2357 if (bnd_at_zero(cur_bnd)) {
2358 #if SHOW_BOUNDARY_ANALYSIS_DETAILS
2359 FPRINTF(ASCERR,"boundary at zero = %d\n",bndcrossed[b]);
2360 FPRINTF(lif,"factor = %f\n",factor);
2361 for (inc=0; inc<count; inc++) {
2362 cur_var = incidences[inc];
2363 if(var_apply_filter(cur_var,vfilter)) {
2364 FPRINTF(lif,"Variable ");
2365 print_var_name(lif,sys,cur_var); PUTC('\n',lif);
2366 FPRINTF(lif,"value after factor = %f\n",var_value(cur_var));
2367 }
2368 }
2369 #endif /* SHOW_BOUNDARY_ANALYSIS_DETAILS */
2370 conv_flag = 1;
2371 }
2372 }
2373 }
2374 }
2375 destroy_array(bndcrossed);
2376 destroy_array(inc_vars);
2377 destroy_array(bval.cur_val);
2378 destroy_array(bval.pre_val);
2379 destroy_array(incidences);
2380
2381 return factor;
2382 }
2383
2384
2385 /*
2386 * Storing values of real variables.
2387 * ---------------------------------
2388 *
2389 * We use the master list of variables since its order does not change
2390 * and it is given by the master index. We do not touch the master list,
2391 * we only use its order.
2392 */
2393
2394 /*
2395 * Store the values of the var_variables before a newton-like iteration
2396 */
2397 static void store_real_pre_values(slv_system_t server,
2398 struct real_values *rvalues)
2399 {
2400 struct var_variable **master;
2401 struct var_variable *var;
2402 int v, vlen;
2403
2404 master = slv_get_master_var_list(server);
2405 vlen = slv_get_num_master_vars(server);
2406
2407 rvalues->pre_values = create_array(vlen,real64);
2408
2409 for (v=0; v<vlen; v++) {
2410 var = master[v];
2411 rvalues->pre_values[v] = var_value(var);
2412 }
2413 }
2414
2415 /*
2416 * Store the values of the var_variables after a newton-like iteration
2417 */
2418 static void store_real_cur_values(slv_system_t server,
2419 struct real_values *rvalues)
2420 {
2421 struct var_variable **master;
2422 struct var_variable *var;
2423 int v, vlen;
2424
2425 master = slv_get_master_var_list(server);
2426 vlen = slv_get_num_master_vars(server);
2427
2428 rvalues->cur_values = create_array(vlen,real64);
2429
2430 for (v=0; v<vlen; v++) {
2431 var = master[v];
2432 rvalues->cur_values[v] = var_value(var);
2433 }
2434 }
2435
2436 /*
2437 * After the length of the step has been modified so that the current point
2438 * lies at a boundary, the values of all the variables is updated so that
2439 * they all reduce the length of their step by the same factor.
2440 */
2441 static void update_real_var_values(slv_system_t server,
2442 struct real_values *rvalues,
2443 var_filter_t *vfilter, real64 factor)
2444 {
2445 struct var_variable **master;
2446 struct var_variable *var;
2447 real64 newvalue;
2448 int v, vlen;
2449
2450 master = slv_get_master_var_list(server);
2451 vlen = slv_get_num_master_vars(server);
2452
2453 for (v=0; v<vlen; v++) {
2454 var = master[v];
2455 if(var_apply_filter(var,vfilter)) {
2456 newvalue = rvalues->pre_values[v] +
2457 factor * (rvalues->cur_values[v] - rvalues->pre_values[v]);
2458 var_set_value(var,newvalue);
2459 }
2460 }
2461 destroy_array(rvalues->cur_values);
2462 destroy_array(rvalues->pre_values);
2463 }
2464
2465
2466 /*
2467 * Set the flagbit NONBASIC for all the variables in the list
2468 * to the value passed as argument
2469 */
2470 static void set_nonbasic_status_in_var_list(slv_system_t server,
2471 uint32 value)
2472 {
2473 struct var_variable **master;
2474 struct var_variable *var;
2475 int v, vlen;
2476
2477 master = slv_get_master_var_list(server);
2478 vlen = slv_get_num_master_vars(server);
2479
2480 for (v=0; v<vlen; v++) {
2481 var = master[v];
2482 var_set_nonbasic(var,value);
2483 }
2484 }
2485
2486
2487 /*
2488 * After the length of the step has been modified so that the current point
2489 * lies at a boundary, the residuals of the equations are updated.
2490 */
2491 static void update_relations_residuals(slv_system_t server)
2492 {
2493 struct rel_relation **master;
2494 struct rel_relation *rel;
2495 rel_filter_t rfilter;
2496 real64 resid;
2497 int32 r, rlen, status;
2498
2499 master = slv_get_master_rel_list(server);
2500 rlen = slv_get_num_master_rels(server);
2501
2502 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY);
2503 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY);
2504
2505 #ifdef ASC_SIGNAL_TRAPS
2506 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
2507 #endif
2508
2509 for (r=0; r<rlen; r++) {
2510 rel = master[r];
2511 if(rel_apply_filter(rel,&rfilter)) {
2512 resid = relman_eval(rel,&status,1);
2513 }
2514 }
2515 #ifdef ASC_SIGNAL_TRAPS
2516 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
2517 #endif
2518
2519 }
2520
2521
2522 #ifdef ASC_WITH_CONOPT
2523 /*------------------------------------------------------------------------------
2524 CALLBACK ROUTINES FOR CONOPT
2525 */
2526
2527 /*
2528 * COIRMS Based on the information provided in Coispz, CONOPT will
2529 * allocate the number of vectors into which the user can define
2530 * the details of the model. The details of the model are defined
2531 * here.
2532 *
2533 * COIRMS(lower, curr, upper, vsta, type,rhs, fv, esta, colsta,
2534 * rowno, value, nlflag, n, m, nz, usrmem)
2535 *
2536 * lower - lower bounds on the variables
2537 * curr - intial values of the variables
2538 * upper - upper bounds on the variables
2539 * vsta - initial status of the variable(o nonbasic, 1 basic)
2540 * type - types of equations (0 equality,1 greater,2 less)
2541 * rhs - values of the right hand sides
2542 * fv - sum of the nonlinear terms in the initial point
2543 * esta - initial status of the slack in the constraint (nonbasic,basic)
2544 * colsta- start of column pointers
2545 * rowno - row or equation numbers of the nonzeros
2546 * value - values of the jacobian elements
2547 * nlflag- nonlinearity flags(0 nonzero constant,1 varying)
2548 * n - number of variables
2549 * m - number of constraints
2550 * nz - number of jacobian elements
2551 * usrmem- user memory defined by conopt
2552 */
2553 static int COI_CALL slv9_conopt_readmatrix(
2554 double *lower, double *curr, double *upper
2555 , int *vsta, int *type, double *rhs
2556 , int *esta, int *colsta, int *rowno
2557 , double *value, int *nlflag, int *n, int *m, int *nz
2558 , double *usrmem
2559 ){
2560 slv9_system_t sys;
2561 struct var_variable *var;
2562 struct var_variable **varlist;
2563 struct opt_matrix *coeff_matrix;
2564 static var_filter_t vfilter;
2565 real64 /*obj_val,*/ deriv;
2566 real64 nominal, up, low, uplow;
2567 int32 num_var, n_subregions, c, num_eqns;
2568 int32 numnz, eq;
2569 int32 count, totvar;
2570
2571 UNUSED_PARAMETER(vsta);
2572 UNUSED_PARAMETER(esta);
2573
2574 sys = (slv9_system_t)usrmem;
2575 n_subregions = sys->subregions;
2576 coeff_matrix = sys->coeff_matrix;
2577 num_var = (*n) - n_subregions;
2578 num_eqns = num_var + 1;
2579
2580 vfilter.matchbits = (VAR_ACTIVE_AT_BND | VAR_INCIDENT
2581 | VAR_SVAR | VAR_FIXED);
2582 vfilter.matchvalue = (VAR_ACTIVE_AT_BND | VAR_INCIDENT | VAR_SVAR);
2583 varlist = sys->mvlist;
2584 totvar = sys->mvtot;
2585
2586 /*
2587 * Variables: Current Value, lower value and upper value. Note that for
2588 * this problem the variables are the vector of steps dx. So we "invent"
2589 * approximations for the bounds based on the bounds of the real
2590 * variables of the problem. The use of the parameter ASC_INFINITY is
2591 * Hack to keep CONOPT from complaining for large bounds.
2592 */
2593
2594 count = 0;
2595 for (c=0; c<totvar; c++) {
2596 var = varlist[c];
2597 if (var_apply_filter(var,&vfilter)) {
2598 nominal = var_nominal(var);
2599 low = var_lower_bound(var);
2600 up = var_upper_bound(var);
2601 uplow = fabs( up - low);
2602
2603 if(-uplow > -CONOPT_BOUNDLIMIT){
2604 lower[count] = -uplow;
2605 }else{
2606 lower[count] = -0.5*CONOPT_BOUNDLIMIT;
2607 CONSOLE_DEBUG("REDUCING LOWER BOUND LIMIT FOR VAR %d TO %e",count,lower[count]);
2608 }
2609
2610 if(uplow < CONOPT_BOUNDLIMIT){
2611 upper[count] = uplow;
2612 }else{
2613 upper[count] = 0.5*CONOPT_BOUNDLIMIT;
2614 CONSOLE_DEBUG("REDUCING UPPER BOUND LIMIT FOR VAR %d TO %e",count,upper[count]);
2615 }
2616
2617 curr[count] = 0.5 * nominal;
2618 count++;
2619 }
2620 }
2621 /* alphas for each subregion */
2622 for (c=count; c<(*n); c++) {
2623 lower[c] = 0.0;
2624 upper[c] = 1.0;
2625 curr[c] = 1.0;
2626 }
2627
2628
2629 /*
2630 * vsta not in use since STATOK, ipsz[14], is zero
2631 */
2632
2633 /*
2634 * All the equations, but the last row (which is the objective), are
2635 * equalities.
2636 */
2637 for (c = 0; c < (*m); c++) {
2638 type[c] = 0;
2639 }
2640 type[(*m)-1] = 3;
2641
2642
2643 /*
2644 * RHS. It is zero for all the equations except the summation of
2645 * alphas, whose RHS is one.
2646 */
2647 for (c = 0; c < (*m); c++) {
2648 rhs[c] = 0;
2649 }
2650 rhs[(*m)-2] = 1.0;
2651
2652 #ifdef DISUSED_CONOPT_PARAMETER
2653 /*
2654 * fv =0 for all linear relations. For the objective is the two
2655 * norm
2656 */
2657 for (c = 0; c < (*m); c++) {
2658 fv[c] = 0;
2659 }
2660 obj_val = 0.0;
2661 for (c = 0; c<num_var; c++) {
2662 obj_val = obj_val + (curr[c] * curr[c]);
2663 }
2664 fv[(*m)-1] = obj_val;
2665 #endif
2666
2667 /*
2668 * esta not used since STATOK is zero
2669 */
2670
2671
2672 /*
2673 * For the following parameters, it is important ot see that:
2674 * 1) The values for the rows and nonzeros that conopt wants start
2675 * with 1, not with 0.
2676 * 2) The indeces of the arrays that we are using in the C side start
2677 * with 0.
2678 */
2679
2680 /*
2681 * colsta
2682 */
2683
2684 for (c=0; c<num_var; c++) {
2685 colsta[c] = 2 * c;
2686 }
2687
2688 for (c=num_var; c<(*n); c++) {
2689 colsta[c] = 2 * num_var + num_eqns * (c - num_var);
2690 }
2691
2692 colsta[*n] = *nz; /** @TODO check this */
2693
2694 /*
2695 rowno, value and nlflag can be done in same loop. The use of the
2696 parameter RTMAXJ is really a Hack to keep CONOPT from complaining
2697 about large derivatives
2698 */
2699
2700 numnz = 0;
2701 CONSOLE_DEBUG("NUM_VAR = %d",num_var);
2702 for (c=0; c<num_var; c++) {
2703 CONSOLE_DEBUG("C = %d",c);
2704 rowno[numnz] = c;
2705 nlflag[numnz] = 0;
2706 value[numnz] = -1;
2707 numnz++;
2708 CONSOLE_DEBUG("NUMNZ = %d",numnz);
2709 rowno[numnz] = *m - 1;
2710 nlflag[numnz] = 1;
2711 CONSOLE_DEBUG("FOUND A NONLINEAR NONZERO");
2712 numnz++;
2713 CONSOLE_DEBUG("NUMNZ = %d",numnz);
2714 }
2715
2716 for (c=num_var; c<(*n); c++) {
2717 numnz = 2 * num_var + num_eqns * (c - num_var);
2718 for(eq = 0; eq<num_eqns-1; eq++) {
2719 rowno[numnz] = eq;
2720 nlflag[numnz] = 0;
2721 deriv = -1.0 * (coeff_matrix->cols[c - num_var].element[eq]);
2722 if (deriv > RTMAXJ ) {
2723 deriv = 0.5 * RTMAXJ;
2724 } else {
2725 if (deriv < -RTMAXJ ) {
2726 deriv = -0.5*RTMAXJ;
2727 }
2728 }
2729 value[numnz] = deriv;
2730 numnz++;
2731 }
2732 rowno[numnz] = num_eqns - 1;
2733 nlflag[numnz] = 0;
2734 value[numnz] = 1.0;
2735 }
2736
2737 return 0;
2738 }
2739
2740 #if 0 /* not in API any more */
2741 /*
2742 * COIFBL Defines the nonlinearities of the model by returning
2743 * numerical values. It works on a block of rows during each call.
2744 * COIFBL( x, g, otn, nto, from, to, jac, stcl, rnum, cnum, nl, strw,
2745 * llen, indx, mode, errcnt, n, m, n1, m1, nz, usrmem)
2746 *
2747 * x - punt of evaluation provided by conopt
2748 * g - vector of function values
2749 * otn - old to new permutation vector
2750 * nto - new to old permutation vector
2751 * from - range in permutation
2752 * to - range in permutation
2753 * jac - vector of jacobian values.
2754 * The following are vectors defining the jacobian structure
2755 * stcl - start of column pointers
2756 * rnum - row numbers
2757 * cnum - column numbers
2758 * nl - nonlinearity flags
2759 * strw - start row pointers
2760 * llen - count of linear jacobian elements
2761 * indx - pointers from the row-wise representation
2762 * mode - indicator of mode of evaluation
2763 * errcnt- number of function evaluation errors
2764 * n - umber of variables
2765 * m - number of constraints
2766 * n1 - n+1
2767 * m1 - m+1
2768 * nz - number of jacobian elements
2769 * usrmem- user memory defined by conopt
2770 */
2771 static void slv9_coifbl(real64 *x, real64 *g, int32 *otn, int32 *nto,
2772 int32 *from, int32 *to, real64 *jac, int32 *stcl,
2773 int32 *rnum, int32 *cnum, int32 *nl, int32 *strw,
2774 int32 *llen, int32 *indx, int32 *mode, int32 *errcnt,
2775 int32 *n, int32 *m, int32 *n1, int32 *m1,
2776 int32 *nz, real64 *usrmem)
2777 {
2778 /* non defined for this solver */
2779
2780 /* stop gcc whining about unused parameter */
2781 (void)x; (void)g; (void)otn; (void)nto; (void)from; (void)to;
2782 (void)jac; (void)stcl; (void)rnum; (void)cnum; (void)nl; (void)strw;
2783 (void)llen; (void)indx; (void)mode; (void)errcnt; (void)n; (void)m;
2784 (void)n1; (void)m1; (void)nz;
2785 (void)usrmem;
2786
2787 return;
2788 }
2789 #endif
2790
2791 /*
2792 * COIFDE Defines the nonlinearities of the model by returning
2793 * numerical values. It works on one row or equation at a time
2794 * COIFDE(x, g, jac, rowno, jcnm, mode, errcnt, newpt, n, nj, usrmem)
2795 *
2796 * x - punt of evaluation provided by conopt
2797 * g - function value
2798 * jac - jacobian values
2799 * rowno - number of the row for which nonlinearities will be eval
2800 * jcnm - list of column number fon the NL nonzeros
2801 * mode - indicator of mode of evaluation
2802 * errcnt - sum of number of func evaluation errors thus far
2803 * newpt - new point indicator
2804 * nj - number of nonlinear nonzero jacobian elements
2805 * n - number of variables
2806 * usrmem - user memory
2807 *
2808 * For the optimization problem at a boundary, this subroutine will
2809 * be called only of the objective function, constraint number m.
2810 *
2811 */
2812 static int COI_CALL slv9_conopt_fdeval(
2813 double *x, double *g, double *jac
2814 , int *rowno, int *jcnm, int *mode, int *ignerr
2815 , int *errcnt, int *newpt, int *n, int *nj
2816 , double *usrmem
2817 ){
2818 slv9_system_t sys;
2819 int32 num_vars, v;
2820 real64 obj, deriv;
2821
2822 UNUSED_PARAMETER(jcnm);
2823 UNUSED_PARAMETER(errcnt);
2824 UNUSED_PARAMETER(newpt);
2825 UNUSED_PARAMETER(n);
2826 UNUSED_PARAMETER(nj);
2827
2828 sys = (slv9_system_t)usrmem;
2829 num_vars = sys->con.n - sys->subregions;
2830
2831 if (*mode == 1 || *mode == 3) {
2832 if (*rowno == sys->con.m - 1){
2833 obj = 0.0;
2834 for (v=0; v<num_vars; v++) {
2835 obj = obj + (x[v] * x[v]);
2836 }
2837 *g = obj;
2838 } else {
2839 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Wrong number of constraints");
2840 return 1;
2841 }
2842 }
2843
2844 /*
2845 * The use of the parameter RTMAXJ is really a Hack to keep CONOPT
2846 * from complaining about large derivatives.
2847 */
2848
2849 if (*mode == 2 || *mode == 3) {
2850 if (*rowno == sys->con.m - 1){
2851 for (v=0; v<num_vars; v++) {
2852 deriv = 2.0 * x[v];
2853 if (deriv > RTMAXJ ) {
2854 deriv = 0.5*RTMAXJ;
2855 } else {
2856 if (deriv < -RTMAXJ ) {
2857 deriv = -0.5*RTMAXJ;
2858 }
2859 }
2860 jac[v] = deriv;
2861 }
2862 } else {
2863 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Wrong number of constraints");
2864 return 1;
2865 }
2866 }
2867
2868 return 0;
2869 }
2870
2871
2872 /*
2873 * COISTA Pass the solution from CONOPT to the modeler. It returns
2874 * completion status
2875 * COISTA(modsta, solsts, iter, objval, usrmem)
2876 *
2877 * modsta - model status
2878 * solsta - solver status
2879 * iter - number of iterations
2880 * objval - objective value
2881 * usrmem - user memory
2882 */
2883 static int COI_CALL slv9_conopt_status(int *modsta, int *solsta, int *iter
2884 , double *objval, double *usrmem
2885 ){
2886 slv9_system_t sys;
2887
2888 sys = (slv9_system_t)usrmem;
2889
2890 sys->con.modsta = *modsta;
2891 sys->con.solsta = *solsta;
2892 sys->con.iter = *iter;
2893 sys->con.obj = *objval;
2894
2895 return 0;
2896 }
2897
2898 /**
2899 CONOPT error message reporting
2900 */
2901 int COI_CALL slv9_conopt_errmsg( int* ROWNO, int* COLNO, int* POSNO, int* MSGLEN
2902 , double* USRMEM, char* MSG, int LENMSG
2903 ){
2904 slv9_system_t sys;
2905 char *varname=NULL;
2906 struct var_variable **vp;
2907
2908 sys = (slv9_system_t)USRMEM;
2909
2910
2911 if(*COLNO!=-1){
2912 vp=sys->mvlist;
2913 vp = vp + *COLNO;
2914 assert(*vp!=NULL);
2915 varname= var_make_name(SERVER,*vp);
2916 }
2917
2918 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERR);
2919 if(*ROWNO == -1){
2920 FPRINTF(ASCERR,"Variable %d (Maybe it's '%s'): ",*COLNO,varname);
2921 ASC_FREE(varname);
2922 }else if(*COLNO == -1 ){
2923 FPRINTF(ASCERR,"Relation %d: ",*ROWNO);
2924 }else{
2925 FPRINTF(ASCERR,"Variable %d (Maybe it's '%s') appearing in relation %d: ",*COLNO,varname,*ROWNO);
2926 ASC_FREE(varname);
2927 }
2928 FPRINTF(ASCERR,"%*s", *MSGLEN, MSG);
2929 error_reporter_end_flush();
2930 return 0;
2931 }
2932
2933
2934 /*
2935 * COIRS Pass the solution from CONOPT to the modeler. It returns
2936 * solution values
2937 * COIRS(val, xmar, xbas, xsta, yval, ymar, ybas, ysta, n, m, usrmem)
2938 *
2939 * xval - the solution values of the variables
2940 * xmar - corresponding marginal values
2941 * xbas - basis indicator for column (at bound, basic, nonbasic)
2942 * xsta - status of column (normal, nonoptimal, infeasible,unbounded)
2943 * yval - values of the left hand side in all the rows
2944 * ymar - corresponding marginal values
2945 * ybas - basis indicator for row
2946 * ysta - status of row
2947 * n - number of variables
2948 * m - number of constraints
2949 * usrmem - user memory
2950 */
2951 static int COI_CALL slv9_conopt_solution(double *xval, double *xmar, int *xbas, int *xsta,
2952 double *yval, double *ymar, int *ybas, int * ysta,
2953 int *n, int *m, double *usrmem
2954 ){
2955 slv9_system_t sys;
2956 struct opt_vector *opt_var_values;
2957 int32 c;
2958 real64 value;
2959
2960 UNUSED_PARAMETER(xmar);UNUSED_PARAMETER(xbas);UNUSED_PARAMETER(xsta);
2961 UNUSED_PARAMETER(yval);UNUSED_PARAMETER(ymar);UNUSED_PARAMETER(ybas);
2962 UNUSED_PARAMETER(ysta);UNUSED_PARAMETER(m);
2963
2964 sys = (slv9_system_t)usrmem;
2965 opt_var_values = sys->opt_var_values;
2966
2967 for (c = 0; c < (*n); c++) {
2968 value = xval[c];
2969 opt_var_values->element[c] = value;
2970 }
2971
2972 return 0;
2973 }
2974
2975 #if 0
2976 /*
2977 * COIUSZ communicates and update of an existing model to CONOPT
2978 * COIUSZ(nintg, ipsz, nreal, rpsz, usrmem)
2979 *
2980 * nintg - number of positions in ipsz
2981 * ipsz - array describing problem size and options
2982 * nreal - number of positions in rpsz
2983 * rpsz - array of reals describing problem size and options
2984 * usrmem- user memory
2985 */
2986 static void slv9_coiusz(int32 *nintg, int32 *ipsz, int32 *nreal, real64 *rpsz,
2987 real64 *usrmem)
2988 {
2989 /* non defined for this solver */
2990
2991 /*
2992 * stop gcc whining about unused parameter
2993 */
2994 (void)nintg; (void)ipsz; (void)nreal; (void)rpsz;
2995 (void)usrmem;
2996
2997 return;
2998 }
2999 #endif
3000
3001 /*
3002 * COIOPT communicates non-default option values to CONOPT
3003 * COIOPT(name, rval, ival, lval, usrmem)
3004 * name - the name of a CONOPT CR-cell defined by the modeler
3005 * rval - the value to be assigned to name if the cells contains a real
3006 * ival - the value to be assigned to name if the cells contains an int
3007 * lval - the value to be assigned to name if the cells contains a log value
3008 * usrmem - user memory
3009 */
3010 static int COI_CALL slv9_conopt_option(
3011 int *NCALL, double *rval, int *ival, int *logical
3012 , double *usrmem, char *name, int lenname
3013 ){
3014 slv9_system_t sys;
3015 sys = (slv9_system_t)usrmem;
3016
3017 UNUSED_PARAMETER(logical);
3018
3019 name = memset(name,' ',8);
3020 while (sys->con.opt_count < slv9_PA_SIZE) {
3021 if (strlen(sys->p.parms[sys->con.opt_count].interface_label) == 6) {
3022 if (strncmp(sys->p.parms[sys->con.opt_count].interface_label,
3023 "R",1) == 0) {
3024 name = strncpy(name, sys->p.parms[sys->con.opt_count]. /* . break */
3025 interface_label,6);
3026 *rval = sys->p.parms[sys->con.opt_count].info.r.value;
3027 sys->con.opt_count++;
3028 return 0;
3029 } else if (strncmp(sys->p.parms[sys->con.opt_count]. /* . break */
3030 interface_label,"L",1) == 0) {
3031 name = strncpy(name,sys->p.parms[sys->con.opt_count]. /* . break */
3032 interface_label,6);
3033 *ival = sys->p.parms[sys->con.opt_count].info.i.value;
3034 sys->con.opt_count++;
3035 return 0;
3036 }
3037 }
3038 sys->con.opt_count++;
3039 }
3040
3041 /* sending blank to quit iterative calling */
3042 name = memset(name,' ',8);
3043 return 0;
3044 }
3045
3046 #if 0 /* see slv_conopt_iterate */
3047 /*
3048 * COIPSZ communicates the model size and structure to CONOPT
3049 * COIPSZ(nintgr, ipsz, nreal, rpsz, usrmem)
3050 *
3051 * ningtr - number of positions in ipsz
3052 * ipsz - array describing problem size and options
3053 * nreal - number of positions in rpsz
3054 * rpsz - array of reals describing problem size and options
3055 * usrmem - user memory
3056 */
3057 static void slv9_coipsz(int32 *nintg, int32 *ipsz, int32 *nreal, real64 *rpsz,
3058 real64 *usrmem)
3059 {
3060 slv9_system_t sys;
3061
3062 /*
3063 * stop gcc whining about unused parameter
3064 */
3065 (void)nintg; (void)nreal;
3066
3067 sys = (slv9_system_t)usrmem;
3068 /*
3069 * Integer array
3070 */
3071 ipsz[F2C(1)] = sys->con.n; /* variables */
3072 ipsz[F2C(2)] = sys->con.m; /* constraints including objective */
3073 ipsz[F2C(3)] = sys->con.nz; /* non zeros in Jacobian */
3074 ipsz[F2C(4)] = sys->con.nz - (sys->con.m - 2); /* linear nonzeros */
3075 ipsz[F2C(5)] = sys->con.m - 2; /* nonlinear nonzeros */
3076 ipsz[F2C(6)] = -1; /* direction of optimization min */
3077 ipsz[F2C(7)] = sys->con.m; /* objective will be last row */
3078 ipsz[F2C(8)] = OPT_ITER_LIMIT; /* iteration limit */
3079 ipsz[F2C(9)] = DOMLIM; /* max number of error in func evals */
3080 ipsz[F2C(10)] = 0; /* output to file */
3081 ipsz[F2C(11)] = 1; /* progress info to screen */
3082 ipsz[F2C(12)] = 1; /* correct value of func in coirms */
3083 ipsz[F2C(13)] = 0; /* not correct value of jacs in coirms */
3084 ipsz[F2C(14)] = 0; /* status not known by modeler */
3085 ipsz[F2C(15)] = 0; /* function value include only NL terms */
3086 ipsz[F2C(16)] = 1; /* Objective is a constraint */
3087 ipsz[F2C(17)] = 0; /* sorted order for jacobian */
3088 ipsz[F2C(18)] = 0; /* append the log file after restarts */
3089 ipsz[F2C(19)] = 0; /* one subroutine call to coirms */
3090 ipsz[F2C(20)] = 0; /* eval subroutine is coifde */
3091 ipsz[F2C(21)] = 0; /* no debugging of derivatives */
3092 ipsz[F2C(22)] = 0; /* coifde not called for linear eqns */
3093 /*
3094 * skipping remainder of ipsz which are fortran io parameters
3095 */
3096
3097 /*
3098 * Real array
3099 */
3100 rpsz[F2C(1)] = ASC_INFINITY; /* infinity */
3101 rpsz[F2C(2)] = -ASC_INFINITY; /* -infinity */
3102 rpsz[F2C(3)] = UNDEFINED; /* undefined */
3103 rpsz[F2C(6)] = 0; /* work space allocated by conopt */
3104 rpsz[F2C(7)] = TIME_LIMIT; /* resource limit (time) */
3105 rpsz[F2C(8)] = 1; /* initial value for vars if none given */
3106
3107
3108 }
3109 #endif
3110
3111
3112 /**
3113 Perform CONOPT solution. For the details of what this does in the larger
3114 context of CMSlv, read (???)
3115
3116 @TODO document this.
3117
3118 @see conopt.h
3119 */
3120 static void slv_conopt_iterate(slv9_system_t sys){
3121
3122 if(sys->con.cntvect == NULL){
3123 sys->con.cntvect = ASC_NEW_ARRAY(int,COIDEF_Size());
3124 }
3125
3126 COIDEF_Ini(sys->con.cntvect);
3127
3128 /*
3129 We pass pointer to sys as usrmem data.
3130 Cast back to slv9_system_t to access the information required
3131 */
3132 COIDEF_UsrMem(sys->con.cntvect,(double *)sys);
3133
3134 COIDEF_NumVar(sys->con.cntvect, &(sys->con.n));
3135 COIDEF_NumCon(sys->con.cntvect, &(sys->con.m)); /* include the obj fn */
3136 COIDEF_NumNZ(sys->con.cntvect, &(sys->con.nz));
3137 COIDEF_NumNlNz(sys->con.cntvect, &(sys->con.nlnz));
3138 COIDEF_OptDir(sys->con.cntvect, &(sys->con.optdir));
3139
3140 COIDEF_ObjCon(sys->con.cntvect, &(sys->con.objcon)); /* objective will be last row */
3141 COIDEF_Base(sys->con.cntvect, &(sys->con.base));
3142 COIDEF_ErrLim(sys->con.cntvect, &(DOMLIM));
3143 COIDEF_ItLim(sys->con.cntvect, &(OPT_ITER_LIMIT));
3144
3145 COIDEF_ReadMatrix(sys->con.cntvect, &slv9_conopt_readmatrix);
3146 COIDEF_FDEval(sys->con.cntvect, &slv9_conopt_fdeval);
3147 COIDEF_Option(sys->con.cntvect, &slv9_conopt_option);
3148 COIDEF_Solution(sys->con.cntvect, &slv9_conopt_solution);
3149 COIDEF_Status(sys->con.cntvect, &slv9_conopt_status);
3150 COIDEF_Message(sys->con.cntvect, &asc_conopt_message);
3151 COIDEF_ErrMsg(sys->con.cntvect, &slv9_conopt_errmsg);
3152 COIDEF_Progress(sys->con.cntvect, &asc_conopt_progress);
3153
3154 /** @TODO implement the following options as well... */
3155 #if 0
3156 ipsz[F2C(10)] = 0; /* output to file */
3157 ipsz[F2C(11)] = 1; /* progress info to screen */
3158 ipsz[F2C(12)] = 1; /* correct value of func in coirms */
3159 ipsz[F2C(13)] = 0; /* not correct value of jacs in coirms */
3160 ipsz[F2C(14)] = 0; /* status not known by modeler */
3161 ipsz[F2C(15)] = 0; /* function value include only NL terms */
3162 ipsz[F2C(16)] = 1; /* Objective is a constraint */
3163 ipsz[F2C(17)] = 0; /* sorted order for jacobian */
3164 ipsz[F2C(18)] = 0; /* append the log file after restarts */
3165 ipsz[F2C(19)] = 0; /* one subroutine call to coirms */
3166 ipsz[F2C(20)] = 0; /* eval subroutine is coifde */
3167 ipsz[F2C(21)] = 0; /* no debugging of derivatives */
3168 ipsz[F2C(22)] = 0; /* coifde not called for linear eqns */
3169 /*
3170 * skipping remainder of ipsz which are fortran io parameters
3171 */
3172
3173 /*
3174 * Real array
3175 */
3176 rpsz[F2C(1)] = ASC_INFINITY; /* infinity */
3177 rpsz[F2C(2)] = -ASC_INFINITY; /* -infinity */
3178 rpsz[F2C(3)] = UNDEFINED; /* undefined */
3179 rpsz[F2C(6)] = 0; /* work space allocated by conopt */
3180 rpsz[F2C(7)] = TIME_LIMIT; /* resource limit (time) */
3181 rpsz[F2C(8)] = 1; /* initial value for vars if none given */
3182 #endif
3183
3184 /*
3185 * reset count on coiopt calls
3186 */
3187 sys->con.opt_count = 0;
3188
3189 /*
3190 * do not keep model in memory after solution
3191 */
3192 sys->con.kept = 0;
3193
3194 COI_Solve(sys->con.cntvect);
3195 /* conopt_start(&(sys->con.kept), usrmem, &(sys->con.lwork),
3196 sys->con.work, &(sys->con.maxusd), &(sys->con.curusd)); */
3197
3198 /*
3199 * We assume that we get convergence in optimization problem at
3200 * boundary
3201 */
3202 sys->con.optimized = 1;
3203 }
3204
3205 #endif /* ASC_WITH_CONOPT */
3206
3207 /*-------------------end of conopt callbacks----------------------------------*/
3208
3209
3210 /*
3211 * Creates an array of columns (containing an array of real elements
3212 * each) to storage the linear coefficient matrix of the optimization problem.
3213 * It also creates the arrays of reals required to storage the values
3214 * of the gradients of a subregion, which change depending on whether the
3215 * problem is a simulation or an optimization.
3216 */
3217 static void create_opt_matrix_and_vectors(int32 num_opt_eqns,
3218 int32 n_subregions,
3219 struct opt_matrix *coeff_matrix,
3220 struct opt_vector *opt_var_values,
3221 struct opt_vector *invariant,
3222 struct opt_vector *variant,
3223 struct opt_vector *gradient,
3224 struct opt_matrix *multipliers)
3225 {
3226 int32 c;
3227 int32 num_vars;
3228
3229 num_vars = num_opt_eqns - 1 + n_subregions;
3230
3231 coeff_matrix->cols = ASC_NEW_ARRAY(struct opt_vector *,n_subregions);
3232
3233 if (g_optimizing) {
3234 multipliers->cols = ASC_NEW_ARRAY(struct opt_vector *,n_subregions);
3235 }
3236
3237 for (c=0; c<n_subregions; c++) {
3238 coeff_matrix->cols[c].element = ASC_NEW_ARRAY(real64,num_opt_eqns);
3239 }
3240 opt_var_values->element = ASC_NEW_ARRAY(real64,num_vars);
3241
3242 if (g_optimizing) {
3243 gradient->element = ASC_NEW_ARRAY(real64,num_opt_eqns);
3244 } else {
3245 invariant->element = ASC_NEW_ARRAY(real64,num_opt_eqns);
3246 variant->element = ASC_NEW_ARRAY(real64,num_opt_eqns);
3247 }
3248 }
3249
3250
3251 /*
3252 * destroy the arrays created to storage the gradients for the optimization
3253 * problem
3254 */
3255 static void destroy_opt_matrix_and_vectors(int32 n_subregions,
3256 struct opt_matrix *coeff_matrix,
3257 struct opt_vector *opt_var_values,
3258 struct opt_vector *invariant,
3259 struct opt_vector *variant,
3260 struct opt_vector *gradient,
3261 struct opt_matrix *multipliers)
3262 {
3263 int32 c;
3264 for (c=0; c<n_subregions; c++) {
3265 destroy_array(coeff_matrix->cols[c].element);
3266 if (g_optimizing) {
3267 destroy_array(multipliers->cols[c].element);
3268 }
3269 }
3270 destroy_array(coeff_matrix->cols);
3271 destroy_array(opt_var_values->element);
3272 if (g_optimizing) {
3273 destroy_array(multipliers->cols);
3274 destroy_array(gradient->element);
3275 } else {
3276 destroy_array(invariant->element);
3277 destroy_array(variant->element);
3278 }
3279 }
3280
3281
3282
3283 /*
3284 * Set Factorization Options
3285 */
3286 static void set_factor_options (linsolqr_system_t lsys)
3287 {
3288 linsolqr_prep(lsys,linsolqr_fmethod_to_fclass(ranki_ba2));
3289 linsolqr_set_pivot_zero(lsys, 1e-12);
3290 linsolqr_set_drop_tolerance(lsys,1e-16);
3291 linsolqr_set_pivot_tolerance(lsys, 0.1);
3292 linsolqr_set_condition_tolerance(lsys, 0.1);
3293 }
3294
3295
3296 /*
3297 * Calculating the Lagrange Multipliers for each subregion
3298 *
3299 * We are assuming here that the matrix is structurally nonsingular
3300 * and than the rank is equal to the number of rows in the matrix.
3301 * Much more efficient checking must be done.
3302 */
3303 static void get_multipliers(SlvClientToken asys,
3304 int32 subregion,
3305 int32 nrel,
3306 real64 *grad_obj,
3307 struct opt_matrix *multipliers)
3308 {
3309 slv9_system_t sys;
3310 linsolqr_system_t lsys;
3311 mtx_region_t *newblocks, *oneblock;
3312 int32 rank;
3313 int32 c, cr, len, row;
3314 real64 *weights;
3315 real64 summ;
3316
3317 sys = SLV9(asys);
3318 check_system(sys);
3319
3320 mtx_output_assign(sys->lin_mtx,nrel,nrel);
3321 rank = mtx_symbolic_rank(sys->lin_mtx);
3322 mtx_partition(sys->lin_mtx);
3323 len = mtx_number_of_blocks(sys->lin_mtx);
3324 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
3325 if (newblocks == NULL) {
3326 mtx_destroy(sys->lin_mtx);
3327 return;
3328 }
3329 for (c = 0 ; c < len; c++) {
3330 mtx_block(sys->lin_mtx,c,&(newblocks[c]));
3331 }
3332 for (c = 0 ; c < len; c++) {
3333 mtx_reorder(sys->lin_mtx,&(newblocks[c]),mtx_SPK1);
3334 }
3335
3336 /* unifying block */
3337 oneblock = (mtx_region_t *)ascmalloc(sizeof(mtx_region_t));
3338 oneblock->row.low = oneblock->col.low = 0;
3339 oneblock->row.high = nrel-1;
3340 oneblock->col.high = nrel-1;
3341
3342 /*
3343 * Scaling of the linear system
3344 */
3345
3346 /*
3347 *Calculating weights
3348 */
3349 weights = ASC_NEW_ARRAY(real64,nrel);
3350 for (row=0; row<nrel; row++) {
3351 summ = mtx_sum_sqrs_in_row(sys->lin_mtx,row,&(oneblock->col));
3352 if (summ <= 0.0) {
3353 weights[row] = 1.0;
3354 } else {
3355 weights[row] = 1.0 / sqrt(summ);
3356 }
3357 #if DEBUG
3358 FPRINTF(ASCERR," weight of row %d = %f \n",row,summ);
3359 #endif /* DEBUG */
3360 }
3361
3362 /*
3363 * Dividing rows by weights
3364 */
3365 for (row=0; row<nrel; row++) {
3366 mtx_mult_row(sys->lin_mtx,row,weights[row],&(oneblock->col));
3367 }
3368
3369 /*
3370 * dividing rhs
3371 */
3372 for (row=0; row<nrel; row++) {
3373 grad_obj[mtx_row_to_org(sys->lin_mtx,row)] =
3374 grad_obj[mtx_row_to_org(sys->lin_mtx,row)] * weights[row];
3375 }
3376
3377 /*
3378 * End of scaling
3379 */
3380
3381 lsys = linsolqr_create();
3382 linsolqr_set_matrix(lsys,sys->lin_mtx);
3383
3384 for (cr=0; cr<nrel; cr++) {
3385 multipliers->cols[subregion].element[cr] = 0.0;
3386 }
3387
3388 set_factor_options(lsys);
3389 /* rhs for multipliers */
3390 linsolqr_add_rhs(lsys,grad_obj,FALSE);
3391 linsolqr_set_region(lsys,*oneblock);
3392 linsolqr_factor(lsys,ranki_ba2);
3393 linsolqr_solve(lsys,grad_obj);
3394 for (cr=0; cr<nrel; cr++) {
3395 multipliers->cols[subregion].element[cr] = linsolqr_var_value
3396 (lsys,grad_obj,cr);
3397 #if SHOW_LAGRANGE_DETAILS
3398 FPRINTF(ASCERR, " Row = %d \n",cr);
3399 FPRINTF(ASCERR,
3400 "Multiplier = %f \n",multipliers->cols[subregion].element[cr]);
3401 #endif /* SHOW_LAGRANGE_DETAILS */
3402 }
3403 linsolqr_set_matrix(lsys,NULL);
3404 mtx_destroy(sys->lin_mtx);
3405 linsolqr_destroy(lsys);
3406 destroy_array(newblocks);
3407 destroy_array(weights);
3408 ascfree(oneblock);
3409 }
3410
3411
3412 /*
3413 * Calculate the invariant part of the gradients of the subregions
3414 */
3415 static void get_gradient_in_subregion(slv_system_t server,
3416 SlvClientToken asys,
3417 int32 subregion,
3418 int32 num_opt_eqns,
3419 struct opt_vector *gradient,
3420 struct opt_matrix *multipliers)
3421 {
3422 slv9_system_t sys;
3423 struct rel_relation **rlist;
3424 struct var_variable **vlist;
3425 struct rel_relation *rel;
3426 struct var_variable *var;
3427 var_filter_t vfilter;
3428 rel_filter_t rfilter;
3429 mtx_coord_t coord;
3430 real64 *tmp_value;
3431 real64 *derivatives, resid;
3432 real64 *grad_obj, *func_val;
3433 real64 *f_red_grad;
3434 struct opt_matrix rel_red_grad;
3435 int32 *variables_master, *variables_solver, count;
3436 int32 nvar, nrel, ntotvar, ntotrel;
3437 int32 countrel,countvar,cr,cv,len,vind;
3438 int32 nvnb, countnbv;
3439 FILE *lif;
3440
3441 sys = SLV9(asys);
3442 check_system(sys);
3443 lif = LIF(sys);
3444
3445 rlist = slv_get_master_rel_list(server);
3446 vlist = slv_get_master_var_list(server);
3447 ntotvar = slv_get_num_master_vars(server);
3448 ntotrel = slv_get_num_master_rels(server);
3449 tmp_value = ASC_NEW_ARRAY(real64,ntotvar);
3450
3451 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_NONBASIC
3452 | VAR_SVAR | VAR_FIXED);
3453 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
3454 nvar = slv_count_master_vars(server,&vfilter);
3455
3456 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE );
3457 rfilter.matchvalue =(REL_INCLUDED | REL_EQUALITY | REL_ACTIVE );
3458 nrel = slv_count_master_rels(server,&rfilter);
3459
3460 if (nrel != nvar) {
3461 FPRINTF(ASCERR," nrel = %d\n",nrel);
3462 FPRINTF(ASCERR," nvar = %d\n",nvar);
3463 FPRINTF(ASCERR,
3464 "PANIC: number relations does not match number of variables\n");
3465 }
3466
3467 /*
3468 * residual of the relations in the subregion
3469 */
3470 func_val = ASC_NEW_ARRAY(real64,nrel);
3471
3472 /*
3473 * Lagrange Multipliers of the subregion
3474 */
3475 multipliers->cols[subregion].element =
3476 (real64 *)(ascmalloc(nrel*sizeof(real64)));
3477 /*
3478 * Gradients of the objective function
3479 */
3480 grad_obj = ASC_NEW_ARRAY(real64,nvar);
3481
3482 /*
3483 * Matrix for solving linear system
3484 */
3485 sys->lin_mtx = mtx_create();
3486 mtx_set_order(sys->lin_mtx,nrel);
3487
3488 /*
3489 * Counting nonbasic variables
3490 */
3491 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_NONBASIC
3492 | VAR_SVAR | VAR_FIXED);
3493 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_NONBASIC
3494 | VAR_SVAR);
3495 nvnb = slv_count_master_vars(server,&vfilter);
3496
3497 /*
3498 * Information for reduced gradient
3499 */
3500 f_red_grad = ASC_NEW_ARRAY(real64,nvnb);
3501 rel_red_grad.cols = (struct opt_vector *)
3502 (ascmalloc(nvnb*sizeof(struct opt_vector)));
3503 for (cv=0; cv<nvnb; cv++) {
3504 rel_red_grad.cols[cv].element =
3505 (real64 *)(ascmalloc(nrel*sizeof(real64)));
3506 }
3507
3508 /*
3509 * Setting every sindex to -1
3510 */
3511 for (cv=0; cv<ntotvar; cv++) {
3512 var_set_sindex(vlist[cv],-1);
3513 }
3514 for (cr=0; cr<ntotrel; cr++) {
3515 rel_set_sindex(rlist[cr],-1);
3516 }
3517
3518 /*
3519 * Initializing values
3520 */
3521 for (cv=0; cv<ntotvar;cv++) {
3522 tmp_value[cv] = 0.0;
3523 }
3524 for (cv=0; cv<num_opt_eqns;cv++) {
3525 gradient->element[cv] = 0.0;
3526 }
3527
3528 for (cr=0; cr<nrel; cr++) {
3529 func_val[cr] = 0.0;
3530 multipliers->cols[subregion].element[cr] = 0.0;
3531 }
3532
3533 for (cv=0; cv<nvnb; cv++) {
3534 f_red_grad[cv] = 0.0;
3535 for (cr=0; cr<nrel; cr++) {
3536 rel_red_grad.cols[cv].element[cr] = 0.0;
3537 }
3538 }
3539
3540 for (cv=0; cv<nvar; cv++) {
3541 grad_obj[cv] = 0.0;
3542 }
3543
3544 /*
3545 * Calculate Values
3546 */
3547 vfilter.matchbits = (VAR_ACTIVE_AT_BND | VAR_INCIDENT
3548 | VAR_SVAR | VAR_FIXED);
3549 vfilter.matchvalue = (VAR_ACTIVE_AT_BND | VAR_INCIDENT | VAR_SVAR);
3550
3551 /*
3552 * List of relations
3553 */
3554 countrel = 0;
3555 countvar = 0;
3556 countnbv = 0;
3557 for (cr=0; cr<ntotrel; cr++) {
3558 rel = rlist[cr];
3559 if (rel_apply_filter(rel,&rfilter)) {
3560 rel_set_sindex(rel,countrel);
3561 coord.col = rel_sindex(rel);
3562 len = rel_n_incidences(rel);
3563 variables_master = ASC_NEW_ARRAY(int32,len);
3564 variables_solver = ASC_NEW_ARRAY(int32,len);
3565 derivatives = ASC_NEW_ARRAY(real64,len);
3566 relman_diff_grad(rel,&vfilter,derivatives,variables_master,
3567 variables_solver,&count,&resid,1);
3568 func_val[countrel] = resid;
3569 #if SHOW_LAGRANGE_DETAILS
3570 FPRINTF(ASCERR,"Equation = %d \n",coord.col);
3571 FPRINTF(ASCERR,"Residual = %f \n",resid);
3572 #endif /* SHOW_LAGRANGE_DETAILS */
3573 for (cv=0; cv<count;cv++) {
3574 var = vlist[variables_master[cv]];
3575 if (!var_nonbasic(var)) {
3576 tmp_value[variables_master[cv]] = tmp_value[variables_master[cv]] +
3577 derivatives[cv] * RHO * resid;
3578 if(var_active(var)) {
3579 coord.row = var_sindex(var);
3580 if (coord.row == -1) {
3581 var_set_sindex(var,countvar);
3582 coord.row = countvar;
3583 countvar++;
3584 }
3585 assert(coord.col >= 0 && coord.col < mtx_order(sys->lin_mtx));
3586 #if SHOW_LAGRANGE_DETAILS
3587 FPRINTF(ASCERR,"Coordinate row = %d \n",coord.row);
3588 FPRINTF(ASCERR,"Coordinate col = %d \n",coord.col);
3589 FPRINTF(ASCERR,"Derivative = %f \n",derivatives[cv]);
3590 #endif /* SHOW_LAGRANGE_DETAILS */
3591 mtx_fill_org_value(sys->lin_mtx,&coord,derivatives[cv]);
3592 }
3593 } else {
3594 if (var_sindex(var)== -1) {
3595 var_set_sindex(var,countnbv);
3596 countnbv++;
3597 }
3598 rel_red_grad.cols[var_sindex(var)].element[countrel] =
3599 derivatives[cv];
3600 #if SHOW_LAGRANGE_DETAILS
3601 FPRINTF(lif,"Nonbasic Variable ");
3602 print_var_name(lif,sys,var); PUTC('\n',lif);
3603 FPRINTF(ASCERR,"Derivative = %f \n",derivatives[cv]);
3604 #endif /* SHOW_LAGRANGE_DETAILS */
3605 }
3606 }
3607 destroy_array(variables_master);
3608 destroy_array(variables_solver);
3609 destroy_array(derivatives);
3610 countrel++;
3611 }
3612 }