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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 83089 byte(s)
Setting up web subdirectory in repository
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!");
2433 return;
2434 }
2435
2436 free_vector(sys->xn); /* create solution array */
2437 sys->xn=alloc_zero_vector(sys->nb);
2438 if (!sys->xn) {
2439 FPRINTF(stderr,"Error mallocing variable array!");
2440 return;
2441 }
2442
2443 free_vector(sys->a); /* create jacobian element array */
2444 sys->a=alloc_zero_vector(sys->njac+1);
2445 if (!sys->a) {
2446 FPRINTF(stderr,"Error mallocing jacobian array!");
2447 return;
2448 }
2449
2450 /* set initial i/o unit options and nwcore.*/
2451 if (sys->iarray[SP1_SUMMY])
2452 sys->iprint=6;
2453 else
2454 sys->iprint=0;
2455 if (sys->iarray[SP1_FSUMY])
2456 sys->isumm=9;
2457 else
2458 sys->isumm=0;
2459 sys->inform=0;
2460 sys->ispecs=0;
2461 MISPEC(&(sys->ispecs),
2462 &(sys->iprint),
2463 &(sys->isumm),
2464 &(sys->nwcore),
2465 &(sys->inform));
2466
2467 iteration_ends(sys);
2468 return;
2469 }
2470
2471 boolean slv1_change_basis(slv1_system_t sys,int32 var, mtx_range_t *rng){
2472 return FALSE;
2473 }
2474 /* doesn't work right now */
2475 static void slv1_resolve(slv_system_t server, SlvClientToken asys)
2476 {
2477 struct rel_relation **rp;
2478 slv1_system_t sys;
2479 sys = SLV1(asys);
2480 check_system(sys);
2481
2482 /* iteration_begins() does no damage and is needed so */
2483 /* iteration_ends() can be called afterwards. */
2484 iteration_begins(sys);
2485
2486 /* Reset global iteration count and cpu-elapsed. */
2487 sys->s.iteration = 0;
2488 sys->s.cpu_elapsed = D_ZERO;
2489
2490 for( rp=sys->rlist ; *rp != NULL ; ++rp )
2491 if( rel_included(*rp) && rel_active(*rp))
2492 rel_set_satisfied(*rp,FALSE);
2493
2494 mtx_clear_region(sys->jacobian.mtx,mtx_ENTIRE_MATRIX);
2495 sys->s.over_defined =
2496 sys->s.under_defined =
2497 sys->s.struct_singular =
2498 sys->s.converged =
2499 sys->s.diverged =
2500 sys->s.inconsistent = FALSE;
2501 sys->s.block.number_of = 1;
2502
2503 sys->s.block.previous_total_size = 0;
2504 sys->s.block.current_block = 0;
2505 sys->s.block.current_size = sys->maxndx;
2506 sys->s.block.iteration = 0;
2507 sys->s.block.cpu_elapsed = D_ZERO;
2508 sys->s.calc_ok = TRUE;
2509
2510 sys->basis = FALSE;
2511
2512 iteration_ends(sys);
2513 }
2514
2515 /* the following 2 functions need to be combined (or something) */
2516 static void slv1_solve(slv_system_t server, SlvClientToken asys)
2517 {
2518 slv1_system_t sys;
2519 sys = SLV1(asys);
2520 if (server == NULL || sys==NULL) return;
2521 if (check_system(SLV1(sys))) return;
2522 if( !sys->s.ready_to_solve ) {
2523 FPRINTF(stderr,"ERROR: (slv1) slv1_solve\n");
2524 FPRINTF(stderr," System is not ready to solve.\n");
2525 return;
2526 }
2527 invoke_minos(sys);
2528 }
2529
2530 static void slv1_iterate(slv_system_t server, SlvClientToken asys)
2531 {
2532 slv1_system_t sys;
2533 sys = SLV1(asys);
2534 if (server == NULL || sys==NULL) return;
2535 if (check_system(SLV1(sys))) return;
2536 slv1_solve(server,asys);
2537 }
2538
2539 static mtx_matrix_t slv1_get_jacobian(slv_system_t server, SlvClientToken sys)
2540 {
2541 if (server == NULL || sys==NULL) return NULL;
2542 if (check_system(SLV1(sys))) return NULL;
2543 return SLV1(sys)->jacobian.sys;
2544 }
2545
2546 static int slv1_destroy(slv_system_t server, SlvClientToken asys)
2547 {
2548 slv1_system_t sys;
2549 sys = SLV1(asys);
2550 if (check_system(sys)) return 1;
2551 slv1_set_var_list(sys,NULL); /* also destroys list if internal */
2552 slv1_set_rel_list(sys,NULL);
2553 /* the above lists and obj will be destroyed by destroy_from_instance */
2554
2555 destroy_jacobian(sys);
2556 free_unless_null( (POINTER)(sys->nzlist) );
2557 sys->nzlist=NULL;
2558
2559 /* free int lists */
2560 free_unless_null( (POINTER)(sys->ka) );
2561 free_unless_null( (POINTER)(sys->hs) );
2562 free_unless_null( (POINTER)(sys->ha) );
2563 sys->ha=NULL;
2564 sys->hs=NULL;
2565 sys->ka=NULL;
2566
2567 /* free double lists */
2568 free_vector(sys->z);
2569 free_vector(sys->rc);
2570 free_vector(sys->a);
2571 free_vector(sys->bu);
2572 free_vector(sys->bl);
2573 free_vector(sys->xn);
2574 free_vector(sys->pi);
2575 sys->pi=NULL;
2576 sys->xn=NULL;
2577 sys->bl=NULL;
2578 sys->bu=NULL;
2579 sys->a=NULL;
2580 sys->rc=NULL;
2581 sys->z=NULL;
2582
2583 sys->integrity = DESTROYED;
2584 ascfree( (POINTER)sys );
2585 return 0;
2586 }
2587
2588 int slv1_register(SlvFunctionsT *sft)
2589 {
2590 if (sft==NULL) {
2591 FPRINTF(stderr,"slv1_register called with NULL pointer\n");
2592 return 1;
2593 }
2594
2595 sft->name = "MINOS";
2596 sft->ccreate = slv1_create;
2597 sft->cdestroy = slv1_destroy;
2598 sft->celigible = slv1_eligible_solver;
2599 sft->getparam = slv1_get_parameters;
2600 sft->setparam = slv1_set_parameters;
2601 sft->getstatus = slv1_get_status;
2602 sft->solve = slv1_solve;
2603 sft->presolve = slv1_presolve;
2604 sft->iterate = slv1_iterate;
2605 sft->resolve = slv1_resolve;
2606 sft->getlinsol = slv1_get_linsol_sys;
2607 sft->getlinsys = slv1_get_linsolqr_sys;
2608 sft->getsysmtx = slv1_get_jacobian;
2609 sft->dumpinternals = slv1_dump_internals;
2610 return 0;
2611 }
2612 #endif /* #else clause of DYNAMIC_MINOS */
2613 #endif /* #else clause of !STATIC_MINOS && !DYNAMIC_MINOS */

Properties

Name Value
svn:executable *

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