/[ascend]/branches/shaun/disused/solver/slv1.c
ViewVC logotype

Contents of /branches/shaun/disused/solver/slv1.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2825 - (show annotations) (download) (as text)
Tue Feb 17 06:57:06 2015 UTC (7 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 82426 byte(s)
branch for shaun

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

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