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

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

Parent Directory Parent Directory | Revision Log Revision Log


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