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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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