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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 58663 byte(s)
Setting up web subdirectory in repository
1 /*
2 * MPS: Ascend MPS file generator
3 * by Craig Schmidt
4 * Created: 2/11/95
5 * Version: $Revision: 1.29 $
6 * Version control file: $RCSfile: slv6.c,v $
7 * Date last modified: $Date: 2000/01/25 02:27:38 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1990 Karl Michael Westerberg
13 * Copyright (C) 1993 Joseph Zaher
14 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15 * Copyright (C) 1995 Craig Schmidt
16 *
17 * The SLV solver is free software; you can redistribute
18 * it and/or modify it under the terms of the GNU General Public License as
19 * published by the Free Software Foundation; either version 2 of the
20 * License, or (at your option) any later version.
21 *
22 * The SLV solver is distributed in hope that it will be
23 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25 * General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with the program; if not, write to the Free Software Foundation,
29 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
30 * COPYING. COPYING is found in ../compiler.
31 *
32 */
33 /* known bugs
34 * still uses pl_ functions and assumes the old slv protocol.
35 */
36 #include "utilities/ascConfig.h"
37 #include "utilities/ascMalloc.h"
38 #include "utilities/set.h"
39 #include "utilities/mem.h"
40 #include "general/tm_time.h"
41 #include "general/list.h"
42 #include "general/dstring.h"
43 #include "compiler/module.h"
44 #include "compiler/compiler.h"
45 #include "compiler/library.h"
46 #include "compiler/instance.h"
47 #include "compiler/instance_io.h"
48 #include "solver/mtx.h"
49 #include "solver/linsol.h"
50 #include "solver/linsolqr.h"
51 #include "solver/slv_types.h"
52 #include "solver/var.h"
53 #include "solver/rel.h"
54 #include "solver/discrete.h"
55 #include "solver/conditional.h"
56 #include "solver/logrel.h"
57 #include "solver/bnd.h"
58 #include "solver/calc.h"
59 #include "solver/relman.h"
60 #include "solver/slv_common.h"
61 #include "solver/slv_client.h"
62 #include "solver/slv6.h"
63 #include "solver/mps.h"
64 #include "interface/old_utils.h"
65
66 #if !defined(STATIC_MPS) && !defined(DYNAMIC_MPS)
67 int slv6_register(SlvFunctionsT *f)
68 {
69 (void)f; /* stop gcc whine about unused parameter */
70
71 FPRINTF(stderr,"makeMPS not compiled in this ASCEND IV.\n");
72 return 1;
73 }
74 #else /* either STATIC_MPS or DYNAMIC_MPS is defined */
75 #ifdef DYNAMIC_MPS
76 /* do dynamic loading stuff. yeah, right */
77 #else /* following is used if STATIC_MPS is defined */
78
79 #ifndef KILL
80 #define KILL TRUE
81 #endif
82 #define DEBUG FALSE
83
84 struct slv6_system_structure {
85
86 /**
87 *** Problem definition
88 **/
89 slv_system_t slv; /* slv_system_t back-link */
90 struct rel_relation *obj; /* Objective function: NULL = none */
91 struct var_variable **vlist; /* Variable list (NULL terminated) */
92 struct var_variable **vlist_user; /* User vlist (NULL = determine) */
93 bnd_boundary_t *blist; /* Boundary list (NULL terminated) */
94 bnd_boundary_t *blist_user; /* User blist (NULL = none) */
95 struct rel_relation **rlist; /* Relation list (NULL terminated) */
96 struct rel_relation **rlist_user; /* User rlist (NULL = none) */
97 struct ExtRelCache **erlist; /* External relations cache list */
98 struct ExtRelCache **erlist_user;/* User erlist (NULL = none) */
99
100 /**
101 *** Solver information
102 **/
103 int integrity; /* ? Has the system been created */
104 slv_parameters_t p; /* Parameters */
105 slv_status_t s; /* Status flags */
106 double clock; /* CPU time */
107 int iarray[slv6_IA_SIZE]; /* Integer subparameters */
108 double rarray[slv6_RA_SIZE]; /* Real subparameters */
109 char *carray[slv6_CA_SIZE]; /* Charptr subparameter */
110
111 /**
112 *** Calculated Data
113 ***
114 **/
115 mps_data_t mps; /* the main chunk of data for the problem */
116
117 };
118
119 /* _________________________________________________________________________ */
120
121 /**
122 *** Integrity checks
123 *** ----------------
124 *** check_system(sys)
125 **/
126
127 #define OK ((int)813025392)
128 #define DESTROYED ((int)103289182)
129 static int check_system(slv6_system_t sys)
130 /**
131 *** Checks sys for NULL and for integrity.
132 **/
133 {
134 if( sys == NULL ) {
135 FPRINTF(stderr,"ERROR: (slv6) check_system\n");
136 FPRINTF(stderr," NULL system handle.\n");
137 return 1;
138 }
139
140 switch( sys->integrity ) {
141 case OK:
142 return 0;
143 case DESTROYED:
144 FPRINTF(stderr,"ERROR: (slv6) check_system\n");
145 FPRINTF(stderr," System was recently destroyed.\n");
146 return 1;
147 default:
148 FPRINTF(stderr,"ERROR: (slv6) check_system\n");
149 FPRINTF(stderr," System reused or never allocated.\n");
150 return 1;
151 }
152 }
153
154 /* _________________________________________________________________________ */
155
156 /**
157 *** Array/vector operations
158 *** ----------------------------
159 *** destroy_array(p)
160 *** create_array(len,type)
161 *** zero_array(arr,len,type)
162 *** nuke_pointers(mps_data_t mps) - free allocated memory in mps
163 **/
164
165 #define destroy_array(p) \
166 if( (p) != NULL ) ascfree((p))
167 #define create_array(len,type) \
168 ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
169 #define create_zero_array(len,type) \
170 ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
171 #define zero_array(arr,nelts,type) \
172 mem_zero_byte_cast((arr),0,(nelts)*sizeof(type))
173 /* Zeros an array of nelts objects, each having given type. */
174
175
176 static void nuke_pointers(mps_data_t mps) { /* free all allocated memory in mps data structure */
177
178 if (mps.Ac_mtx != NULL) { /* delete old matrix if the exist */
179 mtx_destroy(mps.Ac_mtx);
180 mps.Ac_mtx = NULL;
181 }
182
183 if (mps.lbrow != NULL) { /* delete old vector if it exists */
184 destroy_array(mps.lbrow);
185 mps.lbrow = NULL;
186 }
187
188 if (mps.ubrow != NULL) { /* delete old vector if it exists */
189 destroy_array(mps.ubrow);
190 mps.ubrow = NULL;
191 }
192
193 if (mps.bcol != NULL) { /* delete old vector if the exist */
194 destroy_array(mps.bcol);
195 mps.bcol = NULL;
196 }
197
198 if (mps.typerow != NULL) { /* delete old vector if it exists */
199 destroy_array(mps.typerow);
200 mps.typerow = NULL;
201 }
202
203 if (mps.relopcol != NULL) { /* delete old vector if it exists */
204 destroy_array(mps.relopcol);
205 mps.relopcol = NULL;
206 }
207 }
208
209 /* _________________________________________________________________________ */
210
211 /**
212 *** General input/output routines
213 *** -----------------------------
214 *** fp = MIF(sys)
215 *** fp = LIF(sys)
216 **/
217
218 static FILE *get_output_file(FILE *fp)
219 /**
220 *** Returns fp if fp!=NULL, or a file pointer
221 *** open to nul device if fp == NULL.
222 **/
223 {
224 static FILE *nuldev = NULL;
225 static char fname[] = "/dev/null";
226
227 if( fp==NULL ) {
228 if(nuldev==NULL)
229 if( (nuldev=fopen(fname,"w")) == NULL ) {
230 FPRINTF(stderr,"ERROR: (slv6) get_output_file\n");
231 FPRINTF(stderr," Unable to open %s.\n",fname);
232 }
233 fp=nuldev;
234 }
235 return(fp);
236 }
237
238 /* #define MIF(sys) get_output_file( (sys)->p.output.more_important )
239 * #define LIF(sys) get_output_file( (sys)->p.output.less_important )
240 */
241
242 /* _________________________________________________________________________ */
243
244 /**
245 *** Routines for common filters
246 *** -----------------
247 *** free_inc_var_filter - true for non-fixed incident variables
248 *** inc_rel_filter - true for incident relations
249 **/
250
251 extern boolean free_inc_var_filter(struct var_variable *var)
252 /**
253 *** I've been calling this particular var filter a lot ,
254 *** so I decided to make it a subroutine. Returns true if
255 *** var is not fixed and incident in something.
256 **/
257 {
258 var_filter_t vfilter;
259 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
260 vfilter.matchvalue = (VAR_INCIDENT | VAR_ACTIVE);
261
262 /* vfilter.fixed = var_false;*/ /* calc for all non-fixed vars */
263 /* vfilter.incident = var_true; */ /* incident vars only */
264 /* vfilter.in_block = var_ignore; */
265
266 return var_apply_filter(var,&vfilter);
267 }
268
269 static boolean inc_rel_filter(struct rel_relation *rel)
270 /**
271 *** Returns true if rel is an incident relation.
272 **/
273 {
274 rel_filter_t rfilter; /* filter for included rels */
275 rfilter.matchbits = (REL_INCLUDED | REL_ACTIVE);
276 rfilter.matchvalue = (REL_INCLUDED| REL_ACTIVE );
277 /* rfilter.included = rel_true;
278 rfilter.equality = rel_ignore;
279 rfilter.in_block = rel_ignore;
280 rfilter.in_subregion = rel_ignore; */
281
282 return rel_apply_filter(rel,&rfilter);
283 }
284
285
286 /* _________________________________________________________________________ */
287
288 /**
289 *** Routines to calculate the mps problem representation
290 *** --------------------
291 *** var_relaxed - is the variable relaxed or not?
292 *** calc_c - calculate c vector, coefficients of objective
293 (called by calc_matrix)
294 *** calc_bounds - convert var bounds to an array of numbers
295 *** calc_reloplist - convert relation operators <=, >=, = to array of numbers
296 *** calc_svtlist - create array of numbers containing var types
297 *** calc_matrix - compute entire matrix representaion of problem
298 **/
299
300 static boolean calc_c(mtx_matrix_t mtx, /* matrix to store derivs */
301 int32 org_row, /* original number of row to store them */
302 struct rel_relation *obj) /* expression to diffs */
303 /**
304 *** Calculate gradient of the objective function. (or any expression, for that matter)
305 *** On the linear system we should have, this is the c vector of
306 *** our problem max/min {cx: Ax<=b}.
307 *** On nonlinear problems is the linearization of problem at current point
308 **/
309 {
310 var_filter_t vfilter;
311 int32 row;
312
313 if ((mtx == NULL) || (obj == NULL)) { /* got a bad pointer */
314 FPRINTF(stderr,"ERROR: (slv6) calc_c\n");
315 FPRINTF(stderr," Routine was passed a NULL pointer!\n");
316 return FALSE;
317 }
318
319 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
320 vfilter.matchvalue = (VAR_INCIDENT | VAR_ACTIVE);
321 /* vfilter.fixed = var_false;
322 vfilter.incident = var_true;
323 vfilter.in_block = var_ignore; */
324
325 row = mtx_org_to_row(mtx,org_row); /* convert from original numbering to current */
326
327 #ifndef KILL
328 exprman_diffs(obj, &vfilter, row, mtx); /* find the deriviative, sets calc_ok as result */
329 #else
330 FPRINTF(stderr," function calc_c is slv6 is broken.\n");
331 FPRINTF(stderr," It calls exprman_diffs with wrong number of args.\n");
332 exit;
333 #endif
334 return calc_ok; /* did things work out ?, calc_ok in calc.h */
335 }
336
337
338 static real64 *calc_bounds(struct var_variable **vlist, /* variable list to get bounds */
339 int32 vinc, /* number of incident variables */
340 boolean upper) /* do upper, else lower */
341
342 /**
343 ** Stores the upper or lower bounds of all non-fixed, incident vars
344 ** in a new array, which is returned by the routine
345 **/
346 {
347 real64 *tmp,*tmp_array_origin; /* temporary storage for our bounds data */
348
349 if (vlist == NULL) { /* got a bad pointer */
350 FPRINTF(stderr,"ERROR: (slv6) calc_bounds\n");
351 FPRINTF(stderr," Routine was passed a NULL variable list pointer!\n");
352 return FALSE;
353 }
354
355 tmp_array_origin = create_array(vinc,real64);
356 if (tmp_array_origin == NULL) {
357 FPRINTF(stderr,"ERROR: (slv6) calc_bounds\n");
358 FPRINTF(stderr," Memory allocation failed!\n");
359 return FALSE;
360 }
361
362 tmp = tmp_array_origin;
363 for( ; *vlist != NULL ; ++vlist )
364 if( free_inc_var_filter(*vlist) )
365 if (upper) { *tmp=var_upper_bound(*vlist); tmp++; }
366 else { *tmp=var_lower_bound(*vlist); tmp++; }
367
368 return tmp_array_origin;
369 }
370
371
372 static char *calc_reloplist(struct rel_relation **rlist,
373 int32 rused) /* entry for each relation */
374 /**
375 *** This function constructs the a list of relational operators: <=, >=, =
376 *** from the relations list. The values for each relational operator
377 *** corresponds to rel_TOK_less, rel_TOK_equal, and rel_TOK_greater.
378 *** (Defined in rel.h)
379 *** Or rel_TOK_nonincident for relations that aren't incident (see slv6.h)
380 ***
381 *** Note: the rel_less, rel_equal, and rel_greater routines in rel.h don't
382 *** behave as you'd expect, and should _not_ be used. Use rel_type instead.
383 **/
384 {
385 char *reloplist;
386 char *tmp; /* pointer for storage */
387
388 reloplist = create_array(rused,char); /* see macro, over all incident relations */
389 if (reloplist == NULL) { /* memory allocation failed */
390 FPRINTF(stderr,"ERROR: (slv6) calc_reloplist\n");
391 FPRINTF(stderr," Memory allocation failed!\n");
392 return NULL;
393 }
394
395 tmp = reloplist;
396 for ( ;*rlist != NULL; rlist++)
397 {
398 if (inc_rel_filter(*rlist)) /* is an incident var */
399 switch (rel_type(*rlist)) {
400 case e_less:
401 case e_lesseq:
402 *tmp = rel_TOK_less;
403 break;
404
405 case e_equal:
406 *tmp = rel_TOK_equal;
407 break;
408 case e_greater:
409 case e_greatereq:
410 *tmp = rel_TOK_greater;
411 break;
412 default:
413 FPRINTF(stderr,"ERROR: (slv6) calc_reloplist\n");
414 FPRINTF(stderr," Unknown relation type (not greater, less, or equal)\n");
415 return NULL;
416 }
417 else
418 *tmp = rel_TOK_nonincident;
419
420 tmp++; /* increment to new location in array */
421
422 }
423
424 return reloplist;
425 }
426
427
428 static char *calc_svtlist( struct var_variable **vlist, /* input, not modified */
429 int32 vused, /* number of vars (incident or nonincident, free or fixed) */
430 int *solver_var_used, /* number of each type of var are cached */
431 int *solver_relaxed_used,
432 int *solver_int_used,
433 int *solver_binary_used,
434 int *solver_semi_used,
435 int *solver_other_used,
436 int *solver_fixed)
437 /**
438 *** This function constructs the solver var list from the variable list.
439 ***
440 *** WARNING: This routine assumes that a struct var_variable *is an Instance.
441 *** In the future this is going to change, and this routine will break.
442 ***
443 **/
444 {
445 struct TypeDescription *type; /* type of the current var */
446 struct TypeDescription *solver_var_type; /* type of the standard types */
447 struct TypeDescription *solver_int_type;
448 struct TypeDescription *solver_binary_type;
449 struct TypeDescription *solver_semi_type;
450
451 char *svtlist; /* pointer for storage */
452 char *tmp; /* temp pointer to set values */
453
454 /* get the types for variable definitions */
455
456 if( (solver_var_type = FindType(SOLVER_VAR_STR)) == NULL ) {
457 FPRINTF(stderr,"ERROR: (slv6.c) get_solver_var_type\n");
458 FPRINTF(stderr," Type solver_var not defined.\n");
459 FPRINTF(stderr," MPS will not work.\n");
460 return NULL;
461 }
462 if( (solver_int_type = FindType(SOLVER_INT_STR)) == NULL ) {
463 FPRINTF(stderr,"ERROR: (slv6.c) get_solver_var_type\n");
464 FPRINTF(stderr," Type solver_int not defined.\n");
465 FPRINTF(stderr," MPS will not work.\n");
466 return NULL;
467 }
468 if( (solver_binary_type = FindType(SOLVER_BINARY_STR)) == NULL ) {
469 FPRINTF(stderr,"ERROR: (slv6.c) get_solver_var_type\n");
470 FPRINTF(stderr," Type solver_binary not defined.\n");
471 FPRINTF(stderr," MPS will not work.\n");
472 return NULL;
473 }
474 if( (solver_semi_type = FindType(SOLVER_SEMI_STR)) == NULL ) {
475 FPRINTF(stderr,"ERROR: (slv6.c) get_solver_var_type\n");
476 FPRINTF(stderr," Type solver_semi not defined.\n");
477 FPRINTF(stderr," MPS will not work.\n");
478 return NULL;
479 }
480
481 /* allocate memory and initialize stuff */
482 svtlist = create_array(vused,char); /* see macro */
483 if (svtlist == NULL) { /* memory allocation failed */
484 FPRINTF(stderr,"ERROR: (slv6) calc_svtlist\n");
485 FPRINTF(stderr," Memory allocation failed for solver var type list!\n");
486 return NULL;
487 }
488
489 *solver_var_used = 0;
490 *solver_relaxed_used = 0;
491 *solver_int_used = 0;
492 *solver_binary_used = 0;
493 *solver_semi_used = 0;
494 *solver_other_used = 0;
495 *solver_fixed = 0;
496 tmp = svtlist;
497
498 /* loop over all vars */
499
500 for( ; *vlist != NULL ; ++vlist ) {
501 if( free_inc_var_filter(*vlist) )
502 {
503 type = InstanceTypeDesc( (struct Instance *) *vlist);
504
505 if (type == MoreRefined(type,solver_binary_type) )
506 {
507 if (var_relaxed(*vlist))
508 {
509 *tmp = SOLVER_RELAXED;
510 *solver_relaxed_used++;
511 }
512 else
513 {
514 *tmp = SOLVER_BINARY;
515 *solver_binary_used++;
516 }
517 }
518 else
519 {
520 if (type == MoreRefined(type,solver_int_type) )
521 {
522 if (var_relaxed(*vlist))
523 {
524 *tmp = SOLVER_RELAXED;
525 *solver_relaxed_used++;
526 }
527 else
528 {
529 *tmp = SOLVER_INT;
530 *solver_int_used++;
531 }
532 }
533 else
534 {
535 if (type == MoreRefined(type,solver_semi_type) )
536 {
537 if (var_relaxed(*vlist))
538 {
539 *tmp = SOLVER_RELAXED;
540 *solver_relaxed_used++;
541 }
542 else
543 {
544 *tmp = SOLVER_SEMI;
545 *solver_semi_used++;
546 }
547 }
548 else
549 {
550 if (type == MoreRefined(type,solver_var_type) )
551 { /* either solver var or some refinement */
552 *tmp = SOLVER_VAR;
553 *solver_var_used++;
554 }
555 else
556 {
557 FPRINTF(stderr,"ERROR: (slv6) determine_svtlist\n");
558 FPRINTF(stderr," Unknown solver_var type encountered.\n");
559 /* should never get to here */
560 }
561 } /* if semi */
562 } /* if int */
563 } /* if binary */
564 } /* if free inc var */
565 else
566 {
567 *tmp = SOLVER_FIXED;
568 *solver_fixed++;
569 }
570
571 tmp++;
572
573 } /* for */
574
575 return svtlist;
576
577 }
578
579 static mtx_matrix_t calc_matrix(int32 cap,
580 int32 rused,
581 int32 vused,
582 struct rel_relation **rlist,
583 struct rel_relation *obj,
584 int32 crow,
585 slv_status_t *s,
586 int32 *rank,
587 real64 **rhs_orig)
588 /**
589 int32 cap, in: capacity of matrix
590 int32 rused, in: total number of relations used
591 int32 vused, in: total number of variables used
592 struct rel_relation **rlist, in: Relation list (NULL terminated)
593 struct rel_relation *obj, in: objective function
594 int32 crow, in: row to store objective row
595 slv_status_t *s, out: s.block.jactime, and s.calc_ok
596 int32 *rank, out
597 real64 **rhs_orig out: rhs array origin
598 **/
599 /**
600 *** Creates and calculates a matrix representation of the A matrix, c row,
601 *** and the RHS or b column (which is stored in rhs_orig).
602 *** On nonlinear problems is the linearization of problem at current point.
603 ***
604 *** Note: the residual stored in the rhs array is not the real right hand side.
605 *** The residual returned by the diffs call is just the value of
606 ** (lhs expr) - (rhs expr)
607 *** we want the residual excluding the current variables.
608 *** At the moment there isn't a
609 *** clean way to do this. It will be calculated in the real_rhs routine.
610 *** This routine
611 *** will take for each relation:
612 *** sum(i, (Jacobian value var[i])*(Variable value var[i])) - rhs[i]
613 *** which is the real rhs.
614 *** However, this will _not_ be valid for the nonlinear case.
615 *** Other than this, the mps file will be the linearization of a nonlinear
616 *** system at the
617 *** current point. If you are adding an MINLP feature,
618 *** you'll need to come up with a better way.
619 ***
620 ***
621 *** MPS matrix strucutre
622 *** v
623 *** min/max cx: u
624 *** Ax (<= = >=) b s
625 *** e
626 *** 1 d
627 ***
628 *** | |
629 *** | |
630 *** \ / \ /
631 ***
632 *** +- -+
633 *** 1 -> | |
634 *** | |
635 *** | A |
636 *** | |
637 *** rused -> | |
638 *** +- -+
639 ***
640 *** crow -> [ c ]
641 ***
642 ***
643 *** rused row of last incident relation
644 *** crow = rused + 1, row of cost vector
645 *** vused column of last incident variable
646 ***
647 *** cap = max(vused+1,rused), size of sparse square matrix
648 *** cap = N ---> row/column 0 to N-1 exist
649 ***
650 **/
651 {
652 mtx_matrix_t mtx; /* main data structure */
653 var_filter_t vfilter; /* checks for free incident vars in relman_diffs */
654 double time0; /* time of Jacobian calculation, among other things */
655 struct rel_relation **rp; /* relation pointer */
656 real64 *rhs; /* temporary pointer for our RHS data */
657
658 if(obj == NULL) { /* a little preflight checking */
659 FPRINTF(stderr,"ERROR: (slv6) calc_matrix\n");
660 FPRINTF(stderr," System must have an objective!\n");
661 return NULL;
662 }
663
664 time0=tm_cpu_time(); /* start timing */
665 s->calc_ok = TRUE; /* no errors yet */
666
667 mtx = mtx_create(); /* create 0 order matrix */
668 mtx_set_order(mtx,cap); /* adjust size of square matrix to size of cap
669 (cap set in presolve.) The relman_diffs
670 routine returns values in the matrix */
671 /* these routines don't return success/fail */
672
673 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
674 vfilter.matchvalue = (VAR_INCIDENT | VAR_ACTIVE);
675 /* vfilter.fixed = var_false;
676 vfilter.incident = var_true;
677 vfilter.in_block = var_ignore; */
678
679 /* want to save column of residuals as they come along from relman_diffs */
680 *rhs_orig = create_array(rused,real64);
681 if(*rhs_orig == NULL) { /* memory allocation failed */
682 FPRINTF(stderr,"ERROR: (slv6) calc_matrix\n");
683 FPRINTF(stderr," Memory allocation for right hand side failed!\n");
684 return NULL;
685 }
686 rhs = *rhs_orig;
687
688
689 /* note: the rhs array is the residual at the current point, not what we want!
690 see further comments at the start of this routine */
691
692 for( rp = rlist ; *rp != NULL ; ++rp ) {
693 /* fill out A matrix only for used elements */
694 if( inc_rel_filter(*rp) ) {
695 *rhs = relman_diffs(*rp,&vfilter,mtx);
696 /*calculate each row of A matrix here! */
697 rhs++;
698 if( ! calc_ok ) { /* error check each call, calc_ok is in calc.h */
699 s->calc_ok = FALSE; /* error in diffs ! */
700 FPRINTF(stderr,"ERROR: (slv6) calc_matrix\n");
701 FPRINTF(stderr," Error in calculating A matrix.\n");
702 destroy_array(rhs_orig); /* clean up house, then die */
703 mtx_destroy(mtx); /* zap all alocated memory */
704 return NULL;
705 }
706 }
707 }
708 /* Calculate the rank of the matrix, before we add extra rows/cols */
709 mtx_output_assign(mtx, crow, vused);
710 if(! mtx_output_assigned(mtx)) { /* output assignment failed */
711 FPRINTF(stderr,"ERROR: (slv6) calc_matrix\n");
712 FPRINTF(stderr," Output assignment to calculate rank of problem failed.\n");
713 mtx_destroy(mtx); /* zap all alocated memory */
714 destroy_array(rhs_orig); /* clean up house, then die */
715 return NULL;
716 }
717 *rank = mtx_symbolic_rank(mtx);
718
719 if( *rank < 0 ) {
720 FPRINTF(stderr,"ERROR: (slv6) calc_matrix\n");
721 FPRINTF(stderr," Symbolic rank calculation failed, matrix may be bad.\n");
722 return mtx;
723 }
724
725 /* calculate the c vector and save it to the matrix */
726 if( ! calc_c(mtx, mtx_org_to_row(mtx,crow), obj) ) {
727 s->calc_ok = FALSE; /* error in diffs ! */
728 FPRINTF(stderr,"ERROR: (slv6) calc_matrix\n");
729 FPRINTF(stderr," Error in calculating objective coefficients.\n");
730 mtx_destroy(mtx); /* commit suicide */
731 destroy_array(rhs_orig); /* clean up house, then die */
732 return NULL;
733 }
734
735 s->block.jactime = tm_cpu_time() - time0; /* set overall jacobian time */
736
737 return mtx;
738 }
739
740
741 static void real_rhs(mtx_matrix_t Ac_mtx, /* Matrix representation of problem */
742 char relopcol[], /* is it incident? */
743 struct var_variable **vlist, /* Variable list (NULL terminated) */
744 int32 rused, /* in: total number of relations used */
745 real64 rhs[]) /* out: rhs array origin */
746 /**
747 *** Takes the residuals stored in rhs, and converts them into the actual right
748 *** hand sides we want.
749 ***
750 *** Note: the residual stored in the rhs array is not the real right hand side.
751 *** The residual returned by the diffs call is just the value of
752 ** (lhs expr) - (rhs expr)
753 *** we want the residual excluding the current variables. At the moment there isn't a
754 *** clean way to do this. It will be calculated in the real_rhs routine. This routine
755 *** will take for each relation:
756 *** sum(i, (Jacobian value var[i])*(Variable value var[i])) - rhs[i]
757 *** which is the real rhs. However, this will _not_ be valid for the nonlinear case.
758 *** Other than this, the mps file will be the linearization of a nonlinear system at the
759 *** current point. If you are adding an MINLP feature, you'll need to come up with a better way.
760 ***
761 **/
762 {
763 real64 a; /* value of mtx element */
764 mtx_coord_t nz; /* coordinate of row/column in A matrix */
765 mtx_range_t range; /* storage for range of A matrix, run down a column */
766 int32 currow; /* counter for current row */
767 int orgrow; /* original row number */
768 int orgcol; /* orignal col number */
769 double rowval; /* the sum of a[i]*x[i] in the row */
770
771 if(rhs == NULL) { /* a little preflight checking */
772 FPRINTF(stderr,"ERROR: (slv6) real_rhs\n");
773 FPRINTF(stderr," The routine was passed a NULL rhs pointer!\n");
774 return;
775 }
776
777 for(currow = 0; currow < rused; currow++) { /* loop over all rows, is _current_ column number */
778 orgrow = mtx_row_to_org(Ac_mtx, currow);
779 if (relopcol[orgrow] != rel_TOK_nonincident) { /* if it is incident row */
780
781 nz.col = mtx_FIRST; /* first nonzero col */
782 nz.row = currow; /* current row */
783 rowval = 0.0; /* accumulate value here */
784
785 a = mtx_next_in_row(Ac_mtx,&nz,mtx_range(&range,0,rused));
786
787 do { orgcol = mtx_col_to_org(Ac_mtx, nz.col);
788 rowval += a*var_value(*(vlist+orgcol));
789 a = mtx_next_in_row(Ac_mtx,&nz,mtx_range(&range,0,rused));
790
791 } while (nz.col != mtx_LAST);
792
793 rhs[orgrow] = rowval - rhs[orgrow]; /* set real value of right hand side */
794
795 }
796 }
797
798 }
799
800 /* _________________________________________________________________________ */
801
802 /**
803 *** Routines used by presolve
804 *** --------------------
805 *** insure_bounds - fix inconsistent bounds
806 *** update_vlist - add vars to vlist
807 *** determine_vlist - build new vlist
808 **/
809
810
811 static void insure_bounds(FILE *mif,slv6_system_t sys, struct var_variable *var)
812 /**
813 *** Insures that the variable value is within its bounds.
814 **/
815 {
816 real64 val,low,high;
817
818 low = var_lower_bound(var);
819 high = var_upper_bound(var);
820 val = var_value(var);
821 if( low > high ) {
822 FPRINTF(mif,"Bounds for variable ");
823 slv_print_var_name(mif,sys->slv,var);
824 FPRINTF(mif," are inconsistent [%g,%g].\n",low,high);
825 FPRINTF(mif,"Bounds will be swapped.\n");
826 var_set_upper_bound(var, low);
827 var_set_lower_bound(var, high);
828 low = var_lower_bound(var);
829 high = var_upper_bound(var);
830 }
831
832 if( low > val ) {
833 FPRINTF(mif,"Variable ");
834 slv_print_var_name(mif,sys->slv,var);
835 FPRINTF(mif," was initialized below its lower bound.\n");
836 FPRINTF(mif,"It will be moved to its lower bound.\n");
837 var_set_value(var, low);
838 } else if( val > high ) {
839 FPRINTF(mif,"Variable ");
840 slv_print_var_name(mif,sys->slv,var);
841 FPRINTF(mif," was initialized above its upper bound.\n");
842 FPRINTF(mif,"It will be moved to its upper bound.\n");
843 var_set_value(var, high);
844 }
845 }
846
847 #ifndef KILL
848
849 static struct var_variable **update_vlist(struct var_variable * *vlist, expr_t expr)
850 /**
851 *** Updates vlist, adding variables incident on given expression. The
852 *** old list is destroyed and the new list is returned.
853 **/
854 {
855 struct var_variable **newlist,**elist;
856 long nelts;
857
858 elist = expr_incidence_list(expr,NULL);
859 if( vlist == NULL ) return(elist);
860 if( (nelts=pl_length(elist)) == 0 ) {
861 ascfree( (POINTER)elist );
862 return(vlist);
863 }
864
865 nelts += pl_length(vlist) + 1;
866 newlist = (struct var_variable **)ascmalloc( sizeof(struct var_variable *) * (int)nelts );
867 pl_merge_0(newlist,vlist,elist);
868 ascfree( (POINTER)vlist );
869 ascfree( (POINTER)elist );
870 return(newlist);
871 }
872
873 #endif
874
875 #ifndef KILL
876
877 static void determine_vlist(slv6_system_t sys)
878 /**
879 *** This function constructs the variable list from the relation list. It
880 *** is assumed that sys->vlist_user == NULL so that this is necessary.
881 **/
882 {
883 bnd_boundary_t *bp;
884 struct rel_relation **rp;
885
886 if( pl_length(sys->vlist) > 0 )
887 ascfree( (POINTER)sys->vlist );
888 sys->vlist = NULL;
889 for( bp=sys->blist ; *bp != NULL ; ++bp ) {
890 sys->vlist = update_vlist(sys->vlist,bnd_lhs(*bp));
891 sys->vlist = update_vlist(sys->vlist,bnd_rhs(*bp));
892 }
893 for( rp=sys->rlist ; *rp != NULL ; ++rp ) {
894 sys->vlist = update_vlist(sys->vlist,rel_lhs(*rp));
895 sys->vlist = update_vlist(sys->vlist,rel_rhs(*rp));
896 }
897 if( sys->obj )
898 sys->vlist = update_vlist(sys->vlist,sys->obj);
899
900 if( sys->vlist == NULL )
901 slv6_set_var_list(sys,NULL);
902 }
903
904 #endif
905
906 /* _________________________________________________________________________ */
907
908 /**
909 *** External routines used from slv0 without modificiation
910 ***
911 *** slv6_set_var_list(sys,vlist)
912 *** slv6_get_var_list(sys)
913 *** slv6_set_bnd_list(sys,blist)
914 *** slv6_get_bnd_list(sys)
915 *** slv6_set_rel_list(sys,rlist)
916 *** slv6_get_rel_list(sys)
917 *** slv6_set_extrel_list(sys,erlist)
918 *** slv6_get_extrel_list(sys)
919 *** slv6_count_vars(sys,vfilter)
920 *** slv6_count_bnds(sys,bfilter)
921 *** slv6_count_rels(sys,rfilter)
922 *** slv6_set_obj_function(sys,obj)
923 *** slv6_get_obj_function(sys)
924 *** slv6_get_parameters(sys,parameters)
925 *** slv6_set_parameters(sys,parameters)
926 *** slv6_get_status(sys,status)
927 *** slv6_dump_internals(sys,level)
928 **/
929
930
931 void slv6_set_var_list(sys,vlist)
932 slv6_system_t sys;
933 struct var_variable **vlist;
934 {
935 static struct var_variable *empty_list[] = {NULL};
936 check_system(sys);
937 if( sys->vlist_user == NULL )
938 if( sys->vlist != NULL && pl_length(sys->vlist) > 0 )
939 ascfree( (POINTER)(sys->vlist) );
940 sys->vlist_user = vlist;
941 sys->vlist = (vlist==NULL ? empty_list : vlist);
942 sys->s.ready_to_solve = FALSE;
943 }
944
945 struct var_variable **slv6_get_var_list(sys)
946 slv6_system_t sys;
947 {
948 check_system(sys);
949 return( sys->vlist_user );
950 }
951
952 void slv6_set_bnd_list(sys,blist)
953 slv6_system_t sys;
954 bnd_boundary_t *blist;
955 {
956 static bnd_boundary_t empty_list[] = {NULL};
957 check_system(sys);
958 sys->blist_user = blist;
959 sys->blist = (blist==NULL ? empty_list : blist);
960 sys->s.ready_to_solve = FALSE;
961 }
962
963 bnd_boundary_t *slv6_get_bnd_list(sys)
964 slv6_system_t sys;
965 {
966 check_system(sys);
967 return( sys->blist_user );
968 }
969
970 void slv6_set_rel_list(sys,rlist)
971 slv6_system_t sys;
972 struct rel_relation **rlist;
973 {
974 static struct rel_relation *empty_list[] = {NULL};
975 check_system(sys);
976 sys->rlist_user = rlist;
977 sys->rlist = (rlist==NULL ? empty_list : rlist);
978 sys->s.ready_to_solve = FALSE;
979 }
980
981 struct rel_relation **slv6_get_rel_list(sys)
982 slv6_system_t sys;
983 {
984 check_system(sys);
985 return( sys->rlist_user );
986 }
987
988 void slv6_set_extrel_list(sys,erlist)
989 slv6_system_t sys;
990 struct ExtRelCache **erlist;
991 {
992 static struct ExtRelCache *empty_list[] = {NULL};
993 check_system(sys);
994 sys->erlist_user = erlist;
995 sys->erlist = (erlist==NULL ? empty_list : erlist);
996 sys->s.ready_to_solve = FALSE;
997 }
998
999 struct ExtRelCache **slv6_get_extrel_list(sys)
1000 slv6_system_t sys;
1001 {
1002 check_system(sys);
1003 return( sys->erlist_user );
1004 }
1005
1006 int slv6_count_vars(sys,vfilter)
1007 slv6_system_t sys;
1008 var_filter_t *vfilter;
1009 {
1010 struct var_variable **vp;
1011 int32 count = 0;
1012 check_system(sys);
1013 for( vp=sys->vlist; *vp != NULL; vp++ )
1014 if( var_apply_filter(*vp,vfilter) ) ++count;
1015 return( count );
1016 }
1017
1018 int slv6_count_bnds(sys,bfilter)
1019 slv6_system_t sys;
1020 bnd_filter_t *bfilter;
1021 {
1022 bnd_boundary_t *bp;
1023 int32 count = 0;
1024 check_system(sys);
1025 for( bp=sys->blist; *bp != NULL; bp++ )
1026 if( bnd_apply_filter(*bp,bfilter) ) ++count;
1027 return( count );
1028 }
1029
1030 int slv6_count_rels(sys,rfilter)
1031 slv6_system_t sys;
1032 rel_filter_t *rfilter;
1033 {
1034 struct rel_relation **rp;
1035 int32 count = 0;
1036 check_system(sys);
1037 for( rp=sys->rlist; *rp != NULL; rp++ )
1038 if( rel_apply_filter(*rp,rfilter) ) ++count;
1039 return( count );
1040 }
1041
1042 void slv6_set_obj_relation(slv6_system_t sys,struct rel_relation *obj)
1043 {
1044 check_system(sys);
1045 sys->obj = obj;
1046 sys->s.ready_to_solve = FALSE;
1047 }
1048
1049 struct rel_relation *slv6_get_obj_relation(slv6_system_t sys)
1050 {
1051 check_system(sys);
1052 return(sys->obj);
1053 }
1054
1055 void slv6_get_parameters(sys,parameters)
1056 slv6_system_t sys;
1057 slv_parameters_t *parameters;
1058 {
1059 check_system(sys);
1060 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
1061 }
1062
1063 void slv6_set_parameters(sys,parameters)
1064 slv6_system_t sys;
1065 slv_parameters_t *parameters;
1066 {
1067 check_system(sys);
1068 if (parameters->whose==slv6_solver_number)
1069 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
1070 }
1071
1072 void slv6_get_status(sys,status)
1073 slv6_system_t sys;
1074 slv_status_t *status;
1075 {
1076 check_system(sys);
1077 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
1078 }
1079
1080 void slv6_dump_internals(sys,level)
1081 slv6_system_t sys;
1082 int level;
1083 {
1084 check_system(sys);
1085 if (level > 0) {
1086 FPRINTF(stderr,"ERROR: (slv6) slv6_dump_internals\n");
1087 FPRINTF(stderr," slv6 does not dump its internals.\n");
1088 }
1089 }
1090
1091
1092 /* _________________________________________________________________________ */
1093
1094 /**
1095 *** External routines with minor modifications
1096 *** -----------------
1097 *** slv6_get_linsol_sys(sys) just returns NULL & error msg
1098 *** slv6_change_basis just return FALSE & error msg
1099 **/
1100
1101
1102 linsol_system_t slv6_get_linsol_sys(sys) /* just returns NULL & error msg */
1103 slv6_system_t sys;
1104 {
1105
1106 /* In the MPS file maker, there is no linsol_sys to return.
1107 So I just write out an error message, and return a NULL pointer */
1108
1109 FPRINTF(stderr,"ERROR: (slv6) slv6_change_basis\n");
1110 FPRINTF(stderr," This solver does not support changing the basis.\n");
1111
1112 return NULL;
1113 }
1114
1115
1116 boolean slv6_change_basis(slv6_system_t sys,int32 var, mtx_range_t *rng){
1117
1118 /* In the MPS file maker, changing the basis doesn't make any sense.
1119 Nor, for that matter, is there a basis in the first place.
1120 So I just write out an error message, and return FALSE */
1121
1122 FPRINTF(stderr,"ERROR: (slv6) slv6_change_basis\n");
1123 FPRINTF(stderr," This solver does not support changing the basis.\n");
1124
1125 return FALSE;
1126 }
1127
1128
1129 /* _________________________________________________________________________ */
1130
1131 /**
1132 *** External routines unique to slv6 (Based on routines from slv0)
1133 *** -----------------
1134 *** slv6_create() added solver specific initialization
1135 *** slv6_destroy(sys) added solver specific dealocation
1136 *** slv6_eligible_solver(sys) see if solver can do the current problem
1137 *** slv6_presolve(sys) set up system and create matrix/vectors
1138 *** slv6_solve(sys) call MPS routines
1139 *** slv6_iterate(sys) just calls slv6_solve
1140 *** slv6_resolve(sys) just calls slv6_solve
1141 **/
1142
1143
1144 slv6_system_t slv6_create() /* added mps initialization */
1145 {
1146 slv6_system_t sys;
1147
1148 /*** This routine allocates memory and initializes all data structures
1149 *** It should be a good source of comments on the system parameters and
1150 *** status flags used in slv6
1151 **/
1152
1153 /*** Allocate main system memory ***/
1154
1155 sys = (slv6_system_t)ascmalloc( sizeof(struct slv6_system_structure) );
1156 mem_zero_byte_cast(sys,0,sizeof(struct slv6_system_structure));
1157 sys->integrity = OK;
1158
1159
1160 /*** Initialize system parameters ***/
1161
1162 sys->p.output.more_important = stdout; /* used in MIF macro */
1163 sys->p.output.less_important = NULL; /* used in LIF macro (which is not used) */
1164
1165 sys->p.tolerance.pivot = 0.1; /* these tolerances are never used */
1166 sys->p.tolerance.singular = 1e-12;
1167 sys->p.tolerance.feasible = 1e-8;
1168 sys->p.tolerance.stationary = 1e-8;
1169 sys->p.tolerance.termination = 1e-12;
1170 sys->p.time_limit = 1500.0; /* never used */
1171 sys->p.iteration_limit = 100; /* never used */
1172 sys->p.partition = FALSE; /* never used, but don't want partitioning */
1173 sys->p.ignore_bounds = FALSE; /* never used, but must satisfy bounds */
1174 sys->p.whose = slv6_solver_number; /* read in slv6_set_parameters */
1175 sys->p.rho = 1.0;
1176 sys->p.sp.iap=&(sys->iarray[0]); /* all defaults in iarray are 0 */
1177 sys->p.sp.rap=&(sys->rarray[0]); /* all defaults in rarray are 0 */
1178 sys->p.sp.cap=&(sys->carray[0]); /* all defaults in carray are NULL */
1179 sys->p.sp.vap=NULL; /* not currently used */
1180
1181
1182 /*** Initialize mps data structure ***/
1183
1184 sys->mps.Ac_mtx = NULL; /* set all pointers to NULL - will all be set in presolve */
1185 sys->mps.lbrow = NULL; /* all other data in mps structure is 0 */
1186 sys->mps.ubrow = NULL;
1187 sys->mps.bcol = NULL;
1188 sys->mps.typerow = NULL;
1189 sys->mps.relopcol = NULL;
1190
1191
1192 /*** Initialize status flags ***/
1193
1194 sys->s.over_defined = FALSE; /* set to (sys->mps.rinc > sys->mps.vinc) in slv6_presolve */
1195 sys->s.under_defined = FALSE; /* set to (sys->mps.rinc < sys->mps.vinc) in slv6_presolve */
1196 sys->s.struct_singular = FALSE; /* set to (sys->mps.rank < sys->mps.rinc) in slv6_presolve */
1197 sys->s.calc_ok = TRUE; /* set in calc_matrix (FALSE if error occurs with diffs calc) */
1198 sys->s.ok = TRUE; /* set to (sys->s.calc_ok && !sys->s.struct_singular) in slv6_presolve */
1199 sys->s.ready_to_solve = FALSE; /* set to (sys->.ok) after slv6_presolve,
1200 set FALSE after: slv6_set_var_list, slv6_set_bnd_list,
1201 slv6_set_rel_list, slv6_set_extrel_list, slv6_set_obj_function
1202 tested in slv6_solve */
1203 sys->s.converged = FALSE; /* set FALSE after slv6_presolve; set TRUE after slv6_solve */
1204 sys->s.diverged = FALSE; /* always FALSE, never used */
1205 sys->s.inconsistent = FALSE; /* always FALSE, never used */
1206 sys->s.iteration_limit_exceeded = FALSE; /* always FALSE, never used */
1207 sys->s.time_limit_exceeded = FALSE; /* always FALSE, never used */
1208
1209 sys->s.block.number_of = 1; /* always 1, just have 1 block */
1210 sys->s.block.current_block = 0; /* always 1, start in first and only block */
1211 sys->s.block.current_size = 0; /* set to sys->mps.vused in slv6_presolve */
1212 sys->s.block.previous_total_size = 0; /* always 0, never used */
1213
1214 /* same : */
1215 sys->s.block.iteration = 0; /* set to 0 after slv6_presolve; set to 1 after slv6_solve */
1216 sys->s.iteration = 0; /* set to 0 after slv6_presolve; set to 1 after slv6_solve */
1217
1218 /* same : */
1219 sys->s.block.cpu_elapsed = 0.0; /* set to time taken by slv6_presolve and slv6_solve */
1220 sys->s.cpu_elapsed = 0.0; /* set to time taken by slv6_presolve and slv6_solve */
1221
1222 sys->s.block.functime = 0.0; /* always 0.0 since no function evaluation, never used */
1223 sys->s.block.residual = 0.0; /* always 0.0 since not iterating, never used */
1224 sys->s.block.jactime = 0.0; /* calculated in slv6_presolve, time for jacobian eval */
1225
1226 sys->s.costsize = sys->s.block.number_of; /* just one cost block, which will be set in */
1227
1228 sys->s.cost=create_zero_array(sys->s.costsize,struct slv_block_cost); /* allocate memory */
1229
1230
1231 /* Note: the cost vars are equivalent to other sys->s.* vars
1232
1233 sys->s.cost->size = sys->s.block.current_size
1234 sys->s.cost->iterations = sys->s.block.iteration
1235 sys->s.cost->jacs = sys->s.block.iteration
1236 sys->s.cost->funcs = always 0 since no function evals needed
1237 sys->s.cost->time = sys->s.block.cpu_elapsed
1238 sys->s.cost->resid = 0.0 whatever this is ?
1239 sys->s.cost->functime = 0.0 since no function evals needed
1240 sys->s.cost->jactime = sys->s.block.jactime
1241
1242 */
1243
1244 return(sys);
1245 }
1246
1247
1248 int slv6_destroy(sys)
1249 slv6_system_t sys;
1250 {
1251 int i;
1252
1253 if (check_system(sys)) return 1;
1254 slv6_set_var_list(sys,(struct var_variable **)NULL);
1255 slv6_set_obj_function(sys,NULL);
1256 slv6_set_bnd_list(sys,NULL);
1257 slv6_set_rel_list(sys,NULL);
1258 slv6_set_extrel_list(sys,NULL);
1259 sys->integrity = DESTROYED;
1260 if (sys->s.cost) ascfree(sys->s.cost); /* deallocate cost array */
1261
1262 /* deallocate strings here */
1263 for (i=0; i< slv6_CA_SIZE; i++) {
1264 if (sys->p.sp.cap[i] != NULL) ascfree(sys->p.sp.cap[i]); /* deallocate old, if any */
1265 }
1266
1267 nuke_pointers(sys->mps); /* free memory, and set all pointers to NULL */
1268 ascfree( (POINTER)sys );
1269
1270
1271 return 0;
1272 }
1273
1274
1275 boolean slv6_eligible_solver(sys)
1276
1277 /*** The system must have a relation list and objective before
1278 *** slv6_eligible_solver will return true
1279 **/
1280
1281 slv6_system_t sys;
1282 {
1283 struct rel_relation **rp;
1284 var_filter_t vfilter;
1285
1286 check_system(sys);
1287 if( sys->rlist == NULL ) {
1288 FPRINTF(stderr,"ERROR: (slv6) slv6_eligible_solver\n");
1289 FPRINTF(stderr," Relation list was never set.\n");
1290 return (FALSE);
1291 }
1292 if( sys->obj == NULL ) {
1293 FPRINTF(stderr,"ERROR: (slv6) slv6_eligible_solver\n");
1294 FPRINTF(stderr," No objective in problem.\n");
1295 return (FALSE);
1296 }
1297
1298 /* To Do: External Relations are currently being ingored. Is that proper?
1299 What if they're nonlinear */
1300
1301 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
1302 vfilter.matchvalue = (VAR_INCIDENT | VAR_ACTIVE);
1303 /*
1304 vfilter.fixed = var_false;
1305 vfilter.incident = var_true;
1306 vfilter.in_block = var_ignore; */
1307
1308 /* Check that the system is linear if iarray[SP6_NONLIN] == 0 */
1309 if (sys->iarray[SP6_NONLIN] == 0){
1310 for( rp=sys->rlist ; *rp != NULL ; ++rp ) /* check relations */
1311 if(!relman_is_linear(*rp,&vfilter)) {
1312 FPRINTF(MIF(sys), "ERROR: With the current settings, the MPS generator can only\n");
1313 FPRINTF(MIF(sys), " handle linear models. Nonlinearity in constraint:\n");
1314 slv_print_rel_name(MIF(sys),sys->slv, *rp);
1315 return(FALSE); /* don't do nonlinearities */
1316 }
1317 #ifndef KILL
1318 if (!exprman_is_linear(sys->obj,&vfilter)){
1319 FPRINTF(MIF(sys), "ERROR: With the current settings, the MPS generator can only\n");
1320 FPRINTF(MIF(sys), " handle linear models. Nonlinearity in objective.\n");
1321 return(FALSE); /* don't do nonlinearities */
1322 }
1323 #else
1324 FPRINTF(stderr," function slv6_elegible_solver is slv6 is broken.\n");
1325 FPRINTF(stderr," It calls exprman_is_linear with wrong number of args.\n");
1326 exit;
1327 #endif
1328 }
1329
1330 /* Note: initially I had this routine check to see if solver could handle
1331 binary, integer, semicontinuous, etc. Now I just convert types automatically.
1332 Binary vars become integer vars if a solver can do int but not binary.
1333 Semicontinuous vars become regular solver vars, with warnings about the
1334 conversion. If it can't handle integer vars, they are treated as regular
1335 solver vars, with warnings.
1336
1337 These conversions take place in the MPS.c file.*/
1338
1339 return TRUE;
1340 }
1341
1342 void slv6_presolve(sys)
1343 slv6_system_t sys;
1344 {
1345 struct var_variable **vp;
1346 struct rel_relation **rp;
1347 bnd_boundary_t *bp;
1348 int32 cap;
1349
1350 bnd_filter_t bfilter;
1351
1352 /* Check if necessary pointers are non-NULL */
1353 check_system(sys);
1354 if( sys->vlist == NULL ) {
1355 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1356 FPRINTF(stderr," Variable list was never set.\n");
1357 return;
1358 }
1359 if( sys->blist == NULL ) {
1360 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1361 FPRINTF(stderr," Boundary list was never set.\n");
1362 return;
1363 }
1364 if( sys->rlist == NULL ) {
1365 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1366 FPRINTF(stderr," Relation list was never set.\n");
1367 return;
1368 }
1369
1370 /* time presolve */
1371 sys->clock = tm_cpu_time(); /* record start time */
1372
1373 /* set up vlist, if necessary, and set all vars, rels, and boundary's
1374 to being nonincident, set up index scheme */
1375
1376 #ifndef KILL
1377 if( sys->vlist_user == NULL ) determine_vlist(sys);
1378 #else
1379 if( sys->vlist_user == NULL ){
1380 FPRINTF(stderr," function determine_vlist is slv6 is broken.\n");
1381 exit;
1382 }
1383 #endif
1384 sys->mps.cap = 0;
1385 for( vp=sys->vlist,cap=0 ; *vp != NULL ; ++vp ) {
1386 var_set_sindex(*vp,cap++);
1387 var_set_in_block(*vp,FALSE);
1388 }
1389 sys->mps.cap = cap;
1390 for( rp=sys->rlist,cap=0 ; *rp != NULL ; ++rp ) {
1391 rel_set_index(*rp,cap++);
1392 rel_set_in_block(*rp,FALSE);
1393 rel_set_satisfied(*rp,FALSE);
1394 }
1395 sys->mps.cap = MAX(sys->mps.cap,cap+1); /* allow an extra relation for crow,
1396 cap = N --> row/col 0 to N-1 exist */
1397 for( bp = sys->blist ; *bp != NULL ; ++bp ) {
1398 bnd_set_in_block(*bp,FALSE);
1399 bnd_set_active(*bp,FALSE);
1400 }
1401
1402 /**
1403 *** Now mark all variables appearing in the objective function,
1404 *** the boundaries and the relations as incident.
1405 **/
1406
1407 /* Mark all variables appearing in the objective function as incident */
1408 if( sys->obj ) exprman_decide_incidence_obj(sys->obj);
1409
1410 /* Set incidence of included vars in included bounds, calc bused, set all bounds inactive */
1411 sys->mps.bused = 0;
1412 bfilter.included = bnd_true;
1413 bfilter.in_block = bnd_ignore;
1414 for( bp = sys->blist ; *bp != NULL ; bp++ ) {
1415 if( bnd_apply_filter(*bp,&bfilter) ) {
1416 bndman_decide_incidence(*bp); /* mark incident variables */
1417 sys->mps.bused++;
1418 }
1419 }
1420
1421 /* Count the incident relations in rused */
1422 sys->mps.rused = 0;
1423 sys->mps.rinc = 0;
1424 for( rp = sys->rlist ; *rp != NULL ; rp++ ) {
1425 rel_set_satisfied(*rp,FALSE);
1426 if( inc_rel_filter(*rp) ) {
1427 relman_decide_incidence(*rp); /* mark incident variables in included rels */
1428 sys->mps.rinc++;
1429 }
1430 sys->mps.rused++;
1431 }
1432
1433 /* compute info for variables */
1434 sys->mps.vused = 0; /* number starting at 0 */
1435 sys->mps.vinc = 0;
1436 for( vp = sys->vlist ; *vp != NULL ; vp++ ) {
1437 if( free_inc_var_filter(*vp) )
1438 sys->mps.vinc++;
1439 sys->mps.vused++; /* count up incident, non-fixed vars */
1440 }
1441
1442 /* calculate values for other index_mps_t vars */
1443 sys->mps.crow = sys->mps.rused; /* note rused = N means rows 0 to N-1, exist,
1444 the next one will be numbered rused */
1445 /* calculate rank later */
1446
1447 /* Call slv6_elgibile_solver to see if the solver has a chance */
1448 /* If not bail now ... requires the incidence values of prev section be set */
1449 if(! slv6_eligible_solver(sys)) return;
1450
1451 /* Make sure that at least one incident variable and at least one incident
1452 relation exist, else bail */
1453 if ((sys->mps.rinc == 0) || (sys->mps.vinc == 0)) {
1454 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1455 FPRINTF(stderr," Your model must have at least one incident variable and equation.\n");
1456 FPRINTF(stderr," Incident variables: %d\n", sys->mps.vinc);
1457 FPRINTF(stderr," Incident equations: %d\n", sys->mps.rinc);
1458 return;
1459 }
1460
1461 /* free memory, and set all pointers to NULL */
1462 nuke_pointers(sys->mps);
1463
1464 /* setup matrix representaion of problem */
1465 sys->mps.Ac_mtx = calc_matrix(sys->mps.cap,
1466 sys->mps.rused,
1467 sys->mps.vused,
1468 sys->rlist,
1469 sys->obj,
1470 sys->mps.crow,
1471 &sys->s, /* how long for the jacobian calcs and any errors */
1472 &sys->mps.rank,
1473 &sys->mps.bcol);
1474 if( sys->mps.Ac_mtx == NULL ) {
1475 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1476 FPRINTF(stderr," Call to calc_matrix failed.\n");
1477 nuke_pointers(sys->mps);
1478 return;
1479 }
1480
1481 /* get upper bound row */
1482 sys->mps.ubrow = calc_bounds(sys->vlist, sys->mps.vinc, TRUE);
1483 if (sys->mps.ubrow == NULL) {
1484 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1485 FPRINTF(stderr," Error in calculating variable upper bounds.\n");
1486 nuke_pointers(sys->mps);
1487 return;
1488 }
1489
1490 /* get lower bound row */
1491 sys->mps.lbrow = calc_bounds(sys->vlist, sys->mps.vinc, FALSE);
1492 if (sys->mps.lbrow == NULL) {
1493 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1494 FPRINTF(stderr," Error in calculating variable lower bounds.\n");
1495 nuke_pointers(sys->mps);
1496 return;
1497 }
1498
1499 /* Call calc_svtlist to allocate array of variable types */
1500 sys->mps.typerow = calc_svtlist(sys->vlist,
1501 sys->mps.vinc,
1502 &sys->mps.solver_var_used, /* output */
1503 &sys->mps.solver_relaxed_used, /* output */
1504 &sys->mps.solver_int_used, /* output */
1505 &sys->mps.solver_binary_used, /* output */
1506 &sys->mps.solver_semi_used, /* output */
1507 &sys->mps.solver_other_used, /* output */
1508 &sys->mps.solver_fixed); /* output */
1509 if(sys->mps.typerow == NULL) { /* allocation failed */
1510 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1511 FPRINTF(stderr," Error in calculating the variable type list!\n");
1512 nuke_pointers(sys->mps);
1513 return;
1514 }
1515
1516 /* Call calc_reloplist here, to calculate the relational operators >=, <=, = */
1517 sys->mps.relopcol = calc_reloplist(sys->rlist, sys->mps.rinc);
1518 if(sys->mps.relopcol == NULL) { /* allocation failed */
1519 FPRINTF(stderr,"ERROR: (slv6) slv6_presolve\n");
1520 FPRINTF(stderr," Error in calculating the relational operators!\n");
1521 nuke_pointers(sys->mps);
1522 return;
1523 }
1524
1525 /* adjust the rhs vector so it actually contains the rhs */
1526 real_rhs(sys->mps.Ac_mtx, /* Matrix representation of problem */
1527 sys->mps.relopcol, /* is it incident? */
1528 sys->vlist, /* in: Variable list (NULL terminated) */
1529 sys->mps.rused, /* in: total number of relations used */
1530 sys->mps.bcol); /* out: rhs array origin */
1531
1532
1533 /* Call insure_bounds over all vars to make bounds self-consistent */
1534 for( vp=sys->vlist; *vp != NULL ; ++vp )
1535 insure_bounds(MIF(sys),sys, *vp);
1536
1537 /* Reset status flags */
1538 sys->s.over_defined = (sys->mps.rinc > sys->mps.vinc);
1539 sys->s.under_defined = (sys->mps.rinc < sys->mps.vinc);
1540 sys->s.struct_singular = (sys->mps.rank < sys->mps.rinc);
1541 sys->s.ok = sys->s.calc_ok && !sys->s.struct_singular;
1542 sys->s.ready_to_solve = sys->s.ok;
1543
1544 sys->s.converged = FALSE; /* changes to true after slv6_solve */
1545 sys->s.block.current_size = sys->mps.vused;
1546 sys->s.cost->size = sys->s.block.current_size;
1547
1548 sys->s.cpu_elapsed = (double)(tm_cpu_time() - sys->clock); /* record times */
1549 sys->s.block.cpu_elapsed = sys->s.cpu_elapsed;
1550 sys->s.cost->time = sys->s.cpu_elapsed;
1551 sys->s.cost->jactime = sys->s.block.jactime; /* from calc_matrix */
1552
1553 sys->s.block.iteration = 0; /* reset iteration "count", changes to 1 after slv6_solve */
1554 sys->s.iteration = 0;
1555 sys->s.cost->iterations = 0;
1556 sys->s.cost->jacs = 0;
1557
1558 }
1559
1560 void slv6_solve(sys)
1561 slv6_system_t sys;
1562 {
1563
1564 /* make sure none of the mps pointers are NULL */
1565 if ((sys->mps.Ac_mtx == NULL) ||
1566 (sys->mps.lbrow == NULL) ||
1567 (sys->mps.ubrow == NULL) ||
1568 (sys->mps.bcol == NULL) ||
1569 (sys->mps.typerow == NULL) ||
1570 (sys->mps.relopcol == NULL)) {
1571 FPRINTF(MIF(sys),"ERROR: Matrix representation of problem is not available.\n");
1572 FPRINTF(MIF(sys)," Perhaps the presolve routine was not called before slv6_solve.\n");
1573 return;
1574 }
1575
1576 /* Check system to see if it can be solved */
1577 check_system(sys);
1578 if( !sys->s.ready_to_solve ) {
1579 FPRINTF(stderr,"ERROR: (slv6) slv6_solve\n");
1580 FPRINTF(stderr," Not ready to solve.\n");
1581 return;
1582 }
1583
1584 sys->clock = tm_cpu_time(); /* record start time for solve */
1585
1586
1587 /* FPRINTF(MIF(sys),"_________________________________________\n");
1588 mtx_write_region_human(MIF(sys), sys->mps.Ac_mtx, mtx_ENTIRE_MATRIX);
1589 FPRINTF(MIF(sys),"_________________________________________\n");
1590 */
1591
1592 /* Call write_mps to create the mps file */
1593 write_MPS(sys->carray[SP6_FILENAME], /* filename for output */
1594 sys->mps, /* main chunk of data */
1595 sys->iarray, /* Integer subparameters */
1596 sys->rarray); /* Real subparameters */
1597
1598 /* replace .mps with .map at end of filename */
1599 *(sys->carray[SP6_FILENAME]+strlen(sys->carray[SP6_FILENAME])-2) = 'a';
1600 *(sys->carray[SP6_FILENAME]+strlen(sys->carray[SP6_FILENAME])-1) = 'p';
1601
1602 /* writes out a file mapping the CXXXXXXX variable names with the actual ASCEND names */
1603 write_name_map(sys->carray[SP6_FILENAME], /* user-specified filename */
1604 sys->vlist);
1605
1606
1607
1608 sys->s.cpu_elapsed += (double)(tm_cpu_time() - sys->clock);
1609 /* compute total elapsed time */
1610 sys->s.block.cpu_elapsed = sys->s.cpu_elapsed;
1611 sys->s.cost->time = sys->s.cpu_elapsed;
1612
1613 sys->s.converged = TRUE;
1614 sys->s.ready_to_solve = FALSE; /* !sys->s.converged */
1615
1616 sys->s.block.iteration = 1; /* change iteration "count", goes to 0 after slv6_presolve */
1617 sys->s.iteration = 1;
1618 sys->s.cost->iterations = 1;
1619 sys->s.cost->jacs = 1;
1620 sys->s.ready_to_solve = FALSE;
1621
1622 }
1623
1624
1625 void slv6_iterate(sys)
1626 slv6_system_t sys;
1627 {
1628 /* Writing an MPS file is a one shot deal. Thus, an interation
1629 is equivalent to solving the problem. So we just call
1630 slv6_solve */
1631
1632 check_system(sys);
1633 slv6_solve(sys);
1634
1635 }
1636
1637
1638 void slv6_resolve(sys)
1639 slv6_system_t sys;
1640 {
1641
1642 /* This routine is meant to be called when the following parts of
1643 the system change:
1644 - any parameter except "partition".
1645 - variable values.
1646 - variable nominal values.
1647 - variable bounds.
1648 However, if var values or bounds change, we need a new MPS file,
1649 so there is no way to use the previous solution.
1650 Just call slv6_solve, and do it the normal way.
1651 */
1652
1653 check_system(sys);
1654 slv6_solve(sys);
1655 }
1656
1657
1658 int slv6_register(SlvFunctionsT *sft)
1659 {
1660 if (sft==NULL) {
1661 FPRINTF(stderr,"slv6_register called with NULL pointer\n");
1662 return 1;
1663 }
1664
1665 sft->name = "makeMPS";
1666 sft->ccreate = slv6_create;
1667 sft->cdestroy = slv6_destroy;
1668 sft->celigible = slv6_eligible_solver;
1669 sft->getparam = slv6_get_parameters;
1670 sft->setparam = slv6_set_parameters;
1671 sft->getstatus = slv6_get_status;
1672 sft->solve = slv6_solve;
1673 sft->presolve = slv6_presolve;
1674 sft->iterate = slv6_iterate;
1675 sft->resolve = slv6_resolve;
1676 sft->getlinsol = slv6_get_linsol_sys;
1677 sft->getlinsys = slv6_get_linsolqr_sys;
1678 sft->getsysmtx = slv6_get_jacobian;
1679 sft->dumpinternals = slv6_dump_internals;
1680 return 0;
1681 }
1682 #endif /* #else clause of DYNAMIC_MPS */
1683 #endif /* #else clause of !STATIC_MPS && !DYNAMIC_MPS */

Properties

Name Value
svn:executable *

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22