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