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