/[ascend]/trunk/solvers/cmslv/cmslv.c
ViewVC logotype

Contents of /trunk/solvers/cmslv/cmslv.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2036 - (show annotations) (download) (as text)
Mon May 18 15:03:16 2009 UTC (13 years, 6 months ago) by ballan
File MIME type: text/x-csrc
File size: 163313 byte(s)
issues resolved:
295
390
301
cmslv.c: unused var cleanup.
lsode/SConscript: fortran flags bugs-- may not work with 0.9x scons.
works with 1.2+. when adding -w, or any special flags, be sure to add
them and not replace the original flag.
system/var.c: 64bit clarity.
system/discrete.c: 64bit clarity.
system/analyze*: g_reuse declared in wrong place. 64bit clarity
system/diffvars: missing prototype function, 64bit clarity.
compiler/numlist.*: changed from int to glint.
compiler/simlist.c: missing includes needed for 64bit clarity.
compiler/instance_io.c: missing includes needed for 64bit clarity.
compiler/initialize.[ch]: const clarifications.
compiler/packages.c: const clarifications.
compiler/module.c: const clarifications.
compiler/statio.c: unused var cleanup.
compiler/procframe; const clarification. memory deallocation bugs.
compiler/notequery.c: repaired multiple casting and 64bit issues.
compiler/importhandler.c: const and free issues fixed.
compiler/type_desc.c: ridiculous if constructs clarified.
compiler/createinst.c: casting stupidity repaired.
linear/ranki2.c: missing includes needed for 64bitness.
solver/solver.c: const issues clarified.
utilities/ascConfig.h: added GLint typedefs for dealing with gllist
64bit portability.
utilities/ascPanic.c: removed extraneous const.
general/ospath.c: safer,quieter handling for string pointer difference.
integrator/integrator.c: const issues clarified.
packages/sensitivity.c: missing includes needed fo 64bit sanity.
tcltk/interface/Integrators.c: 64 bitness.
tcltk/interface/SimsProc.c: const errors.
tcltk/interface/Driver.c: fixed env var handling wrt ascend-config (295)
models/test/z*a4c: fixed meters -> m conversion; someone never ran the
test suite after teasing the default units to ambiguous abbreviations.
SConstruct: added sizeof checks; output might be better put in a ascend
system-wide header.


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