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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 58663 byte(s)
Setting up web subdirectory in repository
1 aw0a 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