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

Annotation of /trunk/ascend4/solver/slv1.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: 83089 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 /*
2     * SLV: Ascend Numeric Solver
3     * by Karl Michael Westerberg
4     * Created: 2/6/90
5     * Version: $Revision: 1.30 $
6     * Version control file: $RCSfile: slv1.c,v $
7     * Date last modified: $Date: 2000/01/25 02:27:21 $
8     * Last modified by: $Author: ballan $
9     *
10     *
11     * This file is part of the SLV solver.
12     *
13     * Copyright (C) 1990 Karl Michael Westerberg
14     * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15     *
16     * The SLV solver is free software; you can redistribute
17     * it and/or modify it under the terms of the GNU General Public License as
18     * published by the Free Software Foundation; either version 2 of the
19     * License, or (at your option) any later version.
20     *
21     * The SLV solver is distributed in hope that it will be
22     * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
23     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24     * General Public License for more details.
25     *
26     * You should have received a copy of the GNU General Public License
27     * along with the program; if not, write to the Free Software Foundation,
28     * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
29     * COPYING. COPYING is found in ../compiler.
30     *
31     */
32    
33     /* known bugs:
34     * still uses pl_ functions.
35     * assumes old solver protocol.
36     * doesn't follow relman_eval signal trapping protocol.
37     */
38     #include <stdarg.h>
39     #include "utilities/ascConfig.h"
40     #include "utilities/ascSignal.h"
41     #include "utilities/ascMalloc.h"
42     #include "general/tm_time.h"
43     #include "utilities/set.h"
44     #include "utilities/mem.h"
45     #include "compiler/compiler.h"
46     #include "compiler/fractions.h"
47     #include "utilities/ascPanic.h"
48     #include "general/list.h"
49     #include "compiler/extfunc.h"
50     #include "compiler/dimen.h"
51     #include "compiler/types.h"
52     #include "compiler/find.h"
53     #include "compiler/relation_type.h"
54     #include "compiler/relation.h" /* relation enum */
55     #include "solver/mtx.h"
56     #include "solver/linsol.h"
57     #include "solver/linsolqr.h"
58     #include "solver/slv_types.h"
59     #include "solver/var.h"
60     #include "solver/rel.h"
61     #include "solver/discrete.h"
62     #include "solver/conditional.h"
63     #include "solver/logrel.h"
64     #include "solver/bnd.h"
65     #include "solver/calc.h"
66     #include "solver/relman.h"
67     #include "solver/slv_common.h"
68     #include "solver/slv_client.h"
69     #include "solver/slv1.h"
70     #if !defined(STATIC_MINOS) && !defined(DYNAMIC_MINOS)
71     /* do nothing */
72     int slv1_register(SlvFunctionsT *f)
73     {
74     (void)f; /* stop gcc whine about unused parameter */
75    
76     FPRINTF(stderr,"MINOS not compiled in this ASCEND IV.\n");
77     return 1;
78     }
79     #else /* either STATIC_MINOS or DYNAMIC_MINOS is defined */
80     #ifdef DYNAMIC_MINOS
81     /* do dynamic loading stuff. yeah, right */
82     #else /* following is used if STATIC_MINOS is defined */
83    
84     /********************************************************************\
85     * minos C interface
86     * ASCEND
87     * (C) Ben Allan, March 28, 1994
88     * $Revision: 1.30 $
89     * $Date: 2000/01/25 02:27:21 $
90     *
91     * MINOS 5.4 is proprietary software sitelicensed to Carnegie Mellon.
92     * Others who wish to use minos with ASCEND must get their own license
93     * and MINOS 5.4 sources. We provide only interface code to feed problems
94     * to MINOS 5.4.
95     *
96     * Notes: MINOS gets a problem scaled by the ascend nominals, not the
97     * straight equations. If you change nominals you must (p?)resolve again;
98     * minos will get very confused if you change nominals on the fly.
99     *
100     * minoss call assumptions:
101     * -nnobj=0 (obj linear) or nnobj=v.nonlinear.
102     * -all variables in a nonlinear rel/obj are nonlinear.( need better
103     * rel_linear primitive.)
104     *
105     \********************************************************************/
106    
107     /*********************************************************************
108     Ben Allan 5-8-94
109     The MINOS interface parameters, and it's relation to the slv0 shoebox:
110     (compare slv.h and MINOS 5.4 User's Guide Chapter 3, and App. A
111     Technical Report SOL 83-20R)
112    
113     OUTPUT.MORE_IMPORTANT:
114     OUTPUT.LESS_IMPORTANT: These will be used by the C interface and the
115     jacobian and function calls. The control of MINOS noise is done via
116     sub_parameter.
117     TOLERANCE.PIVOT: minos has no semantic equivalent.
118     TOLERANCE.SINGULAR: as minos PIVOT TOLERANCE
119     TOLERANCE.FEASIBLE: as minos ROW TOLERANCE & FEASIBILITY TOLERANCE
120     TOLERANCE.STATIONARY: as minos OPTIMALITY TOLERANCE
121     TOLERANCE.TERMINATION: no semantic equivalent.
122     TIME_LIMIT: use in C same as for slv0
123     ITERATION_LIMIT: as minos MAJOR ITERATIONS
124     PARTITION: minos has no semantic equivalent.
125     IGNORE_BOUNDS: minos disallows this entirely.
126     RHO: as minos PENALTY PARAMETER
127    
128     Subparameters implemented: (value/meaning)
129     [Except as noted, real/integer parameters given a value of 0 by the user will
130     default to MINOS values or last legal value set, since 0 is not generally
131     legal.]
132     sp.ia[SP1_COMPLETION] 0=>PARTIAL, 1=>FULL
133     sp.ia[SP1_MINITS] as minos MINOR ITERATIONS
134     sp.ia[SP1_CRASH] as minos CRASH OPTION (range 0-3)
135     sp.ia[SP1_DERIV] as minos DERIVATIVE LEVEL (range 0-3)
136     sp.ia[SP1_CFREQ] as minos CHECK FREQUENCY
137     sp.ia[SP1_FFREQ] as minos FACTORIZATION FREQUENCY
138     sp.ia[SP1_USELG] as minos LAGRANGIAN (0=>no, 1=> yes)
139     sp.ia[SP1_LFREQ] as minos LOG FREQUENCY
140     sp.ia[SP1_MULPR] as minos MULTIPLE PRICE
141     sp.ia[SP1_PARPR] as minos PARTIAL PRICE
142     sp.ia[SP1_JFLXB] as minos PRINT LEVEL (5 on/off bits, 0 legal)
143     sp.ia[SP1_SCALE] let minos SCALE (0 no, 1 yes, we handle scaling type)
144     sp.ia[SP1_SOLN] as minos SOLUTION except only (0=>no, 1=>yes)
145     sp.ia[SP1_PARAM] as minos SUPPRESS PARAMETERS (0 => no, 1=>yes)
146     sp.ia[SP1_VERIFY] as minos VERIFY LEVEL (range -1 - 3)
147     sp.ia[SP1_EFREQ] as minos EXPAND FREQUENCY
148     sp.ia[SP1_SUMMY] turn on unit 6, if not on, no print options will work.
149     sp.ia[SP1_FSUMY] turn on unit 9 for summary. minos.summary appears in pwd
150     sp.ia[SP1_LCONS] 0 => check for linearity, 1 => assume nonlinearity
151    
152     sp.ra[SP1_DAMP] as minos MAJOR DAMPING PARAMETER
153     sp.ra[SP1_FDIFF] as minos DIFFERENCE INTERVAL
154     sp.ra[SP1_CDIFF] as minos CENTRAL DIFFERENCE INTERVAL
155     sp.ra[SP1_FPREC] as minos FUNCTION PRECISION
156     sp.ra[SP1_LSTOL] as minos LINE SEARCH TOLERANCE
157     sp.ra[SP1_LUFTO] as minos LU FACTOR TOLERANCE (>=1.0)
158     sp.ra[SP1_LUUTO] as minos LU UPDATE TOLERANCE (>=1.0)
159     sp.ra[SP1_RADIUS] as minos RADIUS OF CONVERGENCE
160     sp.ra[SP1_SUBSP] as minos SUBSPACE TOLERANCE (range 0.0 - 1.0)
161     sp.ra[SP1_OBJLIM] as minos UNBOUNDED OBJECTIVE VALUE
162     sp.ra[SP1_STEPLM] as minos UNBOUNDED STEP SIZE
163     sp.ra[SP1_LOBJWT] as minos WEIGHT ON LINEAR OBJECTIVE
164     sp.ra[SP1_LUDTO] as minos LU DENSITY TOLERANCE (range 0.0 - 1.0)
165     sp.ra[SP1_LUSTO] as minos LU SINGULARITY TOLERANCE
166     sp.ra[SP1_LUWTO] as minos LU SWAP TOLERANCE
167     sp.ra[SP1_MINDAMP] as minos MINOR DAMPING PARAMETER
168    
169     Notes:
170     minos ITERATIONS will be set to ITERATION_LIMIT * sp.ia[SP1_MINITS]
171    
172     Status flags not implemented:
173     sys->s.:
174     over_defined
175     under_defined
176     struct_singular
177     Since minos adds variables to the system (slacks) to handle inequalities,
178     the ASCEND dof analysis is inaccurate. This could be fixed by anyone who
179     is good at counting beans in C.
180     *********************************************************************/
181    
182     #define KILL 0
183     #define slv1_solver_name "MINOS" /* Solver's name */
184     #define slv1_solver_number 1 /* Solver's number */
185    
186     #define SLV1(s) ((slv1_system_t)(s))
187    
188     #define slv1_IA_SIZE 19
189     #define slv1_RA_SIZE 16
190     #define slv1_CA_SIZE 0
191     #define slv1_VA_SIZE 0
192     /* subscripts for ia */
193     #define SP1_COMPLETION 0
194     #define SP1_MINITS 1
195     #define SP1_CRASH 2
196     #define SP1_DERIV 3
197     #define SP1_CFREQ 4
198     #define SP1_FFREQ 5
199     #define SP1_USELG 6
200     #define SP1_LFREQ 7
201     #define SP1_MULPR 8
202     #define SP1_PARPR 9
203     #define SP1_JFLXB 10
204     #define SP1_SCALE 11
205     #define SP1_SOLN 12
206     #define SP1_PARAM 13
207     #define SP1_VERIFY 14
208     #define SP1_EFREQ 15
209     #define SP1_SUMMY 16
210     #define SP1_FSUMY 17
211     #define SP1_LCONS 18
212     /* subscripts for ra. */
213     #define SP1_DAMP 0
214     #define SP1_FDIFF 1
215     #define SP1_CDIFF 2
216     #define SP1_FPREC 3
217     #define SP1_LSTOL 4
218     #define SP1_LUFTO 5
219     #define SP1_LUUTO 6
220     #define SP1_RADIUS 7
221     #define SP1_SUBSP 8
222     #define SP1_OBJLIM 9
223     #define SP1_STEPLM 10
224     #define SP1_LOBJWT 11
225     #define SP1_MINDAMP 12
226     #define SP1_LUDTO 13
227     #define SP1_LUSTO 14
228     #define SP1_LUWTO 15
229     /* subscripts for ca */
230     /* subscripts for va */
231    
232     /**********************************\
233     * NOUNDERBARS --> FORTRAN compiler *
234     * naming convention for subroutine *
235     * is wierd. *
236     * CRAY is treated as special case *
237     \**********************************/
238     #ifdef APOLLO
239     #define NOUNDERBARS TRUE
240     #endif
241     #ifdef _HPUX_SOURCE
242     #define NOUNDERBARS TRUE
243     #endif
244    
245    
246     #ifdef NOUNDERBARS
247     #define FUNCON_SUB_NAME funcon
248     #define FUNOBJ_SUB_NAME funobj
249     #define MATMOD_SUB_NAME matmod
250     #define MISPEC mispec
251     #define MIOPT miopt
252     #define MIOPTR mioptr
253     #define MIOPTI miopti
254     #define MINOSS minoss
255     #define GETCOMMON get_minos_common
256     #else
257     #define FUNCON_SUB_NAME funcon_
258     #define FUNOBJ_SUB_NAME funobj_
259     #define MATMOD_SUB_NAME matmod_
260     #define MISPEC mispec_
261     #define MIOPT miopt_
262     #define MIOPTR mioptr_
263     #define MIOPTI miopti_
264     #define MINOSS minoss_
265     #define GETCOMMON get_minos_common_
266     #endif
267    
268     /*
269     * Linux sources are typically f2c'ed. The
270     * under_bar in get_minos_common causes f2c to add
271     * another !!.
272     */
273     #ifdef linux
274     #undef GETCOMMON
275     #define GETCOMMON get_minos_common__
276     #endif
277    
278     /* CRAY compiler are a warped puppy: */
279     #ifdef CRAY
280     #define FUNCON_SUB_NAME FUNCON
281     #define FUNOBJ_SUB_NAME FUNOBJ
282     #define MATMOD_SUB_NAME MATMOD
283     #define MISPEC MISPEC
284     #define MIOPT MIOPT
285     #define MIOPTR MIOPTR
286     #define MIOPTI MIOPTI
287     #define MINOSS MINOSS
288     #define GETCOMMON GET_MINOS_COMMON
289     #endif
290    
291     /* ALL MINOS ints below are expected to be 32 bits (INTEGER*4) */
292     /* len on char * are length expected not counting null */
293     /* see minos docs and sources for details. */
294     extern MISPEC( int *, /* ispecs */
295     int *, /* iprint */
296     int *, /* isumm */
297     int *, /* nwcore */
298     int * /* inform */
299     );
300     extern MIOPT ( char *, /* buffer len 72 */
301     int *, /* iprint */
302     int *, /* isumm */
303     int * /* inform */
304     );
305     extern MIOPTR( char *, /* buffer len 56 */
306     double *, /* rvalue */
307     int *, /* iprint */
308     int *, /* isumm */
309     int * /* inform */
310     );
311     extern MIOPTI( char *, /* buffer len 56 */
312     int *, /* ivalue */
313     int *, /* iprint */
314     int *, /* isumm */
315     int * /* inform */
316     );
317     extern MINOSS( char *, /* start len 10 */
318     int *, /* m */
319     int *, /* n */
320     int *, /* nb */
321     int *, /* ne */
322     int *, /* nname */
323     int *, /* nncon */
324     int *, /* nnobj */
325     int *, /* nnjac */
326     int *, /* iobj */
327     double *, /* objadd */
328     char *, /* names len 40 */
329     double *, /* a */
330     int *, /* ha */
331     int *, /* ka */
332     double *, /* bl */
333     double *, /* bu */
334     int *, /* name1 */
335     int *, /* name2 */
336     int *, /* hs */
337     double *, /* xn */
338     double *, /* pi */
339     double *, /* rc */
340     int *, /* inform */
341     int *, /* mincor */
342     int *, /* ns */
343     int *, /* ninf */
344     double *, /* sinf */
345     double *, /* obj */
346     double *, /* z */
347     int * /* nwcore */
348     );
349    
350     extern GETCOMMON(int * /* major iterations */, int * /* major iterations */);
351    
352     #define MINOS_DEBUG FALSE
353     /* if MINOS_DEBUG then some sanity checking is done in calls from f77 to C */
354     #define D_ZERO (double)0.0
355     #define D_ONE (double)1.0
356    
357     struct count_data {
358     int32 used; /* Number of solver vars/included rels */
359     int32 nonlinear; /* Number of non-linear (used) vars/rels
360     non-linear vars are those nonlinear
361     in either rels or objective */
362     };
363    
364     struct jacobian_data {
365     linsol_system_t sys; /* Jacobian linear system */
366     mtx_matrix_t mtx; /* Matrix from funcon jacobian system */
367     mtx_matrix_t objmtx; /* Matrix from funobj */
368     real64 *rhs; /* RHS from jacobian system */
369     };
370     /* NULL ==> not created. sys, mtxs, and rhs created and destroyed together */
371    
372    
373     struct slv1_system_structure {
374     int integrity;
375     slv_parameters_t p; /* the standard slv parameters and status */
376     int iarray[slv1_IA_SIZE];
377     double rarray[slv1_RA_SIZE];
378     slv_status_t s;
379     int panic; /* bailout called noninteractively in FUNCON/OBJ */
380     double clock; /* cumulative cpu */
381    
382     /* Problem definition */
383     slv_system_t slv; /* slv_system_t back-link */
384     struct rel_relation *obj; /* Objective function: NULL = none */
385     struct var_variable **vlist; /* Variable list (NULL terminated) */
386     struct var_variable **vlist_user; /* User vlist (NULL = determine) */
387     struct rel_relation **rlist; /* Relation list (NULL terminated) */
388     struct rel_relation **rlist_user; /* User rlist (NULL = none) */
389    
390     struct count_data v,r; /* Variable and relation data */
391     int32 maxndx; /* Maximum index (capacity) */
392     int32 mrows; /* minos idea of r.used. note that
393     mrows =1 even if r.used =0 */
394    
395     /* Arrays allocated at presolve */
396     struct jacobian_data jacobian;
397     mtx_coord_t *nzlist; /* Complete list of jacobian elements, C indexed */
398     real64
399     *z, /* minos workspace, includes hessian, etc */
400     *rc, /* reduced costs as known to minos */
401     *pi, /* lagrange multipliers as known to minos for nl constraints */
402     *xn, /* variable values (scaled) */
403     *bl, /* lower bound values (scaled) */
404     *bu, /* upper bound values (scaled) */
405     *a; /* sparse jacobian (scaled) */
406     int32
407     *ha, /* row fortran indices for a */
408     *ka, /* column fortran indices for a, ha */
409     *hs; /* variable status (type) vector */
410    
411     /* scalars */
412     int32 itcount; /* Iteration countdown */
413     boolean basis; /* basis exists */
414     /* minos parms */
415     /* fortran files normally 5=stdin, 6=stdout, 0=no file, 9=`pwd`/minos.summary*/
416     int32 ispecs, /* fortran input unit number */
417     iprint, /* fortran output unit number */
418     isumm, /* fortran output unit number */
419     inform, /* return flag from minos */
420     nb, /* # free vars + # included rels */
421     nwcore, /* minos workspace size */
422     njac, /* # jacobian elements (linear+nonlinear) */
423     nnobj; /* number of nonlinear vars in objective */
424     };
425    
426     static slv1_system_t g_sys; /* Current system: used by FUNOBJ & FUNCON */
427    
428    
429     /**
430     *** Integrity checks
431     *** ----------------
432     *** check_system(sys)
433     **/
434    
435     #define OK ((int)491031596)
436     #define DESTROYED ((int)729104829)
437     /* note deaddead would be int 3735936685 */
438     static int check_system(slv1_system_t sys)
439     /* Checks sys for NULL and for integrity. */
440     {
441    
442     if( sys == NULL ) {
443     FPRINTF(stderr,"ERROR: (slv1) check_system\n");
444     FPRINTF(stderr," NULL system handle.\n");
445     return 1;
446     }
447    
448     switch( sys->integrity ) {
449     case OK:
450     return 0;
451     case DESTROYED:
452     FPRINTF(stderr,"ERROR: (slv1) check_system\n");
453     FPRINTF(stderr," System was recently destroyed.\n");
454     return 1;
455     default:
456     FPRINTF(stderr,"ERROR: (slv1) check_system\n");
457     FPRINTF(stderr," System reused or never allocated.\n");
458     return 1;
459     }
460     }
461    
462     /*********************************************************************
463     *** Memory management routines ***
464     free_unless_null(ptr)
465     zero(arr,nelts,type)
466     alloc_vector(len)
467     copy_vector(from,too,len)
468     free_vector(vec)
469     make_jacobian(sys)
470     destroy_jacobian(sys)
471     *********************************************************************/
472    
473     static void free_unless_null(POINTER ptr)
474     {
475     if( ptr != NULL )
476     ascfree(ptr);
477     }
478    
479    
480     #define zero(arr,nelts,type) \
481     mem_zero_byte_cast((arr),0,(nelts)*sizeof(type))
482     /* Zeros an array of nelts objects, each having given type. */
483    
484     #define alloc_vector(len) ((real64 *)ascmalloc((len)*sizeof(real64)))
485     #define alloc_zero_vector(len) ((real64 *)asccalloc((len),sizeof(real64)))
486     #define copy_vector(from,too,len) \
487     mem_move_cast((from),(too),(len)*sizeof(real64))
488     #define free_vector(vec) \
489     free_unless_null((POINTER)(vec))
490    
491     static void destroy_jacobian(slv1_system_t sys)
492     /* destroys, deallocates jacobian (sys,mtx, all rhs vectors attached to sys)*/
493     {
494     check_system(sys);
495     if( sys->jacobian.sys ) {
496     int count = linsol_number_of_rhs(sys->jacobian.sys)-1;
497     for( ; count >= 0; count-- )
498     free_vector(linsol_get_rhs(sys->jacobian.sys,count));
499     mtx_destroy(linsol_get_matrix(sys->jacobian.sys));
500     mtx_destroy(sys->jacobian.objmtx);
501     linsol_destroy(sys->jacobian.sys);
502     sys->jacobian.sys=NULL;
503     sys->jacobian.mtx=NULL;
504     sys->jacobian.objmtx=NULL;
505     sys->jacobian.rhs=NULL;
506     }
507     }
508    
509     static void make_jacobian(slv1_system_t sys)
510     /**
511     *** fills linsol, mtx, and one rhs for size sys->maxndx,
512     *** which is assumed set.
513     **/
514     {
515     destroy_jacobian(sys);
516     sys->jacobian.sys = linsol_create();
517     sys->jacobian.mtx = mtx_create();
518     sys->jacobian.objmtx = mtx_create();
519     mtx_set_order (sys->jacobian.mtx, sys->maxndx );
520     mtx_set_order (sys->jacobian.objmtx, sys->maxndx );
521     sys->jacobian.rhs = alloc_vector(sys->maxndx);
522     linsol_set_matrix(sys->jacobian.sys,sys->jacobian.mtx);
523     linsol_add_rhs(sys->jacobian.sys,sys->jacobian.rhs,0);
524     }
525     /*********************************************************************
526     *** General I/O routines ***
527     macros:
528     print_var_name(out,sys,var)
529     print_rel_name(out,sys,rel)
530     *********************************************************************/
531    
532     static void print_vector (slv1_system_t sys, double * vec,
533     int len, char * name)
534     {
535     int i;
536     for (i=0; i<len; i++)
537     FPRINTF(LIF(sys),"%s(%d) == %20.15g\n",name,i,vec[i]);
538     FFLUSH(LIF(sys));
539     }
540    
541     static void print_ivector (slv1_system_t sys, int * vec,int len, char * name) {
542     int i;
543     for (i=0; i<len; i++)
544     FPRINTF(LIF(sys),"%s(%d) == %d\n",name,i,vec[i]);
545     FFLUSH(LIF(sys));
546     }
547    
548     #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
549     #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
550    
551     static void print_output(FILE *out,slv1_system_t sys)
552     {
553     struct rel_relation **rp;
554     struct var_variable **vp;
555     int mrows, nvars,c;
556     real64 low,val,high,multiplier;
557     int32 index;
558    
559     c=0;
560     nvars = sys->v.used;
561     FPRINTF(out,"%-6s %-12s %-12s %-12s %-12s %-12s\n",
562     "INDEX","LOWER","LEVEL","UPPER","MULTIPLIER","NAME");
563     for (rp = sys->rlist;*rp != NULL; rp++) {
564     if (rel_included(*rp) && rel_active(*rp) ) {
565     index = rel_sindex(*rp);
566     low = sys->bl[nvars+c];
567     high = sys->bu[nvars+c];
568     val = rel_residual(*rp);
569     multiplier = sys->pi[c]; /* rel_multiplier(*rp); */
570     FPRINTF(out," % -6d % -8.4e % -8.4e %- 8.4e %- 8.4e ",
571     index,low,val,high,multiplier);
572     print_rel_name(out,sys,*rp);
573     PUTC('\n',out);
574     c++;
575     }
576     }
577     }
578    
579     /*****/
580    
581     static void make_nominal_positive(slv1_system_t sys,struct var_variable *var)
582     /* Makes nominal value of var > 0 */
583     /* if 0, make it var value unless value 0, in which case, make it 1 */
584     {
585     real64 n = var_nominal(var);
586    
587     if( n <= D_ZERO ) {
588     FILE *fp = MIF(sys);
589     if( n == D_ZERO ) {
590     if ( (n=fabs(var_value(var))) > D_ZERO)
591     var_set_nominal(var,n);
592     else
593     var_set_nominal(var,n = D_ONE);
594     FPRINTF(fp,"ERROR: (slv1) make_nominal_positive\n");
595     FPRINTF(fp," Variable ");
596     print_var_name(fp,sys,var);
597     FPRINTF(fp," \nhas nominal value of zero.\n");
598     FPRINTF(fp," Resetting to %g\n",n);
599     } else {
600     n = -n;
601     FPRINTF(fp,"ERROR: (slv1) make_nominal_positive\n");
602     FPRINTF(fp," Variable ");
603     print_var_name(fp,sys,var);
604     FPRINTF(fp," \nhas negative nominal value.\n");
605     FPRINTF(fp," Resetting to %g.\n",n);
606     var_set_nominal(var,n);
607     }
608     }
609     }
610    
611     /*********************************************************************
612     iteration_begins(sys)
613     iteration_ends(sys)
614     *********************************************************************/
615    
616     static void iteration_begins(slv1_system_t sys)
617     /* Prepares sys for an iteration, increasing the iteration counts
618     and starting the clock. */
619     {
620     sys->clock = tm_cpu_time();
621     }
622    
623     static void iteration_ends(slv1_system_t sys)
624     /* Prepares sys for exiting an iteration, stopping the clock and recording
625     the cpu time, as well as updating the status. */
626     {
627     double cpu_elapsed; /* elapsed this iteration */
628     boolean unsuccessful;
629    
630     cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
631     sys->s.cpu_elapsed += cpu_elapsed;
632    
633     if( !sys->s.converged ) {
634     sys->s.block.cpu_elapsed += cpu_elapsed;
635     sys->s.time_limit_exceeded =
636     (sys->s.block.cpu_elapsed > sys->p.time_limit);
637     sys->s.iteration_limit_exceeded =
638     (sys->s.block.iteration > sys->p.iteration_limit);
639     }
640    
641     unsuccessful = sys->s.diverged ||
642     sys->s.inconsistent ||
643     sys->s.iteration_limit_exceeded ||
644     sys->panic ||
645     sys->s.time_limit_exceeded;
646     sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
647     sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
648    
649     }
650    
651     static void install_nlvars(slv1_system_t sys,real64 *values)
652     /* *values Indexed by column: scaled by var nominal */
653     /*********************************************************************
654     Moves nonlinear variables from given array to the value field of
655     corresponding free variable. unscale in process. (mult x*varnom ->var)
656     *********************************************************************/
657     {
658     int32 col;
659     struct var_variable *var;
660     for( col=0 ; col < sys->v.nonlinear ; ++col ) {
661     var = sys->vlist[mtx_col_to_org(sys->jacobian.mtx,col)];
662     /* var_set_value(var,values[col]*var_nominal(var)); */
663     var_set_value(var,values[col]);
664     }
665     }
666    
667     static void install_lvars(slv1_system_t sys,real64 *values)
668     /*********************************************************************
669     Moves linear variables from given array to the value field
670     of corresponding free variable. unscale in process. (mult x*varnom ->var)
671     *********************************************************************/
672     {
673     int32 col;
674     struct var_variable *var;
675     for( col=sys->v.nonlinear ; col < sys->v.used ; ++col ) {
676     var = sys->vlist[mtx_col_to_org(sys->jacobian.mtx,col)];
677     /* var_set_value(var,values[col]*var_nominal(var)); */
678     var_set_value(var,values[col]);
679     }
680     }
681    
682     static void install_allvars(slv1_system_t sys,real64 *values)
683     /*********************************************************************
684     Moves variables from given array to the value field of each free variable.
685     unscale in process (mult x*varnom ->var)
686     *********************************************************************/
687     {
688     int32 col;
689     struct var_variable *var;
690     for( col=0 ; col < sys->v.used ; ++col ) {
691     var = sys->vlist[mtx_col_to_org(sys->jacobian.mtx,col)];
692     /* var_set_value(var,values[col]*var_nominal(var)); */
693     var_set_value(var,values[col]);
694     }
695     }
696    
697     static void calc_residuals(slv1_system_t sys)
698     /* Residuals are calculated, norm value is stored in status, and
699     sys->s.calc_ok is set accordingly. Knows nothing about scaling */
700     {
701     real64 sum,res;
702     struct rel_relation **rp;
703    
704     calc_ok = TRUE;
705     sum = D_ZERO;
706     for( rp=sys->rlist ; *rp != NULL ; ++rp ) {
707     if (rel_included(*rp) && rel_active(*rp) ) {
708     res = relman_eval(*rp,&calc_ok);
709     if( !relman_calc_satisfied(*rp,sys->p.tolerance.feasible) ||
710     ( !rel_less(*rp) && !rel_greater(*rp) )
711     ) {
712     sum += calc_sqr_D0(res);
713     }
714     }
715     }
716     if (sum > D_ZERO)
717     sys->s.block.residual = calc_sqrt_D0(sum);
718     else sys->s.block.residual=sum;
719     /* no nonlin vars -> funcon, funobj not called and we must set convergence */
720     if (!(sys->v.nonlinear)) {
721     if (sys->inform==0) {
722     sys->s.converged=TRUE;
723     } else {
724     if (sys->inform !=3) {
725     sys->s.diverged=TRUE;
726     }
727     }
728     }
729    
730     if( sys->obj != NULL )
731     relman_eval(sys->obj,&calc_ok); /* Check for OK calculations */
732     sys->s.calc_ok = calc_ok;
733     }
734    
735     /*********************************************************************
736     *** MINOS subroutines ***
737     funobj(...)
738     funcon(...) int are assumed to be FORTRAN i4 in these functions.
739     *********************************************************************/
740    
741     void FUNOBJ_SUB_NAME(int *mode,
742     int *n,
743     double *x,
744     double *f,
745     double *g,
746     int *nstate,
747     int *nprob,
748     double *z,
749     int *nwcore)
750     /*********************************************************************
751     Computes the objective function and its gradient. Called by MINOS.
752     We should be scaling the obj value, I suspect.
753     IN
754     mode Calculate gradient iff mode == 2
755     n # of nonlinear obj variables (which is all nl in sys until smarter)
756     x Variable values (scaled) dimension n
757     nstate >=2 ==> final call. minos return will be nstate -2
758     nprob # of problems (so what?)
759     z minos workspace
760     nwcore minos workspace size
761    
762     OUT
763     mode -1 ==> stop! set if Solv_C_CheckHalt is true.
764     f Value of objective function
765     g Vector of gradients dimension n of objective
766     g scaled (mult by var_nominal)
767     *********************************************************************/
768     {
769     FPRINTF(LIF(g_sys),"FUNOBJ called with nstate = %d\n",*nstate);
770     /* unscale x and shove it back into ASCEND */
771     install_nlvars(g_sys,x);
772    
773     print_vector(g_sys,x,*n,"x");
774     if (*mode!=0)
775     print_vector(g_sys,g,*n,"g");
776     print_vector(g_sys,f,1,"f-in");
777     if( *nstate >= 2 ) {
778     if( *nstate == 2 )
779     g_sys->s.converged = TRUE;
780     else if( *nstate != 5 )
781     g_sys->s.diverged = TRUE;
782     return;
783     }
784     if (*n!=g_sys->v.nonlinear) {
785     FPRINTF(MIF(g_sys),"FUNOBJ called with confusing number of variables(%d)",
786     *n);
787     *mode=-1;
788     g_sys->panic=TRUE;
789     return;
790     }
791     calc_ok = TRUE;
792     if( (*mode == 2) ) {
793     mtx_coord_t nz;
794     real64 value;
795     var_filter_t vfilter;
796     struct var_variable *var;
797    
798     /* can't hurt, del(nothing) =0, shouldn't ever get here */
799     zero(g,*n,double);
800     if( g_sys->obj == NULL ) {
801     *f = D_ZERO;
802     return;
803     }
804    
805     /* calc free & incident dF/dx */
806     vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
807     vfilter.matchvalue = (VAR_INCIDENT )| VAR_ACTIVE;
808     /* vfilter.fixed = var_false;
809     vfilter.incident = var_true;
810     vfilter.in_block = var_ignore; */
811     nz.row=(int32)0;
812     /* this needs to change */
813     mtx_mult_row_zero(g_sys->jacobian.objmtx,nz.row,mtx_ALL_COLS);
814     /* discard status = */
815     relman_diffs(g_sys->obj, &vfilter, g_sys->jacobian.objmtx,f);
816     /* multiply gradient by nominal and put in minos gradient vector */
817     nz.col = mtx_FIRST;
818     while( value = mtx_next_in_row(g_sys->jacobian.objmtx,&nz,mtx_ALL_COLS),
819     nz.col != mtx_LAST ) {
820     var = g_sys->vlist[mtx_col_to_org(g_sys->jacobian.objmtx,nz.col)];
821     #if MINOS_DEBUG
822     if (nz.col >= *n) { /* take this bit out if it works all the time */
823     FPRINTF(MIF(g_sys),
824     "FUNOBJ stuffing a confused objective gradient(%d)", *n);
825     *mode=-1;
826     g_sys->panic=TRUE;
827     return;
828     }
829     #endif
830     /* g[nz.col] = value * var_nominal(var); */
831     g[nz.col] = value;
832     /* *((g_sys->obj->negate==TRUE)? D_ONE : D_ONE) */
833     FPRINTF(LIF(g_sys),"d(obj)/d(");
834     print_var_name(LIF(g_sys),g_sys,var);
835     FPRINTF(LIF(g_sys),")=%20.16g\n",g[nz.col]);
836     FPRINTF(LIF(g_sys),"g subscript stuffed: %d\n",nz.col);
837     } /* for */
838     /* end gradient */
839     } else {
840     *f = (g_sys->obj==NULL) ? D_ZERO : relman_eval(g_sys->obj,&calc_ok);
841     /* ((g_sys->obj->negate==TRUE)? D_ONE:D_ONE)*exprman_eval(g_sys->obj); */
842     }
843    
844     iteration_ends(g_sys); /* timout to check interface */
845     if ( Solv_C_CheckHalt() ) *mode=-1;
846     iteration_begins(g_sys);
847    
848     FPRINTF(LIF(g_sys),"FUNOBJ returning obj= %g\nGradient\n",*f);
849     if (*mode!=0)
850     print_vector(g_sys,g,*n,"gout");
851     else
852     FPRINTF(LIF(g_sys),"no gradient called for.\n");
853     if( !calc_ok ) {
854     FPRINTF(MIF(g_sys),"!!FUNOBJ: Warning: calculation error(s).\n");
855     }
856     }
857    
858     void FUNCON_SUB_NAME( int *mode,
859     int *m,
860     int *n,
861     int *njac,
862     double *x,
863     double *f,
864     double *g,
865     int *nstate,
866     int *nprob,
867     double *z,
868     int *nwcore)
869     /*********************************************************************
870     Computes the jacobian. Called by MINOS. Each call updates the
871     iteration counts.
872    
873     We really should be scaling the nonlinear residuals. this is lunacy.
874    
875     For now, don't do energy balances. :-(
876    
877     Since all variables in non-linear relations/obj are assumed to appear
878     non-linearly, this routine does not require updating the lvars to be
879     correct. It will need to update lvars somehow if we sort variables
880     within a nonlinear equation.
881    
882     IN
883     mode Calculate the gradient iff mode == 2
884     m # of nonlinear constraints (rows)
885     n # of nonlinear variables (cols)
886     njac # of varying jacobian elements as known by minos
887     x nonlinear variable values (index by column)
888     nstate >=2 ==> solution is optimal. (set converged true)
889     nprob # problems
890     z minos workspace
891     nwcore minos workspace size
892    
893     OUT
894     mode -1 ==> stop!
895     f Computed residuals of nonlinear relations (scaled?)
896     g Function gradients (scaled (mult by var_nominal) )
897     it seems to be the case that g is a sans lin col/row elements.
898     *********************************************************************/
899     {
900     int32 row,ndx,gndx,eqn_n;
901     struct rel_relation **rp;
902     rel_filter_t rfilter;
903     var_filter_t vfilter;
904     struct var_variable *var;
905    
906     FPRINTF(LIF(g_sys),"FUNCON was called with nstate = %d\n",*nstate);
907    
908     if (*n!=g_sys->v.nonlinear) {
909     FPRINTF(MIF(g_sys),"FUNCON called with confusing number of variables(%d)",
910     *n);
911     *mode=-1;
912     g_sys->panic=TRUE;
913     return;
914     }
915    
916     print_vector(g_sys,x,*n,"x");
917     if (*mode!=0)
918     print_vector(g_sys,g,*njac,"g");
919     print_vector(g_sys,f,*m,"f-in");
920    
921     if( *nstate >= 2 ) {
922     if( *nstate == 2 )
923     g_sys->s.converged = TRUE;
924     else
925     if( *nstate != 5 )
926     g_sys->s.diverged = TRUE;
927     return;
928     }
929    
930     install_nlvars(g_sys,x);
931     if( *nstate >= 2 )
932     return;
933    
934     rfilter.matchbits = (REL_INCLUDED | REL_ACTIVE);
935     rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
936     vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
937     vfilter.matchvalue = (VAR_INCIDENT | VAR_ACTIVE);
938     /* vfilter.fixed = var_false;
939     vfilter.incident = var_true;
940     vfilter.in_block = var_ignore; */
941     calc_ok = TRUE;
942    
943     if (*mode !=0) { /* skip jacobian if mode 0 */
944     /* don't zero the nl jacobian g as it also may have constants */
945     for( rp=g_sys->rlist; *rp != NULL; rp++ ) {
946     eqn_n = rel_sindex(*rp);
947     if( rel_apply_filter(*rp,&rfilter) ) {
948     /* this needs to change */
949     mtx_mult_row_zero(g_sys->jacobian.mtx,
950     mtx_org_to_row(g_sys->jacobian.mtx,eqn_n),
951     mtx_ALL_COLS);
952     /* discard status = */
953     relman_diffs(*rp,&vfilter,g_sys->jacobian.mtx,
954     &(g_sys->jacobian.rhs[eqn_n]));
955     }
956     }
957    
958     if( !calc_ok ) {
959     FPRINTF(MIF(g_sys),"!!FUNCON: Warning:jacobian calculation error(s).\n");
960     *mode=-1;
961     }
962    
963     FPRINTF(LIF(g_sys),"FUNCON jacobian and residuals calc-ed\n");
964     for( gndx=ndx=0 ; ndx < g_sys->njac ; ++ndx )
965     /* don't restuff lrows, just nl rows*/
966     if( g_sys->nzlist[ndx].row < g_sys->r.nonlinear ) {
967     var = g_sys->vlist[mtx_col_to_org(g_sys->jacobian.mtx,
968     g_sys->nzlist[ndx].col)];
969     #if MINOS_DEBUG
970     if (gndx>*njac){
971     FPRINTF(MIF(g_sys), "FUNCON stuffing a confused jacobian(%d)", *n);
972     g_sys->panic=TRUE;
973     *mode=-1;
974     return;
975     }
976     #endif
977     g[gndx++] = mtx_value(g_sys->jacobian.mtx , g_sys->nzlist+ndx);
978     /* g[gndx++] = var_nominal(var) * mtx_value(g_sys->jacobian.mtx ,
979     g_sys->nzlist+ndx); */
980     }
981     FPRINTF(LIF(g_sys),"FUNCON minos jacobian stuffed \n");
982     } else { /* calc functions only */
983     for( rp=g_sys->rlist; *rp != NULL; rp++ ) {
984     eqn_n = rel_sindex(*rp);
985     if( rel_apply_filter(*rp,&rfilter) ) {
986     g_sys->jacobian.rhs[eqn_n] = relman_eval(*rp,&calc_ok);
987     }
988     }
989     if( !calc_ok ) {
990     FPRINTF(MIF(g_sys),"!!FUNCON: Warning:jacobian calculation error(s).\n");
991     *mode=-1;
992     }
993     }
994     /* stuff function results */
995     zero(f,*m,double);
996     for( row=0 ; row < *m ; ++row )
997     f[row] = g_sys->jacobian.rhs[mtx_row_to_org(g_sys->jacobian.mtx,row)];
998    
999     iteration_ends(g_sys); /* timout to check interface */
1000     if ( Solv_C_CheckHalt() ) *mode=-1;
1001     iteration_begins(g_sys);
1002    
1003     print_vector(g_sys,x,*n,"x");
1004     if (*mode!=0)
1005     print_vector(g_sys,g,*njac,"g");
1006     print_vector(g_sys,f,*m,"f-out");
1007     }
1008    
1009     void MATMOD_SUB_NAME( /* Lots of irrelevent arguments */ )
1010     {
1011     FPRINTF(stderr,"ERROR: (slv1) MATMOD_SUB_NAME\n");
1012     FPRINTF(stderr," This should never be called.\n");
1013     }
1014    
1015     /*********************************************************************
1016     *** MINOS interface ***
1017     make_specs(sys)
1018     insure_bounds(sys,var)
1019     make_bounds(sys)
1020     invoke_minos(sys)
1021     *********************************************************************/
1022    
1023    
1024     #define THEIR_INFINITY (double)1.0e10
1025    
1026     static int make_specs(slv1_system_t sys)
1027     /* stuff the specs into minos */
1028     {
1029     int i=0,j=0,out=0;
1030     double db;
1031     char *s72, *s56;
1032     FILE * fp, *fpl;
1033     check_system(sys);
1034     fp=MIF(sys);
1035     fpl=LIF(sys);
1036     sys->inform=0;
1037     s72=(char *)ascmalloc((73*sizeof(char)));
1038     s56=(char *)ascmalloc((57*sizeof(char)));
1039     if (!s56 || !s72) {
1040     FPRINTF(stderr,"ascmalloc failed in make_specs!");
1041     return (1);
1042     }
1043    
1044     if (!sys->iarray[SP1_PARAM]) out=6;
1045     sprintf(s72,"%-72s","Defaults");
1046     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1047     if ((sys->inform)) {
1048     FPRINTF(fp,"minos: make_specs: Unable to reset Defaults.\n");
1049     sys->inform=0;
1050     j++;
1051     }
1052    
1053     sprintf(s72,"%-72s","Timing 0");
1054     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1055     if ((sys->inform)) {
1056     FPRINTF(fp,"minos: make_specs: Unable to kill clock.\n");
1057     sys->inform=0;
1058     j++;
1059     }
1060    
1061     if (sys->iarray[SP1_COMPLETION])
1062     sprintf(s72,"%-72s","Completion full");
1063     else
1064     sprintf(s72,"%-72s","Completion partial");
1065     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1066     if ((sys->inform)) {
1067     FPRINTF(fp,"minos: make_specs: Unable to set completion.\n");
1068     sys->inform=0;
1069     j++;
1070     }
1071    
1072     sprintf(s56,"%-56s","Crash option");
1073     MIOPTI(s56,&(sys->iarray[SP1_CRASH]),
1074     &out,&(sys->isumm), &(sys->inform));
1075     if ((sys->inform)) {
1076     FPRINTF(fp,"minos: make_specs: Unable to set crash option.\n");
1077     sys->inform=0;
1078     j++;
1079     }
1080    
1081     sprintf(s56,"%-56s","Derivative level");
1082     MIOPTI(s56,&(sys->iarray[SP1_DERIV]),
1083     &out,&(sys->isumm), &(sys->inform));
1084     if ((sys->inform)) {
1085     FPRINTF(fp,"minos: make_specs: Unable to set derivative level.\n");
1086     sys->inform=0;
1087     j++;
1088     }
1089    
1090     sprintf(s72,"%-72s","Jacobian Sparse");
1091     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1092     if ((sys->inform)) {
1093     FPRINTF(fp,"minos: make_specs: Unable to specify jacobian is sparse.\n");
1094     sys->inform=0;
1095     j++;
1096     }
1097    
1098     db=-THEIR_INFINITY;
1099     sprintf(s56,"%-56s","Lower bound");
1100     MIOPTR(s56,&db,
1101     &out,&(sys->isumm), &(sys->inform));
1102     if ((sys->inform)) {
1103     FPRINTF(fp,
1104     "minos: make_specs: Unable to set lower bound to -their_infinity.\n");
1105     sys->inform=0;
1106     j++;
1107     }
1108    
1109     i=(int)sys->p.iteration_limit;
1110     sprintf(s56,"%-56s","Major iterations");
1111     MIOPTI(s56,&i,
1112     &out,&(sys->isumm), &(sys->inform));
1113     if ((sys->inform)) {
1114     FPRINTF(fp,"minos: make_specs: Unable to set major iteration limit.\n");
1115     sys->inform=0;
1116     j++;
1117     }
1118    
1119     i=sys->iarray[SP1_MINITS];
1120     sprintf(s56,"%-56s","Minor iterations");
1121     MIOPTI(s56,&i,
1122     &out,&(sys->isumm), &(sys->inform));
1123     if ((sys->inform)) {
1124     FPRINTF(fp,"minos: make_specs: Unable to set minor iteration limit.\n");
1125     sys->inform=0;
1126     j++;
1127     }
1128    
1129     i=sys->iarray[SP1_MINITS] * (int)sys->p.iteration_limit;
1130     sprintf(s56,"%-56s","Iterations limit");
1131     MIOPTI(s56,&i,
1132     &out,&(sys->isumm), &(sys->inform));
1133     if ((sys->inform)) {
1134     FPRINTF(fp,"minos: make_specs: Unable to set Iterations limit.\n");
1135     sys->inform=0;
1136     j++;
1137     }
1138    
1139     if (sys->p.tolerance.stationary > D_ZERO) {
1140     sprintf(s56,"%-56s","Optimality tolerance");
1141     MIOPTR(s56,&(sys->p.tolerance.stationary),
1142     &out,&(sys->isumm), &(sys->inform));
1143     if ((sys->inform)) {
1144     FPRINTF(fp,"minos: make_specs: Unable to set Optimality tolerance.\n");
1145     sys->inform=0;
1146     j++;
1147     }
1148     }
1149    
1150     if (sys->p.tolerance.singular > D_ZERO) {
1151     sprintf(s56,"%-56s","Pivot tolerance");
1152     MIOPTR(s56,&(sys->p.tolerance.singular),
1153     &out,&(sys->isumm), &(sys->inform));
1154     if ((sys->inform)) {
1155     FPRINTF(fp,"minos: make_specs: Unable to set Pivot tolerance.\n");
1156     sys->inform=0;
1157     j++;
1158     }
1159     }
1160    
1161     if (sys->p.tolerance.feasible > D_ZERO) {
1162     sprintf(s56,"%-56s","Feasibility tolerance");
1163     MIOPTR(s56,&(sys->p.tolerance.feasible),
1164     &out,&(sys->isumm), &(sys->inform));
1165     if ((sys->inform)) {
1166     FPRINTF(fp,"minos: make_specs: Unable to set feasibility tolerance.\n");
1167     sys->inform=0;
1168     j++;
1169     }
1170     sprintf(s56,"%-56s","Row tolerance");
1171     MIOPTR(s56,&(sys->p.tolerance.feasible),
1172     &out,&(sys->isumm), &(sys->inform));
1173     if ((sys->inform)) {
1174     FPRINTF(fp,"minos: make_specs: Unable to set row tolerance.\n");
1175     sys->inform=0;
1176     j++;
1177     }
1178     }
1179    
1180     sprintf(s56,"%-56s","Major damping parameter");
1181     MIOPTR(s56,&(sys->rarray[SP1_DAMP]),
1182     &out,&(sys->isumm), &(sys->inform));
1183     if ((sys->inform)) {
1184     FPRINTF(fp,"minos: make_specs: Unable to set major damping parameter %g.\n",
1185     sys->rarray[SP1_DAMP]);
1186     sys->inform=0;
1187     j++;
1188     }
1189    
1190     sprintf(s56,"%-56s","Minor damping parameter");
1191     MIOPTR(s56,&(sys->rarray[SP1_MINDAMP]),
1192     &out,&(sys->isumm), &(sys->inform));
1193     if ((sys->inform)) {
1194     FPRINTF(fp,"minos: make_specs: Unable to set minor damping parameter %g.\n",
1195     sys->rarray[SP1_MINDAMP]);
1196     sys->inform=0;
1197     j++;
1198     }
1199    
1200     if (sys->rarray[SP1_FPREC] > D_ZERO) {
1201     sprintf(s56,"%-56s","Function precision");
1202     MIOPTR(s56,&(sys->rarray[SP1_FPREC]),
1203     &out,&(sys->isumm), &(sys->inform));
1204     if ((sys->inform)) {
1205     FPRINTF(fp,"minos: make_specs: Unable to set function precision.\n");
1206     sys->inform=0;
1207     j++;
1208     }
1209     }
1210    
1211     if (sys->rarray[SP1_LSTOL] > D_ZERO) {
1212     sprintf(s56,"%-56s","Linesearch tolerance");
1213     MIOPTR(s56,&(sys->rarray[SP1_LSTOL]),
1214     &out,&(sys->isumm), &(sys->inform));
1215     if ((sys->inform)) {
1216     FPRINTF(fp,"minos: make_specs: Unable to set line search tolerance.\n");
1217     sys->inform=0;
1218     j++;
1219     }
1220     }
1221    
1222     if (sys->rarray[SP1_LUFTO] > D_ZERO) {
1223     sprintf(s56,"%-56s","Lu factor tolerance");
1224     MIOPTR(s56,&(sys->rarray[SP1_LUFTO]),
1225     &out,&(sys->isumm), &(sys->inform));
1226     if ((sys->inform)) {
1227     FPRINTF(fp,"minos: make_specs: Unable to set Lu factor tolerance.\n");
1228     sys->inform=0;
1229     j++;
1230     }
1231     }
1232    
1233     if (sys->rarray[SP1_LUUTO] > D_ZERO) {
1234     sprintf(s56,"%-56s","Lu update tolerance");
1235     MIOPTR(s56,&(sys->rarray[SP1_LUUTO]),
1236     &out,&(sys->isumm), &(sys->inform));
1237     if ((sys->inform)) {
1238     FPRINTF(fp,"minos: make_specs: Unable to set Lu update tolerance.\n");
1239     sys->inform=0;
1240     j++;
1241     }
1242     }
1243    
1244     if (sys->rarray[SP1_LUDTO] > D_ZERO) {
1245     sprintf(s56,"%-56s","Lu density tolerance");
1246     MIOPTR(s56,&(sys->rarray[SP1_LUDTO]),
1247     &out,&(sys->isumm), &(sys->inform));
1248     if ((sys->inform)) {
1249     FPRINTF(fp,"minos: make_specs: Unable to set Lu density tolerance.\n");
1250     sys->inform=0;
1251     j++;
1252     }
1253     }
1254    
1255     if (sys->rarray[SP1_LUSTO] > D_ZERO) {
1256     sprintf(s56,"%-56s","Lu singularity tolerance");
1257     MIOPTR(s56,&(sys->rarray[SP1_LUSTO]),
1258     &out,&(sys->isumm), &(sys->inform));
1259     if ((sys->inform)) {
1260     FPRINTF(fp,
1261     "minos: make_specs: Unable to set Lu singularity tolerance.\n");
1262     sys->inform=0;
1263     j++;
1264     }
1265     }
1266    
1267     if (sys->rarray[SP1_LUWTO] > D_ZERO) {
1268     sprintf(s56,"%-56s","Lu swap tolerance");
1269     MIOPTR(s56,&(sys->rarray[SP1_LUWTO]),
1270     &out,&(sys->isumm), &(sys->inform));
1271     if ((sys->inform)) {
1272     FPRINTF(fp,"minos: make_specs: Unable to set Lu swap tolerance.\n");
1273     sys->inform=0;
1274     j++;
1275     }
1276     }
1277    
1278     if (sys->rarray[SP1_RADIUS] > D_ZERO) {
1279     sprintf(s56,"%-56s","Radius of convergence");
1280     MIOPTR(s56,&(sys->rarray[SP1_RADIUS]),
1281     &out,&(sys->isumm), &(sys->inform));
1282     if ((sys->inform)) {
1283     FPRINTF(fp,"minos: make_specs: Unable to set Radius of convergence.\n");
1284     sys->inform=0;
1285     j++;
1286     }
1287     }
1288    
1289     if (sys->rarray[SP1_SUBSP] > D_ZERO) {
1290     sprintf(s56,"%-56s","Subspace tolerance");
1291     MIOPTR(s56,&(sys->rarray[SP1_SUBSP]),
1292     &out,&(sys->isumm), &(sys->inform));
1293     if ((sys->inform)) {
1294     FPRINTF(fp,"minos: make_specs: Unable to set Subspace tolerance.\n");
1295     sys->inform=0;
1296     j++;
1297     }
1298     }
1299    
1300     if (sys->rarray[SP1_OBJLIM] > D_ZERO) {
1301     sprintf(s56,"%-56s","Unbounded objective value");
1302     MIOPTR(s56,&(sys->rarray[SP1_OBJLIM]),
1303     &out,&(sys->isumm), &(sys->inform));
1304     if ((sys->inform)) {
1305     FPRINTF(fp,
1306     "minos: make_specs: Unable to set Unbounded objective value.\n");
1307     sys->inform=0;
1308     j++;
1309     }
1310     }
1311    
1312     if (sys->rarray[SP1_STEPLM] > D_ZERO) {
1313     sprintf(s56,"%-56s","Unbounded step size");
1314     MIOPTR(s56,&(sys->rarray[SP1_STEPLM]),
1315     &out,&(sys->isumm), &(sys->inform));
1316     if ((sys->inform)) {
1317     FPRINTF(fp,
1318     "minos: make_specs: Unable to set Unbounded step size.\n");
1319     sys->inform=0;
1320     j++;
1321     }
1322     }
1323    
1324     if (sys->rarray[SP1_FDIFF] > D_ZERO) {
1325     sprintf(s56,"%-56s","Difference interval");
1326     MIOPTR(s56,&(sys->rarray[SP1_FDIFF]),
1327     &out,&(sys->isumm), &(sys->inform));
1328     if ((sys->inform)) {
1329     FPRINTF(fp,"minos: make_specs: Unable to set difference interval.\n");
1330     sys->inform=0;
1331     j++;
1332     }
1333     }
1334    
1335     sprintf(s56,"%-56s","Weight on linear objective");
1336     MIOPTR(s56,&(sys->rarray[SP1_FDIFF]),
1337     &out,&(sys->isumm), &(sys->inform));
1338     if ((sys->inform)) {
1339     FPRINTF(fp,
1340     "minos: make_specs: Unable to set Weight on linear objective.\n");
1341     sys->inform=0;
1342     j++;
1343     }
1344    
1345     sprintf(s56,"%-56s","Penalty parameter");
1346     MIOPTR(s56,&(sys->p.rho),
1347     &out,&(sys->isumm), &(sys->inform));
1348     if ((sys->inform)) {
1349     FPRINTF(fp,
1350     "minos: make_specs: Unable to set Penalty parameter(rho)\n");
1351     sys->inform=0;
1352     j++;
1353     }
1354    
1355     if (sys->rarray[SP1_CDIFF] > D_ZERO) {
1356     sprintf(s56,"%-56s","Central difference interval");
1357     MIOPTR(s56,&(sys->rarray[SP1_CDIFF]),
1358     &out,&(sys->isumm), &(sys->inform));
1359     if ((sys->inform)) {
1360     FPRINTF(fp,
1361     "minos: make_specs: Unable to set central difference interval.\n");
1362     sys->inform=0;
1363     j++;
1364     }
1365     }
1366    
1367     sprintf(s72,"%-72s","Minimize");
1368     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1369     if ((sys->inform)) {
1370     FPRINTF(fp,"minos: make_specs: Unable to set objective type\n");
1371     sys->inform=0;
1372     j++;
1373     }
1374    
1375     if (sys->iarray[SP1_SCALE])
1376     sprintf(s72,"%-72s","Scale yes");
1377     else
1378     sprintf(s72,"%-72s","Scale no");
1379     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1380     if ((sys->inform)) {
1381     FPRINTF(fp,"minos: make_specs: Unable to set scaling option\n");
1382     sys->inform=0;
1383     j++;
1384     }
1385    
1386     if (sys->iarray[SP1_USELG])
1387     sprintf(s72,"%-72s","Lagrangian yes");
1388     else
1389     sprintf(s72,"%-72s","Lagrangian no");
1390     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1391     if ((sys->inform)) {
1392     FPRINTF(fp,"minos: make_specs: Unable to set lagrangian option\n");
1393     sys->inform=0;
1394     j++;
1395     }
1396    
1397     sprintf(s72,"%-72s","Debug level 0");
1398     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1399     if ((sys->inform)) {
1400     FPRINTF(fp,"minos: make_specs: Unable to set debug option\n");
1401     sys->inform=0;
1402     j++;
1403     }
1404    
1405     if ( (i=sys->iarray[SP1_MULPR]) ) {
1406     sprintf(s56,"%-56s","Multiple price");
1407     MIOPTI(s56,&i,
1408     &out,&(sys->isumm), &(sys->inform));
1409     if ((sys->inform)) {
1410     FPRINTF(fp,"minos: make_specs: Unable to set multiple price\n");
1411     sys->inform=0;
1412     j++;
1413     }
1414     }
1415    
1416     if ( (i=sys->iarray[SP1_PARPR]) ) {
1417     sprintf(s56,"%-56s","Partial price");
1418     MIOPTI(s56,&i,&out,&(sys->isumm), &(sys->inform));
1419     if ((sys->inform)) {
1420     FPRINTF(fp,"minos: make_specs: Unable to set partial price\n");
1421     sys->inform=0;
1422     j++;
1423     }
1424     }
1425    
1426     sprintf(s56,"%-56s","Print level");
1427     MIOPTI(s56,&(sys->iarray[SP1_JFLXB]),
1428     &out,&(sys->isumm), &(sys->inform));
1429     if ((sys->inform)) {
1430     FPRINTF(fp,"minos: make_specs: Unable to set print level\n");
1431     sys->inform=0;
1432     j++;
1433     }
1434    
1435     if (sys->iarray[SP1_SOLN])
1436     sprintf(s72,"%-72s","Solution yes");
1437     else
1438     sprintf(s72,"%-72s","Solution no");
1439     MIOPT(s72,&out,&(sys->isumm), &(sys->inform));
1440     if ((sys->inform)) {
1441     FPRINTF(fp,"minos: make_specs: Unable to set Solution option\n");
1442     sys->inform=0;
1443     j++;
1444     }
1445    
1446     sprintf(s56,"%-56s","Verify level");
1447     MIOPTI(s56,&(sys->iarray[SP1_VERIFY]),
1448     &out,&(sys->isumm), &(sys->inform));
1449     if ((sys->inform)) {
1450     FPRINTF(fp,"minos: make_specs: Unable to set verify option\n");
1451     sys->inform=0;
1452     j++;
1453     }
1454     if (sys->iarray[SP1_VERIFY]> -1) {
1455     sprintf(s56,"%-56s","Stop constraint check at column");
1456     MIOPTI(s56,&(sys->v.used),&out,&(sys->isumm), &(sys->inform));
1457     if ((sys->inform)) {
1458     FPRINTF(fp,"minos: make_specs: Unable to set verify constraint option\n");
1459     sys->inform=0;
1460     j++;
1461     }
1462     sprintf(s56,"%-56s","Stop objective check at column");
1463     MIOPTI(s56,&(sys->v.used),&out,&(sys->isumm), &(sys->inform));
1464     if ((sys->inform)) {
1465     FPRINTF(fp,"minos: make_specs: Unable to set verify objective option\n");
1466     sys->inform=0;
1467     j++;
1468     }
1469     }
1470    
1471     sprintf(s56,"%-56s","Summary frequency"); /* kaa */
1472     MIOPTI(s56,&(sys->iarray[SP1_LFREQ]),
1473     &out,&(sys->isumm), &(sys->inform));
1474     if ((sys->inform)) {
1475     FPRINTF(fp,"minos: make_specs: Unable to set log freq. option\n");
1476     sys->inform=0;
1477     j++;
1478     }
1479    
1480     sprintf(s56,"%-56s","Check frequency");
1481     MIOPTI(s56,&(sys->iarray[SP1_CFREQ]),
1482     &out,&(sys->isumm), &(sys->inform));
1483     if ((sys->inform)) {
1484     FPRINTF(fp,"minos: make_specs: Unable to set check freq. option\n");
1485     sys->inform=0;
1486     j++;
1487     }
1488    
1489     if (sys->iarray[SP1_EFREQ]>0) {
1490     sprintf(s56,"%-56s","Expand frequency");
1491     MIOPTI(s56,&(sys->iarray[SP1_EFREQ]),
1492     &out,&(sys->isumm), &(sys->inform));
1493     if ((sys->inform)) {
1494     FPRINTF(fp,"minos: make_specs: Unable to set expand frequency\n");
1495     sys->inform=0;
1496     j++;
1497     }
1498     }
1499    
1500     sprintf(s56,"%-56s","Factorization frequency");
1501     MIOPTI(s56,&(sys->iarray[SP1_FFREQ]),
1502     &out,&(sys->isumm), &(sys->inform));
1503     if ((sys->inform)) {
1504     FPRINTF(fp,"minos: make_specs: Unable to set factorization frequency\n");
1505     sys->inform=0;
1506     j++;
1507     }
1508    
1509     ascfree(s56); s56=NULL;
1510     ascfree(s72); s72=NULL;
1511     return (j);
1512     }
1513    
1514     static void insure_bounds(slv1_system_t sys,struct var_variable *var)
1515     /**
1516     *** Insures that the variable value is within ascend bounds.
1517     **/
1518     {
1519     FILE *mif = MIF(sys);
1520     real64 val,low,high;
1521    
1522     if( sys->p.ignore_bounds )
1523     return;
1524    
1525     low = var_lower_bound(var);
1526     high = var_upper_bound(var);
1527     val = var_value(var);
1528     if( low > high ) {
1529     FPRINTF(mif,"Bounds for variable ");
1530     print_var_name(mif,sys,var);
1531     FPRINTF(mif," are inconsistent [%g,%g].\n",low,high);
1532     FPRINTF(mif,"Bounds will be swapped.\n");
1533     var_set_upper_bound(var, low);
1534     var_set_lower_bound(var, high);
1535     low = var_lower_bound(var);
1536     high = var_upper_bound(var);
1537     }
1538    
1539     if( low > val ) {
1540     FPRINTF(mif,"Variable ");
1541     print_var_name(mif,sys,var);
1542     FPRINTF(mif," was found below its lower bound.\n");
1543     FPRINTF(mif,"It will be moved to its lower bound.\n");
1544     var_set_value(var,low);
1545     } else if( val > high ) {
1546     FPRINTF(mif,"Variable ");
1547     print_var_name(mif,sys,var);
1548     FPRINTF(mif," was found above its upper bound.\n");
1549     FPRINTF(mif,"It will be moved to its upper bound.\n");
1550     var_set_value(var, high);
1551     }
1552     }
1553    
1554     static real64 get_linslack(slv1_system_t sys,
1555     int32 row,
1556     struct rel_relation *rel)
1557     /* This returns the slack of a relation as if the relation were a linear
1558     equality. When solving a linear system A.x =b with simplex by rearranging
1559     the system to A.x +s =0, one can either fix RHS=b and drive s->0 or
1560     fix s=-b (by double bounding) and drive RHS->0. Since Minos doesn't
1561     have a _convenient_ way to set RHS, we will use the latter approach.
1562     Simplex lp doesn't notice a difference.
1563     s(x0)=LHS(x0)-RHS(x0)-A.x0 = resid -Ax. Scaling of A,x won't matter.
1564     Does not check sanity of rel,row given. nonlinear rels wil give not
1565     terribly meaningful answers.
1566     Assumes sys->jacobian.mtx, sys->nzlist stuffed.
1567     Ignores any fixed variables that have snuck into the jacobian.
1568     */
1569     {
1570     real64 s=D_ZERO;
1571     int ndx;
1572     struct var_variable *var;
1573    
1574     check_system(sys);
1575     s=relman_eval(rel,&calc_ok);
1576     for( ndx=0 ; ndx < sys->njac ; ndx++ ) {
1577     if( sys->nzlist[ndx].row ==row && sys->nzlist[ndx].col<sys->v.used ) {
1578     var = sys->vlist[mtx_col_to_org(sys->jacobian.mtx,
1579     sys->nzlist[ndx].col)];
1580     s += -var_value(var)* mtx_value(sys->jacobian.mtx ,sys->nzlist+ndx);
1581     }
1582     }
1583     return (s);
1584     }
1585    
1586     static int make_bounds(slv1_system_t sys)
1587     /* Init bl, bu, xn, a. all are scaled by nominals
1588     Note: bounds on slacks must be scaled if we start scaling the
1589     residuals in FUNCON */
1590     {
1591     int32 row,col,ndx,gndx;
1592     struct var_variable **vp;
1593     struct var_variable *var;
1594     var_filter_t vfilter;
1595     struct rel_relation **rp;
1596     rel_filter_t rfilter;
1597     real64 dummy;
1598    
1599     for( vp = sys->vlist ; *vp != NULL ; ++vp )
1600     insure_bounds(sys,*vp);
1601    
1602     /* check variable bounds */
1603     for( col=0 ; col < sys->v.used ; ++col ) {
1604     var = sys->vlist[mtx_col_to_org(sys->jacobian.mtx,col)];
1605    
1606     #if MINOS_DEBUG
1607     FPRINTF(stderr,"setting bounds on ");
1608     print_var_name(MIF(sys),sys,var);
1609     FPRINTF(stderr,"\n");
1610     #endif
1611    
1612     if( (sys->bl[col]=var_lower_bound(var)) <=
1613     -THEIR_INFINITY /* var_nominal(var) */) {
1614     sys->bl[col]=-THEIR_INFINITY;
1615     print_var_name(LIF(sys),sys,var);
1616     FPRINTF(LIF(sys),
1617     " has MINOS lower bound (%g) tighter than ASCEND's.\n",
1618     (-THEIR_INFINITY /* var_nominal(var) */));
1619     }
1620     if( (sys->bu[col]=var_upper_bound(var)) >=
1621     THEIR_INFINITY /* var_nominal(var) */) {
1622     sys->bu[col]=THEIR_INFINITY;
1623     print_var_name(LIF(sys),sys,var);
1624     FPRINTF(LIF(sys)," MINOS upper bound (%g) tighter than ASCEND's.\n",
1625     (THEIR_INFINITY /* var_nominal(var) */));
1626     }
1627     #if MINOS_DEBUG
1628     FPRINTF(stderr,"lo= %g , hi= %g\n", sys->bl[col], sys->bu[col]);
1629     #endif
1630     }
1631    
1632     /* scale variables by nominals */
1633     for( col=0 ; col < sys->v.used ; ++col ) {
1634     var = sys->vlist[mtx_col_to_org(sys->jacobian.mtx,col)];
1635     sys->xn[col]=var_value(var); /* var_nominal(var) */
1636     }
1637    
1638     /* set total scaled jacobian */
1639     calc_ok = TRUE;
1640     rfilter.matchbits = (REL_INCLUDED | REL_ACTIVE);
1641     rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
1642     vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
1643     vfilter.matchvalue = (VAR_INCIDENT | VAR_ACTIVE);
1644     /* vfilter.fixed = var_false;
1645     vfilter.incident = var_true;
1646     vfilter.in_block = var_ignore; */
1647     for( rp=sys->rlist; *rp != NULL; rp++ )
1648     if( rel_apply_filter(*rp,&rfilter) ) {
1649     /* this needs to change */
1650     mtx_mult_row_zero(sys->jacobian.mtx,
1651     mtx_org_to_row(sys->jacobian.mtx,rel_sindex(*rp)),
1652     mtx_ALL_COLS);
1653     /* discard status = */
1654     relman_diffs(*rp,&vfilter,sys->jacobian.mtx,&dummy);
1655     }
1656     if( !calc_ok ) {
1657     FPRINTF(MIF(sys),"!!make_bounds: jacobian calculation error(s).\n");
1658     return (1);
1659     }
1660    
1661     for( gndx=ndx=0 ; ndx < sys->njac ; ++ndx )
1662     if( sys->nzlist[ndx].row < sys->r.used ) {
1663     var = sys->vlist[mtx_col_to_org(sys->jacobian.mtx,
1664     sys->nzlist[ndx].col)];
1665     #if MINOS_DEBUG
1666     /* take this bit out if it works all the time */
1667     if ((gndx+1)>sys->njac){
1668     FPRINTF(MIF(sys),
1669     "make_bounds stuffing a confused gradient(%d)", ndx);
1670     return (1);
1671     }
1672     #endif
1673     sys->a[gndx++] = var_nominal(var) * mtx_value(sys->jacobian.mtx ,
1674     sys->nzlist+ndx);
1675     }
1676    
1677     /* set slack bounds. do after xn, a are set */
1678     for (row=0; row < sys->r.used; row++) {
1679     struct rel_relation *rel = NULL;
1680     rel = sys->rlist[mtx_row_to_org(sys->jacobian.mtx,row)];
1681     ndx=row+sys->v.used; /* move up into slack region of xn, bl, bu */
1682    
1683     switch (rel_relop(rel)) {
1684     case e_equal:
1685     if (row<sys->r.nonlinear)
1686     sys->bl[ndx]=sys->bu[ndx]=D_ZERO;
1687     else
1688     sys->bl[ndx]=sys->bu[ndx]= get_linslack(sys,row,rel);
1689     /* starting residual -A.x0 for linear relations. */
1690     break;
1691     case e_notequal:
1692     sys->bl[ndx]= -THEIR_INFINITY;
1693     sys->bu[ndx]= THEIR_INFINITY;
1694     break;
1695     case e_less:
1696     case e_lesseq:
1697     if (row<sys->r.nonlinear) {
1698     sys->bl[ndx]=D_ZERO;
1699     sys->bu[ndx]=THEIR_INFINITY;
1700     } else {
1701     sys->bl[ndx]= get_linslack(sys,row,rel);
1702     sys->bu[ndx]=THEIR_INFINITY;
1703     }
1704     break;
1705     case e_greater:
1706     case e_greatereq:
1707     if (row<sys->r.nonlinear) {
1708     sys->bu[ndx]=D_ZERO;
1709     sys->bl[ndx]=-THEIR_INFINITY;
1710     } else {
1711     sys->bu[ndx]= get_linslack(sys,row,rel);
1712     sys->bl[ndx]=-THEIR_INFINITY;
1713     }
1714     break;
1715     default:
1716     FPRINTF(stderr,"ERROR: (slv1) make_bounds\n");
1717     FPRINTF(stderr," Undigestible row (%d) passed to MINOS.\n",row);
1718     return (1);
1719     }
1720     } /*for*/
1721     return (0);
1722     }
1723    
1724     static void write_inform(slv1_system_t sys)
1725     /* reports minos exit information */
1726     {
1727     char * res=NULL; /* do not free this ptr ever */
1728     FILE * fp=MIF(sys);
1729     FPRINTF(fp,"\n ***** MINOS inform= %d ***** \n ",sys->inform);
1730     switch (sys->inform) {
1731     case 0 : res="Optimal solution found.";
1732     break;
1733     case 1 : res="The problem is infeasible.";
1734     break;
1735     case 2 : res="The problem is unbounded (or badly scaled).";
1736     break;
1737     case 3 : res="Too many iterations.";
1738     break;
1739     case 4 : res="Stall. The solution has not changed lately.";
1740     break;
1741     case 5 : res="The Superbasics limit is too small.";
1742     break;
1743     case 6 : res="User requested termination in constraint or objective fcn.";
1744     break;
1745     case 7 : res="funobj seems to be giving incorrect gradients.";
1746     break;
1747     case 8 : res="funcon seems to be giving incorrect gradients.";
1748     break;
1749     case 9 : res="The current point cannot be improved.";
1750     break;
1751     case 10 : res="The basis is very ill-conditioned.";
1752     break;
1753     case 11 : res="Cannot find a superbasic to replace a basic variable.";
1754     break;
1755     case 12 : res="Basis factorization requested twice in a row.";
1756     break;
1757     case 13 : res="Near optimal solution found.";
1758     break;
1759     case 20 : res="Not enough storage for basis factorization.";
1760     break;
1761     case 21 : res="Error in basis package.";
1762     break;
1763     case 22 : res=
1764     "The basis is singular after several attempts to factorize it.";
1765     break;
1766     case 30 : /* fall through */
1767     case 31 : res="ASCEND system error.";
1768     break;
1769     case 32 : res="MINOS System error. Wrong number of basic variables.";
1770     break;
1771     case 40 : res="Fatal errors in the MPS file.(What mps file!?!?)";
1772     break;
1773     case 41 : res="Not enough storage to read the MPS file?!?!";
1774     break;
1775     case 42 : res="Not enough storage to solve the problem.";
1776     break;
1777     default : res="Unknown inform returned from MINOSS";
1778     break;
1779     }
1780     FPRINTF(fp,"%s\n",res);
1781     return;
1782     }
1783    
1784     /*
1785     * Sets up the interface to MINOS and then invokes it.
1786     */
1787     static int invoke_minos(slv1_system_t sys)
1788     {
1789     char *start; /* never free this ptr */
1790     static char names[48]=" ";
1791     static int nname=1,iobj=0, name1=0,name2=0,mincor,ns,ninf;
1792     static double objadd=0.0,obj, sinf;
1793     struct var_variable **vp=NULL;
1794     check_system(sys);
1795    
1796     if (!sys->s.ready_to_solve) /* fail if not presolved */
1797     return (1);
1798     /*
1799     * Note: due to the bounding process by which linear relations get their
1800     * if (sys->basis) start="Cold ";
1801     * else start="Warm ";
1802     * proper RHS, all starts are cold for now.
1803     */
1804     start="Cold ";
1805     if (make_specs(sys)>0) {
1806     FPRINTF(MIF(sys),"Unable to invoke MINOS; bad return from make_specs.\n");
1807     sys->s.ready_to_solve=FALSE;
1808     return (1);
1809     }
1810     if (make_bounds(sys)>0) {
1811     FPRINTF(MIF(sys),"Unable to invoke MINOS; bad return from make_bounds.\n");
1812     sys->s.ready_to_solve=FALSE;
1813     return (1);
1814     }
1815    
1816     g_sys = sys;
1817    
1818     FPRINTF(LIF(sys),"About to call minos.\n");
1819     print_vector(sys,sys->xn,sys->v.used,"xn");
1820     print_vector(sys,sys->a,sys->njac,"a");
1821     print_vector(sys,sys->bl,sys->nb,"bl");
1822     print_vector(sys,sys->bu,sys->nb,"bu");
1823     /*
1824     * While !mem ok get more mem. final get will be ok and solve.
1825     */
1826     do {
1827     sys->inform=0;
1828     iteration_begins(sys);
1829     MINOSS(start,
1830     &(sys->mrows), /* m */
1831     &(sys->v.used), /* n */
1832     &(sys->nb),
1833     &(sys->njac), /* ne */
1834     &nname,
1835     &(sys->r.nonlinear), /* nncon */
1836     &(sys->nnobj),
1837     &(sys->v.nonlinear), /* nnjac */
1838     &iobj,
1839     &objadd,
1840     &names[0],
1841     sys->a,
1842     sys->ha,
1843     sys->ka,
1844     sys->bl,
1845     sys->bu,
1846     &name1, /* create name1, name2 if nname !=1 */
1847     &name2,
1848     sys->hs,
1849     sys->xn,
1850     sys->pi,
1851     sys->rc,
1852     &(sys->inform),
1853     &mincor,
1854     &ns,
1855     &ninf,
1856     &sinf,
1857     &obj,
1858     sys->z,
1859     &(sys->nwcore)); /* size of workspace in doubles */
1860     iteration_ends(sys);
1861     if (sys->inform==42) {
1862     free_vector(sys->z);
1863     sys->nwcore=2*mincor;
1864     sys->z=alloc_vector(sys->nwcore);
1865     }
1866     if (sys->inform==20) {
1867     free_vector(sys->z);
1868     sys->nwcore=2*sys->nwcore;
1869     sys->z=alloc_vector(sys->nwcore);
1870     }
1871     FPRINTF(LIF(sys), "MINOS workspace %d bytes.\n",
1872     (sys->nwcore * (int) sizeof(real64) ) );
1873     } while (sys->inform ==20 || sys->inform ==42);
1874    
1875     FPRINTF(LIF(sys),"After calling minos.\n");
1876     print_vector(sys,sys->xn,sys->v.used,"xn");
1877     print_vector(sys,sys->a,sys->njac,"a");
1878     print_vector(sys,sys->bl,sys->nb,"bl");
1879     print_vector(sys,sys->bu,sys->nb,"bu");
1880     print_vector(sys,sys->pi,sys->mrows,"pi");
1881     install_allvars(sys,sys->xn);
1882     sys->basis=TRUE;
1883     /*
1884     * note sys->inform==6 => user interface halt in
1885     * funobj/funcon or eval panic
1886     */
1887     write_inform(sys);
1888     GETCOMMON(&(sys->s.iteration),&(sys->s.block.iteration));
1889     for( vp = sys->vlist ; *vp != NULL ; ++vp )
1890     insure_bounds(sys,*vp);
1891     iteration_begins(sys);
1892     calc_residuals(sys);
1893     if (sys->iarray[SP1_LFREQ]) /* added kaa */
1894     print_output(stdout,sys);
1895     iteration_ends(sys);
1896     return 0;
1897     }
1898    
1899     /*********************************************************************
1900     *** External routines ***
1901     See slv.h
1902     *********************************************************************/
1903    
1904     static struct rel_relation **internal_rlist(struct rel_relation * *rlist)
1905     /* Converts rlist to its non-NULL equivalent. */
1906     {
1907     static struct rel_relation *empty_list[] = {NULL};
1908     return( rlist==NULL ? empty_list : rlist );
1909     }
1910    
1911     static SlvClientToken slv1_create(slv_system_t server, int *statusindex)
1912     {
1913     slv1_system_t sys;
1914    
1915     sys = (slv1_system_t)ascmalloc( sizeof(struct slv1_system_structure) );
1916     mem_zero_byte_cast( sys , 0 , sizeof(struct slv1_system_structure) );
1917     sys->integrity = OK;
1918    
1919     /* Only need to initialize non-zeros */
1920     sys->p.output.more_important = stdout;
1921     sys->p.time_limit = (double)30.0;
1922     sys->p.iteration_limit = 20;
1923     sys->p.tolerance.singular = (double)1e-8;
1924     sys->p.tolerance.feasible = (double)1e-6;
1925     sys->p.tolerance.stationary = (double)1e-8;
1926     sys->p.partition = TRUE;
1927     sys->p.whose = slv1_solver_number;
1928     sys->p.rho = (double)0.1;
1929     sys->s.ok = TRUE;
1930     sys->s.calc_ok = TRUE;
1931     sys->rlist = internal_rlist(NULL);
1932     sys->nwcore=10000;
1933     sys->p.sp.iap=&(sys->iarray[0]);
1934     sys->p.sp.rap=&(sys->rarray[0]);
1935     sys->iarray[SP1_COMPLETION]=0;
1936     sys->iarray[SP1_CRASH]=1;
1937     sys->iarray[SP1_MINITS]=40;
1938     sys->iarray[SP1_DERIV]=3;
1939     sys->iarray[SP1_MULPR]=1;
1940     sys->iarray[SP1_USELG]=1;
1941     sys->iarray[SP1_VERIFY]=-1;
1942     sys->iarray[SP1_CFREQ]=30;
1943     sys->iarray[SP1_LFREQ]=10;
1944     sys->iarray[SP1_FFREQ]=50;
1945     sys->iarray[SP1_EFREQ]=10000;
1946     sys->iarray[SP1_LCONS]=1; /* added by kaa */
1947     sys->rarray[SP1_FPREC]=(double)1.0e-10;
1948     sys->rarray[SP1_LUFTO]=(double)10;
1949     sys->rarray[SP1_LUUTO]=(double)10;
1950     sys->rarray[SP1_LSTOL]=(double)0.1;
1951     sys->rarray[SP1_RADIUS]=(double)0.01;
1952     sys->rarray[SP1_SUBSP]=(double)0.5;
1953     sys->rarray[SP1_OBJLIM]=THEIR_INFINITY;
1954     sys->rarray[SP1_STEPLM]=(double)1.0e10;
1955     return(sys);
1956     }
1957    
1958     static void slv1_set_rel_list( slv1_system_t sys, struct rel_relation **rlist)
1959     {
1960     static struct rel_relation *empty_list[] = {NULL};
1961    
1962     check_system(sys);
1963     sys->rlist_user = rlist;
1964     sys->rlist = (rlist == NULL ? empty_list : rlist);
1965     sys->s.ready_to_solve = FALSE;
1966     }
1967    
1968     static struct rel_relation **slv1_get_rel_list( slv1_system_t sys)
1969     {
1970     check_system(sys);
1971     return( sys->rlist_user );
1972     }
1973    
1974     static void slv1_set_var_list(slv1_system_t sys,struct var_variable **vlist)
1975     {
1976     static struct var_variable *empty_list[] = {NULL};
1977     check_system(sys);
1978     if( sys->vlist_user == NULL )
1979     if( sys->vlist != NULL && pl_length(sys->vlist) > 0 )
1980     ascfree( (POINTER)(sys->vlist) );
1981     sys->vlist_user = vlist;
1982     sys->vlist = (vlist==NULL ? empty_list : vlist);
1983     sys->s.ready_to_solve = FALSE;
1984     }
1985    
1986     static struct var_variable **slv1_get_var_list(slv1_system_t sys)
1987     {
1988     check_system(sys);
1989     return( sys->vlist_user );
1990     }
1991    
1992     static int slv1_eligible_solver(slv_system_t server)
1993     {
1994     struct rel_relation **rp;
1995     for( rp=slv_get_solvers_rel_list(server); *rp != NULL ; ++rp ) {
1996     if( rel_less(*rp) || rel_greater(*rp) ) return(FALSE);
1997     }
1998     return(TRUE);
1999     }
2000    
2001     static void slv1_get_parameters(slv_system_t server, SlvClientToken asys,
2002     slv_parameters_t *parameters)
2003     {
2004     slv1_system_t sys;
2005     sys = SLV1(asys);
2006     if (check_system(sys)) return;
2007     mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
2008     }
2009    
2010     static void slv1_set_parameters(slv_system_t server, SlvClientToken asys,
2011     slv_parameters_t *parameters)
2012     {
2013     slv1_system_t sys;
2014     sys = SLV1(asys);
2015     if (check_system(sys)) return;
2016     mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
2017     }
2018    
2019     static void slv1_get_status(slv_system_t server, SlvClientToken asys,
2020     slv_status_t *status)
2021     {
2022     slv1_system_t sys;
2023     sys = SLV1(asys);
2024     if (check_system(sys)) return;
2025     mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
2026     }
2027    
2028     static linsolqr_system_t slv1_get_linsolqr_sys(slv_system_t server,
2029     SlvClientToken asys)
2030     {
2031     return( NULL );
2032     }
2033    
2034     static linsol_system_t slv1_get_linsol_sys(slv_system_t server,
2035     SlvClientToken asys)
2036     {
2037     slv1_system_t sys;
2038     sys = SLV1(asys);
2039     if (check_system(sys)) return NULL;
2040     return sys->jacobian.sys;
2041     }
2042    
2043     void slv1_dump_internals(slv_system_t server,
2044     SlvClientToken asys,int level)
2045     {
2046     int i;
2047     slv1_system_t sys;
2048     sys = SLV1(asys);
2049     check_system(sys);
2050     /* for reference:
2051     struct jacobian_data jacobian;
2052     real64
2053     *z, minos workspace, includes hessian, etc
2054     *rc, reduced costs as known to minos
2055     *pi, lagrange multipliers as known to minos for nl constraints
2056     *bl, lower bound values (scaled) (nb)
2057     *bu, upper bound values (scaled)(nb)
2058     int32
2059     *hs; variable status (type) vector(nb)
2060     */
2061     if (level==5) { /*paramvec dump*/
2062     for (i=0;i<slv1_IA_SIZE;i++)
2063     FPRINTF(stderr,"iarray[%d]= %d\n",
2064     i,sys->iarray[i]);
2065     for (i=0;i<slv1_RA_SIZE;i++)
2066     FPRINTF(stderr,"rarray[%d]= %.16g\n",
2067     i,sys->rarray[i]);
2068     }
2069    
2070     if (level==4) /*nzlist dump*/
2071     for (i=0;i<sys->njac;i++)
2072     FPRINTF(stderr,"nzl[%d] r%d c%d\n",
2073     i,sys->nzlist[i].row,sys->nzlist[i].col);
2074    
2075     if (level==2) { /* dump jacobian */
2076     FPRINTF(stderr,"jacobian:\n");
2077     for (i=0; i<sys->njac; i++)
2078     FPRINTF(stderr,"nzlist(%d)[%d(ha=%d),col=%d] = %g\n",
2079     i, sys->nzlist[i].row,sys->ha[i],sys->nzlist[i].col, sys->a[i]);
2080     }
2081     if (level==3) { /* dump x, col indexes */
2082     FPRINTF(stderr,"xn vector:\n");
2083     for (i=0; i<sys->nb; i++)
2084     FPRINTF(stderr,"jcol(%d)-ka>%d %g < %g < %g\n",
2085     i,(i<sys->v.used)?sys->ka[i]:-1, sys->bl[i],sys->xn[i],sys->bu[i]);
2086     }
2087     if (level==1) {
2088     FPRINTF(stderr,"file channels:\n");
2089     FPRINTF(stderr,"ispecs = %d\n",sys->ispecs);
2090     FPRINTF(stderr,"iprint = %d\n",sys->iprint);
2091     FPRINTF(stderr,"isumm = %d\n",sys->isumm);
2092     write_inform(sys);
2093     FPRINTF(stderr," maxndx= %d\n",sys->maxndx);
2094     FPRINTF(stderr,"r.used (minos m) = %d (%d)\n",sys->r.used,sys->mrows);
2095     FPRINTF(stderr,"r.nonlinear (minos nncon )= %d\n",sys->r.nonlinear);
2096     FPRINTF(stderr,"v.used (minoss n) = %d\n",sys->v.used);
2097     FPRINTF(stderr,"v.nonlinear (minos nnjac) = %d\n",sys->v.nonlinear);
2098     FPRINTF(stderr,"nb= var+con = %d\n",sys->nb);
2099     FPRINTF(stderr,"nwcore = %d\n",sys->nwcore);
2100     FPRINTF(stderr,"njac (minos ne) = %d\n",sys->njac);
2101     FPRINTF(stderr,"nnobj = %d\n",sys->nnobj);
2102     }
2103     if (level>10)
2104     for (i=1;i<10;i++) slv1_dump_internals(server,asys,i);
2105     }
2106    
2107     /* sort a range (all with the same column index) in nzlist
2108     so that it is increasing row order. predicated on sparsity.
2109     */
2110    
2111     static int nzsort(mtx_coord_t *l, int32 begin, int32 end)
2112     {
2113     int32 t,curi,curj;
2114     if (end<begin) {
2115     #if MINOS_DEBUG
2116     FPRINTF(stderr, "Garbage sent to nzsort: begin %d, end %d.",begin,end);
2117     #endif
2118     return 1;
2119     /* this implies that a var appears in the objective which does not
2120     appear in the constraints. */
2121     }
2122     t=end-begin;
2123     if (!t) return 0; /* 1 element -> do nothing */
2124     if (t==1) { /* if pair, be quick */
2125     if (l[begin].row > l[end].row ) { /* wrong ordered pair */
2126     t=l[end].row;
2127     l[end].row=l[begin].row;
2128     l[begin].row=t;
2129     }
2130     return 0;
2131     }
2132     for (curi=begin; curi<end; curi++) /* yes, it's the evil bubble sort */
2133     for (curj=curi+1; curj<=end; curj++)
2134     if ( l[curi].row > l[curj].row ) {
2135     t=l[curj].row;
2136     l[curj].row=l[curi].row;
2137     l[curi].row=t;
2138     }
2139     return 0;
2140     }
2141     /****
2142     * Presolve creates a sparse jacobian matrix with rows/columns reordered:
2143     * (vars in nonlinear rels)(obj vars not in rels if o nl)(linear vars)(unincid)
2144     * and rows in order (nonlinear rels) (linear rels) (unincluded rels)
2145     *
2146     * The jacobian structure is recorded and manipulated via nzlist, an array
2147     * of matrix coordinates. The coordinates refer to positions in the reordered
2148     * matrix, not the original variable list positions. use row/col_to_org to find
2149     * out what positions in the matrix go with what equations/variables.
2150     * The nzlist structure is used to stuff minos a, the total jacobian, and
2151     * minos g the nonlinear constraint derivatives in FUNCON. The array a
2152     * is indexed by arrays ha, ka which must be correctly set for minos sparse
2153     * jacobian option to work right. Note that the indexes in ha, ka count from
2154     * 1 instead of 0.
2155     *
2156     * Elements of ha must be in ascending minos row order within a column.
2157     * hence, elements of nzlist must also be in increasing mtx column,row
2158     * order.
2159     *
2160     * notes: V.used is the number of vars that would pass filter_var_solver, and
2161     * hence the coord of the first fixed var, if there are any.
2162     * V.nonlinear is the number of vars in nonlinear rels/objective, and
2163     * hence the coord of the first linear var, if there are any.
2164     *
2165     *
2166     ****/
2167     /* THIS USED TO RETURN AN INT ( = 1 for fail ) */
2168     void slv1_presolve(slv_system_t server, SlvClientToken asys)
2169     {
2170     slv1_system_t sys = SLV1(asys);
2171     struct var_variable **vp;
2172     var_filter_t vfilter;
2173     struct rel_relation **rp;
2174     rel_filter_t rfilter;
2175     mtx_range_t nl_rows;
2176     mtx_coord_t nz;
2177     real64 value;
2178     int32 col,el,ello,elhi,cap;
2179    
2180     check_system(sys);
2181     if( sys->vlist == NULL ) {
2182     FPRINTF(stderr,"ERROR: (slv1) slv1_presolve\n");
2183     FPRINTF(stderr," Variable list was never set.\n");
2184     return;
2185     }
2186     if( sys->rlist == NULL ) {
2187     FPRINTF(stderr,"ERROR: (slv1) slv1_presolve\n");
2188     FPRINTF(stderr," Relation list was never set.\n");
2189     return;
2190     }
2191    
2192     /* iteration_begins() does no damage and is needed so */
2193     /* iteration_ends() can be called afterwards. */
2194     iteration_begins(sys);
2195    
2196     /* Reset global iteration count and cpu-elapsed. */
2197     sys->s.iteration = 0;
2198     sys->s.cpu_elapsed = D_ZERO;
2199     sys->panic=FALSE;
2200    
2201     /* Determine system's extent */
2202     if( sys->vlist_user == NULL ) {
2203     Asc_Panic(2, "slv1_presolve",
2204     "the logic in slve presolve needs to be modified\n");
2205     }
2206    
2207     sys->maxndx = 0;
2208     for( vp=sys->vlist,cap=0 ; *vp != NULL ; ++vp ) {
2209     var_set_sindex(*vp,cap++);
2210     var_set_in_block(*vp,FALSE);
2211     }
2212     sys->maxndx = cap;
2213     for( rp=sys->rlist,cap=0 ; *rp != NULL ; ++rp ) {
2214     rel_set_sindex(*rp,cap++);
2215     rel_set_in_block(*rp,FALSE);
2216     rel_set_satisfied(*rp,FALSE);
2217     }
2218     sys->maxndx = MAX(sys->maxndx,cap);
2219    
2220     /* Set up jacobian system. sys,mtx,rhs should be created/destroyed together */
2221     make_jacobian(sys);
2222    
2223     /* This is not the solver's job
2224     if( sys->obj != NULL )
2225     exprman_decide_incidence_obj(sys->obj);
2226    
2227     rfilter.matchbits = (REL_INCLUDED | REL_ACTIVE);
2228     rfilter.matchvalue = (REL_INCLUDED| REL_ACTIVE );
2229     for( rp=sys->rlist; *rp != NULL; ++rp )
2230     if( rel_apply_filter(*rp,&rfilter) )
2231     relman_decide_incidence(*rp);
2232     */
2233    
2234     vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_ACTIVE);
2235     vfilter.matchvalue = (VAR_INCIDENT | VAR _ACTIVE );
2236     /*
2237     vfilter.fixed = var_false;
2238     vfilter.incident = var_true;
2239     vfilter.in_block = var_ignore; */
2240     mtx_clear(sys->jacobian.mtx);
2241     mtx_clear(sys->jacobian.objmtx);
2242     sys->r.used = sys->r.nonlinear = 0;
2243     for( rp=sys->rlist; *rp != NULL; ++rp ) { /* if constraints, nl first */
2244     rel_set_in_block(*rp,FALSE);
2245     if( rel_apply_filter(*rp,&rfilter) ) {
2246     rel_set_in_block(*rp,TRUE);
2247     relman_get_incidence(*rp,&vfilter,sys->jacobian.mtx);
2248     mtx_swap_rows(sys->jacobian.mtx,rel_sindex(*rp),sys->r.used);
2249     if( !relman_is_linear(*rp,&vfilter) || sys->iarray[SP1_LCONS] )
2250     mtx_swap_rows(sys->jacobian.mtx,sys->r.used,sys->r.nonlinear++);
2251     ++(sys->r.used);
2252     }
2253     }
2254    
2255     /* put all the nonlinear rows at the top */
2256     sys->mrows=(sys->r.used==0)? 1:sys->r.used;
2257    
2258     /* This bit assumes all free vars in nonlinear relations occur nonlinearly. */
2259     /* We need badly a var_linear_in_rel/expression boolean */
2260     /* move nl_vars left in constraints and obj mtx */
2261     nl_rows.low = 0;
2262     nl_rows.high = sys->r.nonlinear - 1;
2263     sys->v.used = sys->v.nonlinear = 0;
2264     for( vp=sys->vlist ; *vp != NULL ; vp++ ) {
2265     var_set_in_block(*vp,FALSE);
2266     if( var_apply_filter(*vp,&vfilter) ) {
2267     make_nominal_positive(sys,*vp);
2268     var_set_in_block(*vp,TRUE);
2269     /* since solver variable (else wouldn't be here), move var left to edge of
2270     current solver variable region of jacobian and objective.
2271     */
2272     mtx_swap_cols(sys->jacobian.mtx,var_sindex(*vp),(nz.col=sys->v.used));
2273     mtx_swap_cols(sys->jacobian.objmtx,var_sindex(*vp),(nz.col=sys->v.used));
2274     nz.row = mtx_FIRST;
2275     /* if this column (var) has an incidence in row range 0, r.nonlinear - 1, move
2276     it to the left and move up nonlinear marker.
2277     */
2278     value = mtx_next_in_col(sys->jacobian.mtx,&nz,&nl_rows);
2279     if( nz.row != mtx_LAST ) {
2280     mtx_swap_cols(sys->jacobian.mtx,sys->v.used,sys->v.nonlinear);
2281     mtx_swap_cols(sys->jacobian.objmtx,sys->v.used,sys->v.nonlinear++);
2282     }
2283     /* move up solver variable marker
2284     */
2285     ++(sys->v.used);
2286     }
2287     }
2288     /* go on to the next variable in if */
2289     /* end of for */
2290    
2291     /* this bit assumes all free vars in nonlinear objective occur nonlinearly. */
2292     sys->nnobj=0; /* no nonlinear variables in obj if obj linear or null*/
2293     if( sys->obj != NULL ) {
2294     if (!relman_is_linear(sys->obj,&vfilter) || sys->iarray[SP1_LCONS]) {
2295     /* build a list of incident in obj */
2296     struct var_variable **objvars;
2297     int32 col;
2298     objvars = rel_incidence_list(sys->obj);
2299     for( vp=objvars ; *vp != NULL ; vp++ )
2300     if(!var_fixed(*vp) && var_active(*vp) &&
2301     ( col=mtx_org_to_col(sys->jacobian.objmtx,var_sindex(*vp)) ) >=
2302     sys->v.nonlinear ) {
2303     mtx_swap_cols(sys->jacobian.mtx,col,sys->v.nonlinear);
2304     mtx_swap_cols(sys->jacobian.objmtx,col,sys->v.nonlinear++);
2305     ++(sys->nnobj);
2306     }
2307    
2308     FPRINTF(LIF(sys),"Presolve: nnobj= %d, v.nonlinear= %d\n",
2309     sys->nnobj,sys->v.nonlinear);
2310     sys->nnobj=sys->v.nonlinear=MAX(sys->nnobj,sys->v.nonlinear);
2311     ascfree( (POINTER)objvars );
2312     }
2313     }
2314    
2315     /*************************************************************************
2316     From here on in the MINOS invocation algorithm, don't call clear
2317     on mtx_matrix_t's or you will screw up the nextinrow/col operators
2318     used for stuffing gradients. mtx_mult_row if 0s needed.
2319     *************************************************************************/
2320    
2321     /* count nonzeros */
2322     sys->njac = 0;
2323     for( nz.row=0 ; nz.row < sys->maxndx ; ++nz.row )
2324     sys->njac += mtx_nonzeros_in_row(sys->jacobian.mtx,nz.row,mtx_ALL_COLS);
2325     /* allocate nzlist with at least one element*/
2326     free_unless_null( (POINTER)(sys->nzlist) );
2327     sys->nzlist = (mtx_coord_t *)ascmalloc((sys->njac+1) * sizeof(mtx_coord_t));
2328     sys->nzlist[sys->njac].row=sys->nzlist[sys->njac].col=sys->r.used+5000000;
2329    
2330     /* copy coord to nzlist for use in calc routines */
2331     for( sys->njac=nz.col=0 ; nz.col < sys->maxndx ; ++nz.col ) {
2332     nz.row = mtx_FIRST;
2333     while( value = mtx_next_in_col(sys->jacobian.mtx,&nz,mtx_ALL_ROWS),
2334     nz.row != mtx_LAST )
2335     if( value != D_ZERO )
2336     mem_copy_cast(&nz , &(sys->nzlist[sys->njac++]) ,
2337     sizeof(mtx_coord_t) );
2338     }
2339     /* ascmalloc ha, ka */
2340     free_vector(sys->ha); /* create row index array, and init to 1s */
2341     sys->ha=(int *)ascmalloc((sys->njac+1)*sizeof(int));
2342     for (nz.row=0;nz.row<=sys->njac;nz.row++)
2343     sys->ha[nz.row]=1; /* was nz.col */
2344    
2345     free_vector(sys->ka); /* create col offset array */
2346     sys->ka=(int *)asccalloc((sys->v.used+1),sizeof(int));/*this is not padding*/
2347    
2348     /* secondary rowsort nzlist. already column sorted by creation process */
2349     if (sys->njac > 1) {
2350     ello=elhi=el=0;
2351     for (col=0; col <sys->v.used; col++) {
2352     /* find next colstart in nzlist */
2353     for (el=ello;sys->nzlist[el].col<=col;el++);
2354     elhi=el-1; /* set this col end*/
2355     if (!nzsort(sys->nzlist,ello,elhi) ) ello=el;
2356     /* rowsort nzlist elements from this col */
2357     }
2358     }
2359    
2360     /* translate nzlist into ha, ka */
2361     /* arys number from 1 in fortran, as we have initted ha */
2362     for (nz.row=0; nz.row<sys->njac;nz.row++)
2363     sys->ha[nz.row]+=(sys->nzlist[nz.row].row);
2364    
2365     for (nz.row=nz.col=0; nz.row < sys->njac; nz.row++)
2366     if (sys->nzlist[nz.row].col > nz.col) sys->ka[++nz.col]=nz.row;
2367     /* bump up the ka to fortran */
2368     for (nz.row=0; nz.row<sys->v.used; ++(sys->ka[nz.row++]));
2369     sys->ka[sys->v.used]=sys->njac+1;
2370    
2371    
2372     /* never do this, or it may upset the consistency of next_in_row.
2373     mtx_clear_region(sys->jacobian.mtx,mtx_ENTIRE_MATRIX);
2374     */
2375     sys->basis =
2376     sys->s.over_defined =
2377     sys->s.under_defined =
2378     sys->s.struct_singular =
2379     sys->s.converged =
2380     sys->s.diverged =
2381     sys->s.inconsistent = FALSE;
2382     sys->s.block.number_of = 1;
2383     sys->nb=sys->mrows+sys->v.used;
2384     sys->nwcore=110*(pl_length(sys->rlist)+1);
2385     /* minos guide is 100*r.used. we tend to nonlinearity */
2386    
2387     sys->s.block.previous_total_size = 0;
2388     sys->s.block.current_block = 0;
2389     sys->s.block.current_size = sys->maxndx;
2390     sys->s.block.iteration = 0;
2391     sys->s.block.cpu_elapsed = D_ZERO;
2392     sys->s.calc_ok = TRUE;
2393    
2394     free_vector(sys->z); /* create minos workspace */
2395     sys->z=alloc_zero_vector(sys->nwcore);
2396     if (!sys->z) {
2397     FPRINTF(stderr,"Error mallocing minos workspace array!");
2398     return;
2399     }
2400    
2401     free_vector(sys->bl); /* create lowerbound array */
2402     sys->bl=alloc_zero_vector(sys->nb);
2403     if (!sys->bl) {
2404     FPRINTF(stderr,"Error mallocing lower bound array!");
2405     return;
2406     }
2407    
2408     free_vector(sys->bu); /* create upperbound array */
2409     sys->bu=alloc_zero_vector(sys->nb);
2410     if (!sys->bu) {
2411     FPRINTF(stderr,"Error mallocing upper bound array!");
2412     return;
2413     }
2414    
2415     free_vector(sys->hs); /* create var state flag array */
2416     sys->hs=(int *)asccalloc(sys->nb,sizeof(int));
2417     if (!sys->hs) {
2418     FPRINTF(stderr,"Error mallocing variable state array!");
2419     return;
2420     }
2421    
2422     free_vector(sys->pi); /* create lagrange multiplier array */
2423     sys->pi=alloc_zero_vector(sys->mrows);
2424     if (!sys->pi) {
2425     FPRINTF(stderr,"Error mallocing multiplier array!");
2426     return;
2427     }
2428    
2429     free_vector(sys->rc); /* create reduced cost array */
2430     sys->rc=alloc_zero_vector(sys->nb);
2431     if (!sys->rc) {
2432     FPRINTF(stderr,"Error mallocing reduced cost array!"