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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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