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

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

Parent Directory Parent Directory | Revision Log Revision Log


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