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

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

Parent Directory Parent Directory | Revision Log Revision Log


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