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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 56414 byte(s)
Setting up web subdirectory in repository
1 /*
2 * opt rsqp
3 * by Dave Ternet and Ben Allan
4 * Created: 4/98
5 * Version: $Revision: 1.38 $
6 * Version control file: $RCSfile: slv2.c,v $
7 * Date last modified: $Date: 2000/01/25 02:27:26 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1998 Carnegie Mellon University
13 *
14 * The SLV solver is free software; you can redistribute
15 * it and/or modify it under the terms of the GNU General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * The SLV solver is distributed in hope that it will be
20 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with the program; if not, write to the Free Software Foundation,
26 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
27 * COPYING. COPYING is found in ../compiler.
28 *
29 */
30
31 #include <math.h>
32 #include <stdarg.h>
33 #include "utilities/ascConfig.h"
34 #include "utilities/ascSignal.h"
35 #include "utilities/ascMalloc.h"
36 #include "utilities/ascPanic.h"
37 #include "utilities/set.h"
38 #include "general/tm_time.h"
39 #include "utilities/mem.h"
40 #include "compiler/compiler.h"
41 #include "general/list.h"
42 #include "compiler/fractions.h"
43 #include "compiler/dimen.h"
44 #include "compiler/functype.h"
45 #include "compiler/func.h"
46 #include "solver/mtx.h"
47 #include "solver/linsol.h"
48 #include "solver/linsolqr.h"
49 #include "solver/slv_types.h"
50 #include "solver/var.h"
51 #include "solver/rel.h"
52 #include "solver/discrete.h"
53 #include "solver/conditional.h"
54 #include "solver/logrel.h"
55 #include "solver/bnd.h"
56 #include "solver/calc.h"
57 #include "solver/relman.h"
58 #include "solver/slv_common.h"
59 #include "solver/slv_client.h"
60 #include "solver/slv_stdcalls.h"
61 #include "solver/rsqp.h"
62 #include "solver/slv2.h"
63
64 /*
65 * name of this solver.
66 */
67
68 #define NAME "rSQP"
69
70 /*
71 * var/rel list reorderings to yield nice patterns for ma28
72 * may be indicated by this enum.
73 * Should be translated to a string list parameter later.
74 * not all are yet implemented.
75 */
76 enum slv2_reorderings {
77 s2_natural, /* rank not checked! more or less MODEL ordering. */
78 s2_spk1, /* partition/reorder arbitrary basis. */
79 s2_tspk1, /* partition/reorder arbitrary basis. */
80 s2_teardrop, /* partition/reorder arbitrary basis. */
81 s2_tteardrop, /* partition/reorder arbitrary basis. */
82 s2_teardropfat, /* partition/reorder (double border) arbitrary basis. */
83 s2_tteardropfat, /* partition/reorder (double border) arbitrary basis. */
84 s2_numeric, /* partition/teardrop/unpartitioned lu for pivot size. */
85 s2_random /* yeah, eat that */
86 };
87
88 #if !defined(STATIC_OPTSQP) && !defined(DYNAMIC_OPTSQP)
89 int slv2_register(SlvFunctionsT *f)
90 {
91 (void)f; /* stop gcc whine about unused parameter */
92
93 FPRINTF(stderr,"%s not compiled in this ASCEND IV.\n",NAME);
94 return 1;
95 }
96 #else /* either STATIC_OPTSQP or DYNAMIC_OPTSQP is defined */
97
98 #ifdef DYNAMIC_OPTSQP
99
100 /* do dynamic loading stuff. yeah, right */
101
102 #else /* following is used if STATIC_OPTSQP is defined */
103
104
105
106 /**************** SUBPARAMETERS FOR OPT_ASCEND **************************
107
108 Subparameter implemented: (value/meaning)
109
110 EPS Convergence tolerance on K-T error, default = 10e-7
111 MAXIT number of iterations allowed.
112 KPRINT Print output -> default = 0
113 IDEBUG 0 -> no debug info printed
114 1-9 -> debug info printed
115 IOPTION 1 -> linesearch only
116 2 -> trustregion only
117 3 -> linesearch and trustregion
118 ICHOOSE 0 -> MA28 basis
119 1 -> first m variables are independent
120 2 -> last m variables are independent
121 3 -> get independent vars from d.dat
122 4 -> get dependent vars from d.dat
123 ICORR 0 -> No correction term
124 1 -> Finite difference using multipliers
125 2 -> Finite difference without multipliers
126 3 -> Broyden update (not available)
127 4 -> Limited memory Broyden (not available)
128 5 -> Interactive correction
129 I_MULT_FREE 0 -> with multipliers
130 1 -> without multipliers
131 IIEXACT 0 -> Exact penalty function
132 1 -> Augmented Lagrangian
133 2 -> Watchdog line search (for Maratos effect)
134
135 IVV 0.0 -> use default trust region multipliers, VV=1.0
136 value -> use input trust region multiplier, VV
137 IDD 0 -> use default trust region shape
138 1 -> use input trust region shape, DD
139
140 IPMETHOD 0 -> use active-set method, QPKWIK
141 1 -> use interior-point method, QPIP
142 ILUSOLVE 271 -> use barrier27full()
143 272 -> use barrier27aug()
144 280 -> use barrier28big()
145 281 -> use barrier28full()
146 282 -> use barrier28aug()
147 471 -> use barrier47full()
148 472 -> use barrier47aug()
149 480 -> use barrier48big()
150 481 -> use barrier48full() (default)
151 482 -> use barrier48aug()
152 483 -> use pre_qpip48new() [try 481 first, then 484]
153 484 -> use pre_qpip48infeas()
154 485 -> use barrier48pd()
155 486 -> use barrier48mpc()
156 487 -> use barrier48ob1()
157 488 -> use barrier48gondzio()
158 IWARM 0 -> don't use warm start option for QP solver
159 1 -> use warm start option for QP solver
160
161 CNR 0 -> no output to file
162 1 -> print extra output to file 'journal.num'
163 TIMEMAX number seconds allowed.
164 ***********************************************************************/
165
166 #define SERVER (sys->slv)
167
168 /* updating to Ken's parameter array implementation which is more flexible */
169 #define slv2_PA_SIZE 17
170
171 #define EPS_PTR (sys->parm_array[0])
172 #define EPS ((*(real64 *)EPS_PTR))
173 #define MAXIT_PTR (sys->parm_array[1])
174 #define MAXIT ((*(int32 *)MAXIT_PTR))
175 #define KPRINT_PTR (sys->parm_array[2])
176 #define KPRINT ((*(int32 *)KPRINT_PTR))
177 #define IDEBUG_PTR (sys->parm_array[3])
178 #define IDEBUG ((*(int32 *)IDEBUG_PTR))
179 #define IOPTION_PTR (sys->parm_array[4])
180 #define IOPTION ((*(int32 *)IOPTION_PTR))
181 #define ICHOOSE_PTR (sys->parm_array[5])
182 #define ICHOOSE ((*(int32 *)ICHOOSE_PTR))
183 #define ICORR_PTR (sys->parm_array[6])
184 #define ICORR ((*(int32 *)ICORR_PTR))
185 #define I_MULT_FREE_PTR (sys->parm_array[7])
186 #define I_MULT_FREE ((*(int32 *)I_MULT_FREE_PTR))
187 #define IIEXACT_PTR (sys->parm_array[8])
188 #define IIEXACT ((*(int32 *)IIEXACT_PTR))
189 #define IVV_PTR (sys->parm_array[9])
190 #define IVV ((*(real64 *)IVV_PTR))
191 #define IDD_PTR (sys->parm_array[10])
192 #define IDD ((*(int32 *)IDD_PTR))
193 #define IPMETHOD_PTR (sys->parm_array[11])
194 #define IPMETHOD ((*(int32 *)IPMETHOD_PTR))
195 #define ILUSOLVE_PTR (sys->parm_array[12])
196 #define ILUSOLVE ((*(char **)ILUSOLVE_PTR))
197 #define IWARM_PTR (sys->parm_array[13])
198 #define IWARM ((*(int32 *)IWARM_PTR))
199 #define CNR_PTR (sys->parm_array[14])
200 #define CNR ((*(int32 *)CNR_PTR))
201 #define TIMEMAX_PTR (sys->parm_array[15])
202 #define TIMEMAX ((*(real64 *)TIMEMAX_PTR))
203 /* some ascend side stuff */
204 #define CUTOFF_PTR (sys->parm_array[16])
205 #define CUTOFF ((*(int32 *)CUTOFF_PTR))
206
207
208 #define D_ZERO (double)0.0
209 #define D_ONE (double)1.0
210
211 #define SLV2(a) ((slv2_system_t)(a))
212
213 struct slv2_system_structure {
214
215 /* initialized in slv2_create */
216
217 int integrity;
218 slv_parameters_t p; /* the standard slv parameters and status */
219 void *parm_array[slv2_PA_SIZE];
220 struct slv_parameter padata[slv2_PA_SIZE];
221 enum slv2_reorderings reorderstyle; /* ascend matrix ordering */
222 slv_status_t s;
223 double clock; /* cumulative cpu */
224
225 rel_filter_t relfilter, *rfilter;
226 var_filter_t varfilter, *vfilter;
227
228 /* Problem definition from server */
229 slv_system_t slv; /* slv_system_t back-link */
230 struct rel_relation *obj; /* Objective function: NULL = none */
231 struct var_variable **vlist; /* Variable list */
232 struct rel_relation **rlist; /* Relation list */
233
234 /* computed at presolve */
235 int32 rlen; /* number of relations opt cares about. first in rlist */
236 int32 vlen; /* number of variables opt cares about. first in vlist */
237
238 int32 maxndx;
239
240 /* push space for common block controls. need to match rsqp.h sizes. */
241 real64 rcontrol[RCONTROLSIZE];
242 int32 icontrol[ICONTROLSIZE];
243
244 /* Arrays allocated at presolve */
245 real64 *gsparse; /* space to store obj gradient from relman */
246 int32 *gxsparse; /* C indices of var in gsparse from relman. */
247 int32 gnsparse; /* size of g/gxsparse. */
248
249 int32 *iw; /* integer workspace */
250
251 real64
252 *rw, /* real workspace */
253 *xn, /* variable values scaled by nominals*/
254 *nominals, /* variable scaling from presolve. */
255 *bnds, /* bound limits */
256 *a, /* sparse jacobian, scaled by nominals, weights */
257 *g, /* object gradient scaled by nominals, objwt */
258 *c, /* residual scaled by weights */
259 *weights, /* relation scalings. */
260 *trust;
261 real64 f, objwt; /* objective and its gradient scaling */
262
263 int32 *avar, *acon; /* column indices and row indices of a */
264 int32 *ivec, *jvec; /* jvec and ivec for relman_diff_harwell */
265
266 int32
267 md,
268 nzd,
269 ibnds,
270 nfunc,
271 ngrad,
272 iwsize,
273 rwsize;
274
275 int32 inform, /* opt status? */
276 iluoption, /* ILUSOLVE translated */
277 njac, /* # jacobian elements (linear+nonlinear) */
278 iter; /* number of iterations */
279 };
280
281 /**
282 *** Integrity checks
283 *** ----------------
284 *** check_system(sys)
285 **/
286
287 #define OK ((int)492031596)
288 #define DESTROYED ((int)729104829)
289 /* note deaddead would be int 3735936685 */
290 static int check_system(slv2_system_t sys)
291 /* Checks sys for NULL and for integrity. */
292 {
293
294 if( sys == NULL ) {
295 FPRINTF(stderr,"ERROR: (slv2) check_system\n");
296 FPRINTF(stderr," NULL system handle.\n");
297 return 1;
298 }
299
300 switch( sys->integrity ) {
301 case OK:
302 return 0;
303 case DESTROYED:
304 FPRINTF(stderr,"ERROR: (slv2) check_system\n");
305 FPRINTF(stderr," System was recently destroyed.\n");
306 return 1;
307 default:
308 FPRINTF(stderr,"ERROR: (slv2) check_system\n");
309 FPRINTF(stderr," System reused or never allocated.\n");
310 return 1;
311 }
312 }
313
314 /*********************************************************************
315 *** Memory management routines ***
316 free_unless_null(ptr)
317 zero(arr,nelts,type)
318 alloc_vector(len)
319 alloc_ivector(len)
320 free_vector(vec)
321 *********************************************************************/
322
323 static void free_unless_null(POINTER ptr)
324 {
325 if( ptr != NULL )
326 ascfree(ptr);
327 }
328
329 static void zerodoublesF(real64 *a, int i)
330 {
331 for (; i > 0; ) {
332 i--;
333 a[i] = 0.0;
334 }
335 }
336
337 static void zerointsF(int32 *a, int i)
338 {
339 for (; i > 0; ) {
340 i--;
341 a[i] = 0;
342 }
343 }
344
345 static real64 *alloc_zero_vector(int len)
346 {
347 double *result;
348 result = (real64 *)ascmalloc(len*sizeof(real64));
349 if (result != NULL) {
350 zerodoublesF(result,len);
351 }
352 return result;
353 }
354
355 static int32 *alloc_zero_ivector(int len)
356 {
357 int32 *result;
358 result = (int32 *)ascmalloc(len*sizeof(int32));
359 if (result != NULL) {
360 zerointsF(result,len);
361 }
362 return result;
363 }
364
365 /*
366 * init an array to 0.
367 */
368 #define zerodoubles(arr,nelts,type) zerodoublesF(arr,nelts)
369 /*
370 * create real and integer arrays
371 */
372 real64 *alloc_vector(int len)
373 {
374 real64 *result;
375 result = (real64 *)ascmalloc((len+1)*sizeof(real64));
376 result[len] = OK;
377 return result;
378 }
379 int32 *alloc_ivector(int len)
380 {
381 int32 *result;
382 result = (int32 *)ascmalloc((len+1)*sizeof(int32));
383 result[len] = OK;
384 return result;
385 }
386
387 /* destroy the return of both previous alloc_s */
388 #define free_vector(vec) free_unless_null((POINTER)(vec))
389
390 /*********************************************************************
391 *** General I/O routines ***
392 print_var_name(out,sys,var)
393 print_rel_name(out,sys,rel)
394 *********************************************************************/
395
396 static void print_vector (slv2_system_t sys, double * vec,int len, char *name)
397 {
398 int i;
399 for (i=0; i<len; i++)
400 FPRINTF(LIF(sys),"%s(%d) == %20.15g\n",name,i,vec[i]);
401 FFLUSH(LIF(sys));
402 }
403
404 static void print_ivector (slv2_system_t sys, int * vec,int len, char *name) {
405 int i;
406 for (i=0; i<len; i++)
407 FPRINTF(LIF(sys),"%s(%d) == %d\n",name,i,vec[i]);
408 FFLUSH(LIF(sys));
409 }
410
411 #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
412 #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
413
414 /*********************************************************************
415 ** routines to work around the insanities of common blocks. **
416 *********************************************************************/
417
418 /* Every entry from an interface must do these or there will be
419 * cross chatter among coexisting solvers such as occur in a
420 * qp and an nlp alternating. Don't want the qp confusing the nlp.
421 * The basic idea is that each slv2_system should push the common
422 * info already existing into local storage and then set its own.
423 * Substantial speed increase for using slv2_solve instead of iterate.
424 * Dave alleges there is no internal common other than these which
425 * persist across iterations. In theory this means we only
426 * need to call initoptcommon at register time.
427 */
428 void push_controls(int32 *icontrol, real64 *rcontrol)
429 {
430 /* save the old controls to a local vector */
431 GET_OPT_COMMON(rcontrol,
432 &(icontrol[0]),
433 &(icontrol[1]),
434 &(icontrol[2]),
435 &(icontrol[3]),
436 &(icontrol[4]),
437 &(icontrol[5]),
438 &(icontrol[6]),
439 &(icontrol[7]),
440 &(icontrol[8]),
441 &(icontrol[9]),
442 &(icontrol[10]),
443 &(icontrol[11])
444 );
445 }
446 void pop_controls(int32 *icontrol, real64 *rcontrol)
447 {
448 /* restore the old controls from a local vector */
449 SET_OPT_COMMON(rcontrol,
450 &(icontrol[0]),
451 &(icontrol[1]),
452 &(icontrol[2]),
453 &(icontrol[3]),
454 &(icontrol[4]),
455 &(icontrol[5]),
456 &(icontrol[6]),
457 &(icontrol[7]),
458 &(icontrol[8]),
459 &(icontrol[9]),
460 &(icontrol[10]),
461 &(icontrol[11])
462 );
463 }
464
465
466 /*********************************************************************
467 bookkeeping functions
468 iteration_begins(sys)
469 iteration_ends(sys)
470 *********************************************************************/
471
472 static void iteration_begins(slv2_system_t sys)
473 /* Prepares sys for an iteration, increasing the iteration counts
474 and starting the clock. */
475 {
476 sys->clock = tm_cpu_time();
477 }
478
479 static void iteration_ends(slv2_system_t sys)
480 /* Prepares sys for exiting an iteration, stopping the clock and recording
481 * the cpu time, as well as updating the status.
482 */
483 {
484 double cpu_elapsed; /* elapsed this iteration */
485 boolean unsuccessful;
486
487 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
488 sys->s.cpu_elapsed += cpu_elapsed;
489
490 if( !sys->s.converged ) {
491 sys->s.block.cpu_elapsed += cpu_elapsed;
492 sys->s.time_limit_exceeded =
493 (sys->s.block.cpu_elapsed > sys->p.time_limit);
494 }
495
496 unsuccessful = sys->s.diverged ||
497 sys->s.inconsistent ||
498 sys->s.iteration_limit_exceeded ||
499 sys->s.time_limit_exceeded;
500 /* maybe add this line in a little later
501 * (sys->vlen <= sys->rlen);
502 */
503 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
504 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
505
506 }
507
508
509 /* Residuals are calculated, norm value is stored in status, and
510 * sys->s.calc_ok is set accordingly. Knows nothing about scaling.
511 * this is out of date with slv3 standard, if that's a standard...
512 * This function should not be used for the solution process.
513 */
514 static void calc_residuals(slv2_system_t sys)
515 {
516 real64 sum,res;
517 struct rel_relation **rp;
518
519 calc_ok = TRUE;
520 sum = D_ZERO;
521 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
522 for( rp=sys->rlist ; *rp != NULL ; ++rp ) {
523 if (rel_included(*rp) && rel_active(*rp)) {
524 res = relman_eval(*rp,&calc_ok,1/* safe always*/);
525 if( !relman_calc_satisfied(*rp,sys->p.tolerance.feasible) ||
526 ( !rel_less(*rp) && !rel_greater(*rp) )
527 ) {
528 sum += calc_sqr_D0(res);
529 }
530 }
531 }
532 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
533 if (sum > D_ZERO) {
534 sys->s.block.residual = calc_sqrt_D0(sum);
535 } else {
536 sys->s.block.residual = sum;
537 }
538
539 }
540
541 static void install_allvars(slv2_system_t sys, real64 *values)
542 /*********************************************************************
543 Moves variables from given array to the value field of each variable
544 so that ascend can do evaluation..
545 *********************************************************************/
546 {
547 int32 col;
548 struct var_variable *var;
549 for( col=0 ; col < sys->vlen ; ++col ) {
550 var = sys->vlist[col];
551 var_set_value(var,values[col]*sys->nominals[col]);
552 }
553 }
554
555
556
557 /*********************************************************************
558 *** MODEL server utilities for OPT_ASCEND interface ***
559 *********************************************************************/
560
561 /*
562 * slv2_calc_C -- constraint residuals.
563 * slv2_calc_A -- constraint jacobian.
564 * slv2_calc_F -- objective.
565 * slv2_calc_G -- objective gradient.
566 *
567 * All these functions should automagically scale the problem
568 * from sys->nominals and sys->wts which should
569 * be calculated only at the beginning or after throwing
570 * out the accumulated Hessian following a huge error reduction.
571 *
572 * all return 0 if able to evaluate at sys->xn or nonzero
573 * if some error was seen.
574 */
575
576 static int slv2_calc_C(slv2_system_t sys)
577 {
578 int32 row;
579
580 calc_ok = TRUE;
581 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
582 for( row = 0; row < sys->rlen; row++ ) {
583 sys->c[row] =
584 relman_eval(sys->rlist[row],&calc_ok,1) * sys->weights[row];
585 }
586 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
587
588 if( !calc_ok ) {
589 FPRINTF(MIF(sys),"!Warning:residual calculation error(s).\n");
590 return 1;
591 }
592 return 0;
593 }
594
595 /*
596 * computes scaled gradient a. if firsttime != 0, also
597 * fill sys->acon, sys->avar, sys->ivec, sys->jvec, sys->weights.
598 * return 0 if ok, -1 if severe error, n > 0 equation errors.
599 */
600 static int slv2_calc_A(slv2_system_t sys, int firsttime)
601 {
602 int32 k, row;
603 real64 nominal;
604 int result=0;
605
606 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
607 if (!firsttime) {
608 result =
609 relman_diff_harwell(sys->rlist, sys->vfilter,sys->rfilter,sys->rlen,0,3,
610 sys->a,NULL,NULL);
611 if (result == 1) {
612 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
613 return -1;
614 }
615 if (result == 0) {
616 for (k = sys->njac-1; k >= 0; k--) {
617 /* we computed nominals and weights so this is multiplication */
618 sys->a[k] =
619 sys->a[k] * sys->weights[sys->ivec[k]] * sys->nominals[sys->jvec[k]];
620 }
621 } else {
622 result = 1;
623 FPRINTF(MIF(sys),"!Warning: jacobian calculation error(s).\n");
624 for (k = sys->njac-1; k >= 0; k--) {
625 if (isnan(sys->a[k]) || !finite(sys->a[k])) {
626 sys->a[k] = 0; /* fp error --> singular row probably */
627 }
628 sys->a[k] = sys->a[k]
629 * sys->weights[sys->ivec[k]] * sys->nominals[sys->jvec[k]];
630 }
631 }
632 } else {
633 /* first time through compute weights, ivec, acon, jvec, avar */
634 result =
635 relman_diff_harwell(sys->rlist, sys->vfilter,sys->rfilter,sys->rlen,0,3,
636 sys->a,sys->ivec,sys->jvec);
637 if (result == 1) {
638 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
639 return -1;
640 }
641 if (result < 0) {
642 FPRINTF(MIF(sys),
643 "!Warning: jacobian calculation error(s) at initial point.\n");
644 result = 1;
645 }
646 if (/* 2norm option */ 1) {
647 k = 0;
648 for (row = 0 ; row < sys->rlen; row++) {
649 nominal = 0;
650 while (sys->ivec[k] == row) {
651 nominal += sys->a[k]*sys->a[k];
652 k++;
653 }
654 nominal = sqrt(nominal);
655 if (!isnan(nominal) && finite(nominal) && nominal != 0.0) {
656 sys->weights[row] = 1/nominal;
657 } else {
658 result++;
659 sys->weights[row] = 1;
660 }
661 }
662 } else {
663 /* ascend relation scaling */
664 for (row = 0 ; row < sys->rlen; row++) {
665 nominal = relman_scale(sys->rlist[row]);
666 nominal = sqrt(nominal);
667 if (!isnan(nominal) && finite(nominal) && nominal != 0.0) {
668 sys->weights[row] = 1/nominal;
669 } else {
670 result++;
671 sys->weights[row] = 1;
672 }
673 }
674 }
675 for (k = sys->njac-1; k >= 0; k--) {
676 sys->acon[k] = sys->ivec[k] + 1;
677 sys->avar[k] = sys->jvec[k] + 1;
678 sys->a[k] = sys->a[k]
679 * sys->weights[sys->ivec[k]] * sys->nominals[sys->jvec[k]];
680 }
681 }
682 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
683 FPRINTF(LIF(sys),"jacobian calculated\n");
684 return result;
685 }
686
687 /* return 0 if ok 1 if error */
688 static int slv2_calc_F(slv2_system_t sys)
689 {
690 int result=0;
691 if( sys->obj != NULL ) {
692 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
693 sys->f = relman_eval(sys->obj,&calc_ok,1) * sys->objwt;
694 if (isnan(sys->f) && !finite(sys->f)) {
695 result = 1;
696 }
697 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
698 }
699 sys->s.calc_ok = calc_ok;
700 return result;
701 }
702
703 static int slv2_calc_G(slv2_system_t sys, int firsttime)
704 {
705 int status, ndx, result = 0;
706 real64 nominal;
707
708 if ( sys->obj != NULL ) {
709 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
710 status = relman_diff2(sys->obj, sys->vfilter, sys->gsparse,
711 sys->gxsparse, &(sys->gnsparse),0);
712 if (firsttime) {
713 if (/*2norm option */1) {
714 nominal = 0;
715 for (ndx = 0; ndx < sys->gnsparse; ndx++) {
716 nominal += sys->gsparse[ndx] * sys->gsparse[ndx];
717 }
718 nominal = sqrt(nominal);
719 if (!isnan(nominal) && finite(nominal) && nominal != 0.0) {
720 sys->objwt = 1/nominal;
721 } else {
722 sys->objwt = 1;
723 result++;
724 }
725 } else {
726 nominal = relman_scale(sys->obj);
727 nominal = sqrt(nominal);
728 if (!isnan(nominal) && finite(nominal) && nominal != 0.0) {
729 sys->objwt = 1/nominal;
730 } else {
731 result++;
732 sys->objwt = 1;
733 }
734 }
735 }
736
737 for (ndx = 0; ndx < sys->gnsparse; ndx++) {
738 if (!isnan(sys->gsparse[ndx]) && finite(sys->gsparse[ndx])) {
739 sys->g[sys->gxsparse[ndx]] = sys->gsparse[ndx]
740 * sys->objwt * sys->nominals[sys->jvec[sys->gxsparse[ndx]]];
741 } else {
742 sys->g[sys->gxsparse[ndx]] = 0;
743 }
744
745 FPRINTF(LIF(sys),"d(obj)/d(");
746 print_var_name(LIF(sys),sys,sys->vlist[ndx]);
747 FPRINTF(LIF(sys),") = %20.16g\n",sys->g[ndx]);
748 FPRINTF(LIF(sys),"g (f77) subscript stuffed: %d\n",
749 sys->gxsparse[ndx]+1);
750
751 }
752 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
753 }
754 if( !calc_ok || result) {
755 FPRINTF(MIF(sys),"!Warning: objective gradient calculation error(s).\n");
756 return 1;
757 }
758 return 0;
759 }
760
761
762 /*********************************************************************
763 *** OPT_ASCEND interface ***
764 int scale_and_bounds(sys)
765 invoke_opt(sys)
766 *********************************************************************/
767
768
769 #define THEIR_INFINITY (double)1.0e20
770
771 /* Init bnds, xn, acon, avar, a.
772 * sys->nominals is found.
773 * bounds are scaled by the var nominals after checking ascend sanity.
774 * init point (scaled by nominals) is moved inside theirinfinity bounds.
775 * returns nonzero if gradient errors seen at initial point.
776 */
777 static
778 int scale_and_bounds(slv2_system_t sys)
779 {
780 int32 col,ndx;
781 real64 nominal;
782 int result=0;
783 struct var_variable **vp;
784 struct var_variable *var;
785
786 slv_insure_bounds(SERVER,0,sys->vlen,MIF(sys));
787
788 for( col=0 ; col < sys->vlen ; ++col ) {
789 ndx = 2*col;
790 var = sys->vlist[col];
791
792 sys->xn[col] = var_value(var)/nominal;
793
794 nominal = fabs(var_nominal(var));
795 sys->nominals[ndx] = nominal;
796 #define NEWTONLIMIT 0.1
797 /* Needs to be moved to parameters list.
798 * May want to move trust to solver_var instead of using this calc.
799 */
800 sys->trust[ndx] = nominal * NEWTONLIMIT;
801
802 sys->bnds[ndx] = var_lower_bound(var)/nominal;
803 if (sys->bnds[ndx] <= -THEIR_INFINITY) {
804 sys->bnds[ndx] = -THEIR_INFINITY;
805 print_var_name(LIF(sys),sys,var);
806 FPRINTF(LIF(sys), " has %s lower bound (%g) tighter than ASCEND's.\n",
807 NAME,(-THEIR_INFINITY)*nominal);
808 }
809
810 ndx++;
811 sys->bnds[ndx] = var_upper_bound(var)/nominal;
812 if (sys->bnds[ndx] >= THEIR_INFINITY) {
813 sys->bnds[ndx] = THEIR_INFINITY;
814 print_var_name(LIF(sys),sys,var);
815 FPRINTF(LIF(sys)," %s upper bound (%g) tighter than ASCEND's.\n",
816 NAME, THEIR_INFINITY*nominal);
817 }
818 }
819
820 /* init acon, avar, ivec, jvec, weights */
821 if (slv2_calc_A(sys,1)) {
822 result++;
823 }
824
825 /* init objwt */
826 if (slv2_calc_G(sys,1)) {
827 result++;
828 }
829
830 return result;
831 }
832
833 static void write_inform(slv2_system_t sys)
834 /* reports opt exit information */
835 {
836 char * res=NULL; /* do not free this ptr ever */
837 FILE * fp=MIF(sys);
838 FPRINTF(fp,"\n ***** OPT_ASCEND inform= %d ***** \n ",sys->inform);
839 switch (sys->inform) {
840 case -10: res="End of iteration.";
841 break;
842 case -1 : res="Ready to start a problem.";
843 break;
844 case 0 : res="Starting in OPT_ASCEND";
845 break;
846 case 1 : res="Optimal solution found.";
847 break;
848 case 6 : res="User requested termination?";
849 break;
850 case 10 : res="Optimal solution to square system found.";
851 break;
852 case 42 : res="Need to change the size of iw and/or rw";
853 break;
854 case 55 : res="There is a problem with the bounds on a variable";
855 sys->s.diverged = TRUE;
856 break;
857 case 57 : res="Basis isn't chosen correctly, check ICHOOSE";
858 sys->s.diverged = TRUE;
859 break;
860 case 58 : res="Maximum number of iterations exceeded";
861 sys->s.iteration_limit_exceeded = TRUE;
862 break;
863 case 73 : res="Line search failure";
864 sys->s.diverged = TRUE;
865 break;
866 case 78 : res="singular pivot in correction.f";
867 sys->s.diverged = TRUE;
868 break;
869 case 101 :
870 case 102 :
871 case 103 :
872 case 104 :
873 case 105 :
874 case 106 :
875 case 107 :
876 case 108 :
877 case 109 :
878 case 110 :
879 case 111 :
880 case 112 :
881 case 113 :
882 case 114 :
883 case 115 :
884 case 116 :
885 case 117 :
886 case 118 :
887 case 119 :
888 case 120 : res="Iarray or Rarray sizes too small";
889 sys->s.diverged = TRUE;
890 break;
891 case 177 :
892 case 178 : res="null pointer returned allocing iw or rw";
893 sys->s.diverged = TRUE;
894 break;
895 default : res="Unknown inform returned from OPT_ASCEND";
896 sys->s.diverged = TRUE;
897 break;
898 }
899 FPRINTF(fp,"%s\n",res);
900 return;
901 }
902
903 /*
904 * returns 1 if we have not hit some opt or user stopping condition
905 */
906 static
907 int go_again(slv2_system_t sys)
908 {
909 if (Solv_C_CheckHalt()) {
910 sys->inform = 6;
911 }
912 return ((sys->inform != 1 && /* optimal solution found */
913 sys->inform != 6 && /* user interrupt */
914 sys->inform != 10 && /* opt. sol. to sqare sys. found */
915 sys->inform != -10 && /* end of one iteration */
916 sys->inform != 57 && /* ichoose incorrect */
917 sys->inform != 58 && /* max iter. reached */
918 sys->inform != 73 && /* line search failure */
919 sys->inform <= 100 ) || /* workspace sizes too small */
920 sys->inform == 42); /* reallocate memory */
921 }
922
923
924 /* Sets up the interface to OPT_ASCEND and then invokes it. */
925 static int invoke_opt(slv2_system_t sys)
926 {
927 int i;
928 real64 fout,e1out,e2out,ynout,znout,alfaout;
929 int32 nactout;
930
931
932
933 if (!sys->s.ready_to_solve) return (1); /*fail if not presolved*/
934
935 print_vector(sys,sys->xn,sys->vlen,"xn");
936 print_vector(sys,sys->a,sys->njac,"a");
937 print_vector(sys,sys->bnds,2*sys->vlen,"bnds");
938 print_ivector(sys,sys->avar,sys->njac, "avar");
939 print_ivector(sys,sys->acon,sys->njac, "acon");
940
941 iteration_begins(sys);
942
943 do {
944 if( sys->inform == -1 /* first call from ascend after presolve. */ ||
945 sys->inform == 0 ||
946 sys->inform == -10 ||
947 sys->inform == -15 ||
948 sys->inform == -25 ) {
949 sys->ngrad = sys->ngrad + 1;
950 slv2_calc_G(sys,0);
951 slv2_calc_A(sys,0);
952 }
953
954 if( sys->inform == -1 ||
955 sys->inform == 0 ||
956 sys->inform == -20 ||
957 sys->inform == -25 ) {
958 sys->nfunc = sys->nfunc + 1;
959 slv2_calc_F(sys);
960 slv2_calc_C(sys);
961 }
962
963 if ( IDEBUG != 0 ) {
964 PRINTF("\nn=%3d\tm=%3d\tmd=%3d\n", sys->vlen,sys->rlen,sys->md);
965 if ( sys->vlen <= 100 ) {
966 for(i=0; i<sys->vlen; i++) {
967 PRINTF("\ti=%3d\tlo=%g\tx[i]=%g\thi=%g\n",
968 i,sys->bnds[2*i],sys->xn[i],sys->bnds[2*i+1]);
969 }
970 for(i=0; i<sys->vlen; i++) {
971 PRINTF("\tg[%2d]=%g", i, sys->g[i]);
972 }
973 }
974 PRINTF("\nf=%f\tnjac=%3d\tnzd=%3d\n", sys->f,sys->njac,sys->nzd);
975
976 PRINTF("\n");
977 if ( sys->rlen <= 100 ) {
978 for(i=0; i<sys->rlen; i++) {
979 PRINTF("\tc[%2d]=%g", i, sys->c[i]);
980 }
981 }
982
983 PRINTF("\n");
984 if ( sys->njac <= 100 ) {
985 for(i=0; i<sys->njac; i++) {
986 PRINTF("\ta[%3d]=%g", i, sys->a[i]);
987 }
988 }
989 PRINTF("\nibnds=%2d\tinf=%.8f\n", sys->ibnds,sys->inform);
990 PRINTF("\nliw=%d\tlrw=%d", sys->iwsize,sys->rwsize);
991 }
992
993 OPT_ASCEND(&(sys->vlen), /* n */
994 &(sys->rlen), /* m */
995 &(sys->md), /* md */
996 (sys->xn), /* x */
997 &(sys->f), /* f */
998 (sys->g), /* g */
999 (sys->c), /* c */
1000 (sys->a), /* a */
1001 &(sys->njac), /* nz */
1002 &(sys->nzd), /* nzd */
1003 (sys->avar), /* avar */
1004 (sys->acon), /* acon */
1005 &(sys->ibnds), /* ibnds */
1006 (sys->bnds), /* bnds */
1007 &(sys->inform), /* inf */
1008 (sys->trust), /* trust */
1009 &(sys->iwsize), /* liw */
1010 (sys->iw), /* iw */
1011 &(sys->rwsize), /* lrw */
1012 (sys->rw)); /* rw */
1013
1014 if (sys->inform==42) {
1015 free_vector(sys->iw);
1016 free_vector(sys->rw);
1017 PRINTF("\niwsize = %d, rwsize = %d\n", sys->iwsize, sys->rwsize);
1018 sys->iw = alloc_ivector(sys->iwsize);
1019 if ( sys->iw == NULL ) {
1020 PRINTF("\nsys->iw is null, too much memory required");
1021 sys->inform == 177;
1022 }
1023 sys->rw = alloc_vector(sys->rwsize);
1024 if ( sys->rw == NULL ) {
1025 PRINTF("\nsys->rw is null, too much memory required");
1026 sys->inform == 178;
1027 }
1028 sys->inform = -1;
1029 FPRINTF(LIF(sys),
1030 "OPT_ASCEND workspace %d bytes.\n",(sys->rwsize*sizeof(real64)));
1031 FPRINTF(LIF(sys),
1032 "OPT_ASCEND workspace %d bytes.\n",(sys->iwsize*sizeof(int32)));
1033 continue;
1034 }
1035 GET_OPT_OUTPUT(&fout,&e1out,&e2out,&ynout,&znout,&alfaout,&nactout);
1036 FPRINTF(MIF(sys),"Iter\tObjective\tmax|glag|\tmax|Ci|\tYpYN\tZpZN\tNACT\talfa\n");
1037 FPRINTF(MIF(sys),"%d\t%.12g\t%g\t%g\t%g\t%g\t%d%g\n",
1038 sys->iter,fout,e1out,e2out,ynout,znout,nactout,alfaout);
1039
1040 } while ( go_again(sys) );
1041
1042 iteration_ends(sys);
1043
1044 PRINTF("\niterations=%d\tfunc.evals=%d\tgrad.evals=%d",
1045 sys->iter,sys->nfunc,sys->ngrad);
1046
1047 if (sys->inform == 1 || sys->inform == 10) {
1048 sys->s.converged = 1;
1049 sys->s.block.previous_total_size = sys->vlen;
1050 }
1051
1052 sys->s.iteration = sys->nfunc;
1053 sys->s.block.iteration = sys->ngrad;
1054
1055 print_vector(sys,sys->xn,sys->vlen,"x");
1056 print_vector(sys,sys->a,sys->njac,"a");
1057 print_vector(sys,sys->bnds,2*sys->vlen,"bnds");
1058
1059 install_allvars(sys,sys->xn);
1060
1061 write_inform(sys);
1062
1063 iteration_begins(sys);
1064 calc_residuals(sys);
1065 iteration_ends(sys);
1066 }
1067
1068
1069 /*********************************************************************
1070 *** External routines ***
1071 See slv.h and rsqp.h
1072 *********************************************************************/
1073
1074 static char *ilu_names[] = {
1075 "271-barrier27full()", "272-barrier27aug()", /*0,1*/
1076 "280-barrier28big()", "281-barrier28full()", "282-barrier28aug()", /*2,3,4*/
1077 "471-barrier47full()", "472-barrier27aug()", /*5,6*/
1078 "480-barrier48big()", "481-barrier48full()", "482-barrier48aug()", /*7,8,9*/
1079 "483-barrier48new()", "484-barrier48infeas()", /*10,11*/
1080 "485-barrier48pd()", "486-barrier48mpc()", /*12,13*/
1081 "487-barrier48ob1()", "488-barrier48gondzio()" /*14,15*/
1082 };
1083
1084 /* converts string input to matching integer */
1085 static
1086 int32 ilu_name_to_int(slv2_system_t sys)
1087 {
1088 int option, status;
1089 char *s;
1090
1091 s = ILUSOLVE;
1092
1093 status = sscanf(s,"%d",&option);
1094 if (status == 1) {
1095 /* matched leading integer */
1096 #ifndef __CMU__
1097 /* harwell source for ma27,ma47,ma48 not available outside cmu normally */
1098 if (option < 280 || option > 289) {
1099 if (sys->iluoption < 280 || sys->iluoption >289) {
1100 sys->iluoption = 280;
1101 }
1102 return sys->iluoption;
1103 }
1104 #endif
1105 sys->iluoption = option;
1106 }
1107 /* OTHERWISE leave default alone */
1108 return sys->iluoption;
1109 }
1110
1111 /* converts integer input to matching string index */
1112 /* returns 2 if nothing matched which should be 280 */
1113 static int32 ilu_name_from_int(int32 option)
1114 {
1115 int i;
1116 char buf[20];
1117
1118 sprintf(buf,"%d",option);
1119 for (i = 0; i < (sizeof(ilu_names)/sizeof(char *)); i++) {
1120 if (strncmp(buf,ilu_names[i],3)==0) {
1121 return i;
1122 }
1123 }
1124 return 2;
1125 }
1126
1127 /* this may be called with NULL, NULL, p */
1128 static
1129 int32 slv2_get_default_parameters(slv_system_t server, SlvClientToken asys,
1130 slv_parameters_t *parameters)
1131 {
1132 slv2_system_t sys; int isize; int rsize, dummy;
1133 real64 rcontrol[RCONTROLSIZE];
1134 int32 icontrol[ICONTROLSIZE];
1135 real64 rdefault[RCONTROLSIZE];
1136 int32 idefault[ICONTROLSIZE];
1137 union parm_arg lo,hi,val;
1138 struct slv_parameter *new_parms = NULL;
1139 int32 make_macros = 0;
1140
1141 /* here call init common, call get common, and set our defaults
1142 * from those defaults subject to ascend limits on some.
1143 * reset common block in each major call from the interface
1144 * from our defaults since we don't know if there are two
1145 * systems sharing the same common.
1146 */
1147
1148 push_controls(icontrol,rcontrol);
1149 INIT_OPT_DEFAULTS();
1150
1151 push_controls(idefault,rdefault); /* don't pop this one, stuff defaults. */
1152 if (server != NULL && asys != NULL) {
1153 sys = SLV2(asys);
1154 make_macros = 1;
1155 }
1156
1157 #ifndef NDEBUG /* keep purify from whining on UMR */
1158 lo.argr = hi.argr = val.argr = 0.0;
1159 #endif
1160
1161 if (parameters->parms == NULL) {
1162 new_parms = (struct slv_parameter *)
1163 ascmalloc((slv2_PA_SIZE)*sizeof(struct slv_parameter));
1164 if (new_parms == NULL) {
1165 pop_controls(icontrol,rcontrol);
1166 return -1;
1167 }
1168 parameters->parms = new_parms;
1169 parameters->dynamic_parms = 1;
1170 }
1171 parameters->num_parms = 0;
1172
1173 /*
1174 * define control parameters here.
1175 * BPARM is a boolean (or 0/1 int) parameter.
1176 * RPARM IS_A double parameter.
1177 * IPARM IS_A integer parameter (not 0/1).
1178 * CPARM IS_A array of strings parameter (rarely seen in fortran).
1179 */
1180
1181 slv_define_parm(parameters, real_parm,
1182 "EPS", "convergence tolerance (EPS)",
1183 "Convergence tolerance on K-T error",
1184 U_p_real(val,rdefault[0]),
1185 U_p_real(lo,1e-20),U_p_real(hi,10), 1);
1186 SLV_RPARM_MACRO(EPS_PTR,parameters);
1187
1188 slv_define_parm(parameters, int_parm,
1189 "MAXIT", "iteration limit (MAXIT)",
1190 "iteration limit (MAXIT)",
1191 U_p_int(val,idefault[0]),U_p_int(lo,1),U_p_int(hi,10000),1);
1192 SLV_IPARM_MACRO(MAXIT_PTR,parameters);
1193
1194 #ifndef __WIN32__ /* pc fortran cannot do prints to stdout (file 6) */
1195 slv_define_parm(parameters, bool_parm,
1196 "KPRINT", "print output (KPRINT)", "print output",
1197 U_p_bool(val,idefault[1]),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1198 SLV_BPARM_MACRO(KPRINT_PTR,parameters);
1199 #else
1200 /* force FALSE value */
1201 slv_define_parm(parameters, bool_parm,
1202 "KPRINT", "print output (KPRINT)", "print output",
1203 U_p_bool(val,0),U_p_bool(lo,0),U_p_bool(hi,0), 2);
1204 SLV_BPARM_MACRO(KPRINT_PTR,parameters);
1205 #endif /* __WIN32__ */
1206
1207 #ifndef __WIN32__ /* pc fortran cannot do prints to stdout (file 6) */
1208 slv_define_parm(parameters, bool_parm,
1209 "IDEBUG", "print debug (IDEBUG)", "print debug",
1210 U_p_bool(val,idefault[2]),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1211 SLV_BPARM_MACRO(IDEBUG_PTR,parameters);
1212 #else
1213 /* force FALSE value */
1214 slv_define_parm(parameters, bool_parm,
1215 "IDEBUG", "print debug (IDEBUG)", "print debug",
1216 U_p_bool(val,0),U_p_bool(lo,0),U_p_bool(hi,0), 2);
1217 SLV_BPARM_MACRO(IDEBUG_PTR,parameters);
1218 #endif /* __WIN32__ */
1219
1220 slv_define_parm(parameters, int_parm,
1221 "IOPTION", "convergence option (IOPTION)",
1222 "1-linesearch, 2-trustregion, 3-both",
1223 U_p_int(val,idefault[3]),U_p_int(lo,1),U_p_int(hi,3),1);
1224 SLV_IPARM_MACRO(IOPTION_PTR,parameters);
1225
1226 slv_define_parm(parameters, int_parm,
1227 "ICHOOSE", "basis selection (ICHOOSE)",
1228 "selection: 0-numeric, 1-first columns, 2-last columns",
1229 U_p_int(val,idefault[4]),U_p_int(lo,0),U_p_int(hi,2),1);
1230 SLV_IPARM_MACRO(ICHOOSE_PTR,parameters);
1231
1232 #ifndef __WIN32__
1233 slv_define_parm(parameters, int_parm,
1234 "ICORR", "correction type (ICORR)",
1235 "correction: 0-none, 1-FD with mult., 2-FD only, " /* no comma */
1236 "5-Interactive (FD/none)",
1237 U_p_int(val,idefault[5]),U_p_int(lo,0),U_p_int(hi,5),2);
1238 SLV_IPARM_MACRO(ICORR_PTR,parameters);
1239 #else
1240 /* no interactive input on pcs */
1241 slv_define_parm(parameters, int_parm,
1242 "ICORR", "correction type (ICORR)",
1243 "correction: 0-none, 1-Finite Difference with multiplier,"
1244 " 2-Finite Difference only",
1245 U_p_int(val,idefault[5]),U_p_int(lo,0),U_p_int(hi,2),2);
1246 SLV_IPARM_MACRO(ICORR_PTR,parameters);
1247 #endif
1248
1249 slv_define_parm(parameters, bool_parm,
1250 "I_MULT_FREE", "multiplier option (I_MULT_FREE)",
1251 "0-calculate multipliers, 1-multiplier free",
1252 U_p_bool(val,idefault[6]),U_p_bool(lo,0),U_p_bool(hi,1),2);
1253 SLV_BPARM_MACRO(I_MULT_FREE_PTR,parameters);
1254
1255 slv_define_parm(parameters, int_parm,
1256 "IIEXACT", "merit function for linesearch (IIEXACT)",
1257 "0-none, 1-exact, 2-Watchdog",
1258 U_p_int(val, idefault[7]),U_p_int(lo, 0),U_p_int(hi,2),1);
1259 SLV_IPARM_MACRO(IIEXACT_PTR,parameters);
1260
1261 slv_define_parm(parameters, bool_parm,
1262 "IPMETHOD", "Interior point QP solver (IPMETHOD)",
1263 "0-use active-set method, 1-use interior-point method",
1264 U_p_bool(val,idefault[8]),U_p_bool(lo,0),U_p_bool(hi,1), 1);
1265 SLV_BPARM_MACRO(IPMETHOD_PTR,parameters);
1266
1267
1268 slv_define_parm(parameters, char_parm,
1269 "ILUSOLVE", "Barrier method type for IP method (ILUSOLVE)",
1270 "combinations of many routines",
1271 U_p_string(val,ilu_names[ilu_name_from_int(idefault[9])]),
1272 U_p_strings(lo,ilu_names),
1273 U_p_int(hi,sizeof(ilu_names)/sizeof(char *)),1);
1274 SLV_CPARM_MACRO(ILUSOLVE_PTR,parameters);
1275
1276 slv_define_parm(parameters, bool_parm,
1277 "IWARM", "warm start for QP solver (IWARM)",
1278 "0-don't use warm start strategy, 1-use warm start strategy",
1279 U_p_bool(val,idefault[10]),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1280 SLV_BPARM_MACRO(IWARM_PTR,parameters);
1281
1282 slv_define_parm(parameters, real_parm,
1283 "IVV", "initial trust region multiplier value (IVV)",
1284 "0.0-use default value, value-initial multiplier value",
1285 U_p_real(val,0.0),U_p_real(lo,0.0),U_p_real(hi,100000),2);
1286 SLV_RPARM_MACRO(IVV_PTR,parameters);
1287
1288 slv_define_parm(parameters, bool_parm,
1289 "IDD", "input trust region shape (IDD)",
1290 "0-use default shape, 1-input trust region shape",
1291 U_p_bool(val,0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1292 SLV_BPARM_MACRO(IDD_PTR,parameters);
1293
1294 #ifndef __WIN32__ /* pc fortran cannot do prints to stdout (file 6) */
1295 slv_define_parm(parameters, bool_parm,
1296 "CNR", "file output option (CNR)",
1297 "0-no output to file, 1-print output to journal.num",
1298 U_p_bool(val,idefault[11]),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1299 SLV_BPARM_MACRO(CNR_PTR,parameters);
1300 #else
1301 /* force FALSE value */
1302 slv_define_parm(parameters, bool_parm,
1303 "CNR", "file output option (CNR)",
1304 "0-no output to file, 1-print output to journal.num",
1305 U_p_bool(val,0),U_p_bool(lo,0),U_p_bool(hi,0), 2);
1306 SLV_BPARM_MACRO(CNR_PTR,parameters);
1307 #endif /* __WIN32__ */
1308
1309 slv_define_parm(parameters, real_parm,
1310 "TIMEMAX", "maximum CPU time",
1311 "most CPU seconds allowed without returning to user",
1312 U_p_real(val,180),U_p_real(lo, 1),U_p_real(hi,1800), 1);
1313 SLV_RPARM_MACRO(TIMEMAX_PTR,parameters);
1314
1315 slv_define_parm(parameters, int_parm,
1316 "cutoff", "block size cutoff (MODEL-based)",
1317 "block size cutoff ( MODEL-based)",
1318 U_p_int(val, 500),U_p_int(lo,0),U_p_int(hi,20000), 2);
1319 SLV_IPARM_MACRO(CUTOFF_PTR,parameters);
1320
1321 if (parameters->num_parms != slv2_PA_SIZE) {
1322 FPRINTF(stderr,"%s parameter count (%d) != expected (%d) in %s.\n",
1323 NAME,parameters->num_parms,slv2_PA_SIZE,__FILE__);
1324 }
1325
1326 pop_controls(icontrol,rcontrol);
1327 return 0;
1328 }
1329
1330 /*********************** END lbopt interface ******************************/
1331
1332
1333 /* pretty close. fixme */
1334 static SlvClientToken slv2_create(slv_system_t server, int *statusindex)
1335 {
1336 slv2_system_t sys;
1337
1338 sys = (slv2_system_t)ascmalloc( sizeof(struct slv2_system_structure) );
1339 if (sys==NULL) {
1340 *statusindex = 1;
1341 return sys;
1342 }
1343 SERVER = server;
1344 /* init scalars for slv2 */
1345 sys->integrity = OK;
1346 sys->reorderstyle = s2_natural;
1347
1348 sys->p.parms = sys->padata;
1349 sys->p.dynamic_parms = 0;
1350 sys->p.output.more_important = stdout;
1351 sys->p.output.less_important = stdout;
1352 sys->p.tolerance.feasible = (double)1e-6; /* ? what is rsqp param? */
1353 sys->p.whose = (*statusindex);
1354
1355 sys->s.ok = TRUE;
1356 sys->s.calc_ok = TRUE;
1357 sys->s.costsize = 0;
1358 sys->s.cost = NULL;
1359
1360 /* fetch MODEL server info */
1361
1362 sys->rlist = slv_get_solvers_rel_list(SERVER);
1363 sys->vlist = slv_get_solvers_var_list(SERVER);
1364 sys->obj = slv_get_obj_relation(SERVER);
1365
1366 if (sys->vlist == NULL) {
1367 ascfree(sys);
1368 FPRINTF(stderr,"%s called with no variables.\n",NAME);
1369 *statusindex = -2;
1370 return NULL;
1371 }
1372 if (sys->rlist == NULL && sys->obj == NULL) {
1373 ascfree(sys);
1374 FPRINTF(stderr,"%s called with no relations or objective.\n",NAME);
1375 *statusindex = -1;
1376 return NULL;
1377 }
1378
1379 /* run server checks */
1380 slv_check_var_initialization(SERVER);
1381
1382 /* get control defaults from fortran */
1383 slv2_get_default_parameters(SERVER,(SlvClientToken)sys,&(sys->p));
1384
1385 sys->rfilter = &(sys->relfilter);
1386 sys->vfilter = &(sys->varfilter);
1387 sys->relfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1388 sys->relfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1389 sys->varfilter.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
1390 sys->varfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
1391
1392 /* init to NULL stuff to be used in presolve, etc */
1393 sys->acon = NULL;
1394 sys->avar = NULL;
1395 sys->gxsparse = NULL;
1396 sys->ivec = NULL;
1397 sys->jvec = NULL;
1398 sys->iw = NULL;
1399
1400 sys->a = NULL;
1401 sys->bnds = NULL;
1402 sys->c = NULL;
1403 sys->g = NULL;
1404 sys->gsparse = NULL;
1405 sys->nominals = NULL;
1406 sys->xn = NULL;
1407 sys->trust = NULL;
1408 sys->rw = NULL;
1409 sys->weights = NULL;
1410
1411 /* do this in presolve?? */
1412 sys->iwsize = 10000; /* ??? */
1413 sys->rwsize = 20000;
1414
1415 *statusindex = 0;
1416 return (SlvClientToken)sys;
1417 }
1418
1419
1420 /*
1421 *
1422 */
1423 static
1424 int slv2_destroy(slv_system_t server, SlvClientToken asys)
1425 {
1426 slv2_system_t sys;
1427 (void) server;
1428 sys = SLV2(asys);
1429 if (check_system(sys)) return 1;
1430
1431 slv_destroy_parms(&(sys->p));
1432
1433 /* free int lists */
1434 free_vector(sys->acon );
1435 free_vector(sys->avar);
1436 free_vector(sys->gxsparse);
1437 free_vector(sys->ivec);
1438 free_vector(sys->jvec);
1439 free_vector(sys->iw);
1440
1441 sys->acon = NULL;
1442 sys->avar = NULL;
1443 sys->gxsparse = NULL;
1444 sys->ivec = NULL;
1445 sys->jvec = NULL;
1446 sys->iw = NULL;
1447
1448 /* free double lists */
1449 free_vector(sys->a);
1450 free_vector(sys->bnds);
1451 free_vector(sys->c);
1452 free_vector(sys->g);
1453 free_vector(sys->gsparse);
1454 free_vector(sys->nominals);
1455 free_vector(sys->xn);
1456 free_vector(sys->trust);
1457 free_vector(sys->rw);
1458 free_vector(sys->weights);
1459
1460 sys->a = NULL;
1461 sys->bnds = NULL;
1462 sys->c = NULL;
1463 sys->g = NULL;
1464 sys->gsparse = NULL;
1465 sys->nominals = NULL;
1466 sys->xn = NULL;
1467 sys->trust = NULL;
1468 sys->rw = NULL;
1469 sys->weights = NULL;
1470
1471 sys->integrity = DESTROYED;
1472 ascfree( (POINTER)sys );
1473 return 0;
1474 }
1475
1476 /*************************************************************************
1477 * problem setup functions.
1478 ************************************************************************/
1479
1480 /* clear inblock flags after reordering so next block is happy */
1481 static void clear_inblock(slv2_system_t sys, const mtx_block_t *b, int32 bnum)
1482 {
1483 mtx_region_t reg;
1484 int32 row,col;
1485
1486 reg = b->block[bnum];
1487 for( col = reg.col.low; col <= reg.col.high; col++ ) {
1488 var_set_in_block(sys->vlist[col],FALSE);
1489 }
1490 for( row = reg.row.low; row <= reg.row.high; row++ ) {
1491 rel_set_in_block(sys->rlist[row],FALSE);
1492 }
1493 }
1494
1495 /* reset the block in question to square. since we
1496 * do *nothing* with blocks after the calling function, we don't care.
1497 */
1498 static void square_block(const mtx_block_t *b, int32 bnum)
1499 {
1500 if (b->block[bnum].col.high > b->block[bnum].row.high) {
1501 b->block[bnum].col.high = b->block[bnum].row.high;
1502 return;
1503 }
1504 if (b->block[bnum].row.high > b->block[bnum].col.high) {
1505 b->block[bnum].row.high = b->block[bnum].col.high;
1506 return;
1507 }
1508 }
1509
1510 /* move equalities to rlist [0..rlen-1] and free vars to 0..vlen-1]
1511 * and calculate rlen, vlen.
1512 * Inequalities, if any, are dropped.
1513 * Opt has provisions that allow us to munge inequality to equality
1514 * if desired in future.
1515 * Preconditioning ma28 by reorderings can be done here.
1516 * potentially lots of wasted effort here.
1517 */
1518 static
1519 int slv2_rearrange_lists(slv2_system_t sys)
1520 {
1521 int32 bnum,i;
1522 const mtx_block_t *b;
1523
1524 /* calculate sys->reorderstyle from future string parameter here */
1525
1526 slv_sort_rels_and_vars(SERVER,&(sys->rlen),&(sys->vlen));
1527 switch (sys->reorderstyle) {
1528 case s2_tspk1:
1529 slv_block_partition(SERVER); /* prolly --> rectangle block */
1530 b = slv_get_solvers_blocks(SERVER);
1531 if (b != NULL && b->nblocks > 0) {
1532 bnum = b->nblocks;
1533 for (i=0; i < bnum; i++) {
1534 square_block(b,i);
1535 slv_spk1_reorder_block(SERVER,i,1);
1536 clear_inblock(sys,b,i);
1537 }
1538 }
1539 break;
1540 case s2_spk1:
1541 slv_block_partition(SERVER); /* prolly --> rectangle block */
1542 b = slv_get_solvers_blocks(SERVER);
1543 if (b != NULL && b->nblocks > 0) {
1544 bnum = b->nblocks;
1545 for (i=0; i < bnum; i++) {
1546 square_block(b,i);
1547 slv_spk1_reorder_block(SERVER,i,0);
1548 clear_inblock(sys,b,i);
1549 }
1550 }
1551 break;
1552 case s2_tteardrop:
1553 slv_block_partition(SERVER); /* prolly --> rectangle block */
1554 b = slv_get_solvers_blocks(SERVER);
1555 if (b != NULL && b->nblocks > 0) {
1556 bnum = b->nblocks;
1557 for (i=0; i < bnum; i++) {
1558 slv_tear_drop_reorder_block(SERVER,i,CUTOFF,0,mtx_TSPK1);
1559 clear_inblock(sys,b,i);
1560 }
1561 }
1562 break;
1563 case s2_teardrop:
1564 slv_block_partition(SERVER); /* prolly --> rectangle block */
1565 b = slv_get_solvers_blocks(SERVER);
1566 if (b != NULL && b->nblocks > 0) {
1567 bnum = b->nblocks;
1568 for (i=0; i < bnum; i++) {
1569 slv_tear_drop_reorder_block(SERVER,i,CUTOFF,0,mtx_SPK1);
1570 clear_inblock(sys,b,i);
1571 }
1572 }
1573 break;
1574 case s2_tteardropfat:
1575 slv_block_partition(SERVER); /* prolly --> rectangle block */
1576 b = slv_get_solvers_blocks(SERVER);
1577 if (b != NULL && b->nblocks > 0) {
1578 bnum = b->nblocks;
1579 for (i=0; i < bnum; i++) {
1580 slv_tear_drop_reorder_block(SERVER,i,CUTOFF,1,mtx_TSPK1);
1581 clear_inblock(sys,b,i);
1582 }
1583 }
1584 break;
1585 case s2_teardropfat:
1586 slv_block_partition(SERVER); /* prolly --> rectangle block */
1587 b = slv_get_solvers_blocks(SERVER);
1588 if (b != NULL && b->nblocks > 0) {
1589 bnum = b->nblocks;
1590 for (i=0; i < bnum; i++) {
1591 slv_tear_drop_reorder_block(SERVER,i,CUTOFF,1,mtx_SPK1);
1592 clear_inblock(sys,b,i);
1593 }
1594 }
1595 break;
1596 case s2_random: /* not yet implemented */
1597 case s2_numeric: /* not yet implemented */
1598 case s2_natural: /* natural is the default */
1599 default:
1600 break;
1601 }
1602 return 0;
1603 }
1604
1605 /*
1606 * Presolve allocates all the needed workspaces and initializes
1607 * control parameters.
1608 PRESO
1609 */
1610 static
1611 void slv2_presolve(slv_system_t server, SlvClientToken asys)
1612 {
1613 int32 big;
1614 struct var_variable **vp;
1615 struct rel_relation **rp;
1616 int32 col, row, ind, dummy;
1617 real64 value;
1618
1619 slv2_system_t sys;
1620 (void) server;
1621 sys = SLV2(asys);
1622
1623 big = INT_MAX/16;
1624
1625 sys->s.ready_to_solve = FALSE;
1626 sys->nfunc = 0;
1627 sys->ngrad = 0;
1628
1629 iteration_begins(sys);
1630 check_system(sys);
1631 if( sys->vlist == NULL ) {
1632 FPRINTF(stderr,"ERROR: (slv2) slv2_presolve\n");
1633 FPRINTF(stderr," Variable list was never set.\n");
1634 return;
1635 }
1636
1637 if( sys->rlist == NULL ) {
1638 FPRINTF(stderr,"ERROR: (slv2) slv2_presolve\n");
1639 FPRINTF(stderr," Relation list was never set.\n");
1640 return;
1641 }
1642
1643 /* Reset global iteration count and cpu-elapsed. */
1644 sys->s.iteration = 0L;
1645 sys->s.cpu_elapsed = D_ZERO;
1646
1647 /* Do a variety of things according to options to improve ma28/opt.
1648 * This leaves rlen and vlen set to what opt needs to know about.
1649 * If there are no constraints, does nothing.
1650 */
1651 sys->rlen = slv_get_num_solvers_rels(SERVER);
1652 if (sys->rlen) {
1653 slv2_rearrange_lists(sys);
1654 }
1655
1656 sys->s.block.current_reordered_block = -2;
1657
1658 sys->md = sys->rlen;
1659 sys->ibnds = 1; /* bounded variables problem */
1660
1661
1662 /* count nonzeros */
1663 sys->njac = relman_jacobian_count(sys->rlist,sys->rlen,
1664 sys->vfilter,sys->rfilter,
1665 &dummy);
1666
1667 /* create avar array */
1668 free_vector(sys->avar);
1669 sys->avar = alloc_ivector(sys->njac+1);
1670 free_vector(sys->jvec);
1671 sys->jvec = alloc_ivector(sys->njac+1);
1672
1673 /* create acon array */
1674 free_vector(sys->acon);
1675 sys->acon = alloc_ivector(sys->njac+1);
1676 free_vector(sys->ivec);
1677 sys->ivec = alloc_ivector(sys->njac+1);
1678
1679
1680 /* create objective gradient space */
1681 free_vector(sys->g);
1682 free_vector(sys->gsparse);
1683 free_vector(sys->gxsparse);
1684 sys->g = alloc_vector(sys->vlen+1);
1685 if (sys->obj != NULL) {
1686 sys->gsparse = alloc_vector(rel_n_incidences(sys->obj));
1687 sys->gxsparse = alloc_ivector(rel_n_incidences(sys->obj));
1688 } else {
1689 sys->gsparse = NULL;
1690 sys->gxsparse = NULL;
1691 }
1692
1693 /* create residual space */
1694 free_vector(sys->c);
1695 sys->c = alloc_vector(sys->rlen);
1696
1697 free_vector(sys->weights);
1698 sys->weights = alloc_vector(sys->rlen);
1699
1700 sys->nzd = MAX(1, sys->njac);
1701
1702
1703 CHECK_MEM(&(sys->inform),&big,&big,&(sys->vlen),&(sys->rlen),&(sys->nzd),
1704 &(sys->iwsize),&(sys->rwsize));
1705 free_vector(sys->iw); /* create opt integer workspace */
1706 sys->iw = alloc_ivector(sys->iwsize);
1707 if (!sys->iw) {
1708 FPRINTF(stderr,"Error mallocing opt integer workspace array!");
1709 return;
1710 }
1711
1712 free_vector(sys->rw); /* create opt real workspace */
1713 sys->rw = alloc_vector(sys->rwsize);
1714 if ( !sys->rw) {
1715 FPRINTF(stderr,"Error mallocing opt real workspace array!");
1716 return;
1717 }
1718
1719 free_vector(sys->bnds); /* create bounds array */
1720 sys->bnds = alloc_vector(2*(sys->vlen+1));
1721 if (!sys->bnds) {
1722 FPRINTF(stderr,"Error mallocing bnds array!");
1723 return;
1724 }
1725
1726 free_vector(sys->xn); /* create solution array */
1727 sys->xn = alloc_vector(sys->vlen +1);
1728
1729 free_vector(sys->nominals); /* create xn scaling array */
1730 sys->nominals = alloc_vector(sys->vlen);
1731
1732 free_vector(sys->trust); /* create trust region step limit array */
1733 sys->trust = alloc_vector(sys->vlen +1);
1734
1735 free_vector(sys->a); /* create jacobian element array */
1736 sys->a = alloc_vector(sys->njac+1);
1737 if (!sys->xn) {
1738 FPRINTF(stderr,"Error mallocing variable array!");
1739 return;
1740 }
1741 if (!sys->a) {
1742 FPRINTF(stderr,"Error mallocing jacobian array!");
1743 return;
1744 }
1745
1746 if (scale_and_bounds(sys) > 0) {
1747 FPRINTF(MIF(sys),"Unable to presolve; bad initial point.\n");
1748 sys->s.ready_to_solve = FALSE;
1749 return;
1750 }
1751
1752 sys->iluoption = ilu_name_to_int(sys);
1753
1754 sys->s.over_defined =
1755 sys->s.under_defined =
1756 sys->s.struct_singular =
1757 sys->s.converged =
1758 sys->s.diverged =
1759 sys->s.inconsistent = FALSE;
1760 sys->s.block.number_of = 1;
1761
1762 sys->s.block.previous_total_size = 0;
1763 sys->s.block.current_block = 0;
1764 sys->s.block.current_size = sys->maxndx;
1765 sys->s.block.iteration = 0L;
1766 sys->s.block.cpu_elapsed = D_ZERO;
1767 sys->s.calc_ok = TRUE;
1768 sys->s.ready_to_solve = TRUE;
1769
1770 sys->inform = -1;
1771
1772 iteration_ends(sys);
1773 }
1774
1775 /***************************************************************************
1776 * Below here is plain vanilla ascendism
1777 */
1778
1779 static void
1780 slv2_dump_internals(slv_system_t server, SlvClientToken sys,int level)
1781 {
1782 /* what do we want to put here? */
1783 }
1784
1785 static
1786 int slv2_eligible_solver(slv_system_t server)
1787 {
1788 return 1; /* fixme */
1789 }
1790
1791 static
1792 void slv2_get_parameters(slv_system_t server, SlvClientToken asys,
1793 slv_parameters_t *parameters)
1794 {
1795 slv2_system_t sys;
1796 (void) server;
1797 sys = SLV2(asys);
1798 if (check_system(sys)) return;
1799 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
1800 }
1801
1802 static void slv2_set_parameters(slv_system_t server, SlvClientToken asys,
1803 slv_parameters_t *parameters)
1804 {
1805 slv2_system_t sys;
1806 (void) server;
1807 sys = SLV2(asys);
1808 if (check_system(sys)) return;
1809 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
1810 }
1811
1812 static void slv2_get_status(slv_system_t server, SlvClientToken asys,
1813 slv_status_t *status)
1814 {
1815 slv2_system_t sys;
1816 (void) server;
1817 sys = SLV2(asys);
1818 if (check_system(sys)) return;
1819 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
1820 }
1821
1822 static
1823 void real_iterate(slv_system_t server, SlvClientToken asys,int needpush)
1824 {
1825 slv2_system_t sys;
1826 sys = SLV2(asys);
1827 check_system(sys);
1828 if( !sys->s.ready_to_solve ) {
1829 FPRINTF(stderr,"ERROR: (slv2) slv2_iterate\n");
1830 FPRINTF(stderr," Not ready to solve.\n");
1831 return;
1832 }
1833
1834 if (needpush) {
1835 push_controls(sys->icontrol, sys->rcontrol);
1836 }
1837 invoke_opt(sys);
1838 if (needpush) {
1839 pop_controls(sys->icontrol, sys->rcontrol);
1840 }
1841
1842 }
1843
1844 static
1845 void slv2_iterate(slv_system_t server, SlvClientToken asys)
1846 {
1847 real_iterate(server, asys, 1);
1848 }
1849
1850 static void slv2_solve(slv_system_t server, SlvClientToken asys)
1851 {
1852 slv2_system_t sys;
1853 sys = SLV2(asys);
1854 if (server == NULL || sys==NULL) return;
1855 if (check_system(sys)) return;
1856 push_controls(sys->icontrol, sys->rcontrol);
1857 while( sys->s.ready_to_solve ) {
1858 real_iterate(server,asys,0);
1859 }
1860 pop_controls(sys->icontrol, sys->rcontrol);
1861 }
1862
1863 int slv2_register(SlvFunctionsT *sft)
1864 {
1865 if (sft==NULL) {
1866 FPRINTF(stderr,"slv2_register called with NULL pointer\n");
1867 return 1;
1868 }
1869
1870 sft->name = NAME;
1871 sft->ccreate = slv2_create;
1872 sft->cdestroy = slv2_destroy;
1873 sft->celigible = slv2_eligible_solver;
1874 sft->getdefparam = slv2_get_default_parameters;
1875 sft->getparam = slv2_get_parameters;
1876 sft->setparam = slv2_set_parameters;
1877 sft->getstatus = slv2_get_status;
1878 sft->solve = slv2_solve;
1879 sft->presolve = slv2_presolve;
1880 sft->iterate = slv2_iterate;
1881 sft->resolve = NULL;
1882 sft->getlinsol = NULL;
1883 sft->getlinsys = NULL;
1884 sft->getsysmtx = NULL;
1885 sft->dumpinternals = slv2_dump_internals;
1886 /* it should be safe, but unnecessary, to call this at every iteration
1887 * according to Dave.
1888 */
1889 INIT_OPT_COMMON();
1890 return 0;
1891 }
1892
1893 #endif /* #else clause of DYNAMIC_OPTSQP */
1894 #endif /* #else clause of !STATIC_OPTSQP && !DYNAMIC_OPTSQP */

Properties

Name Value
svn:executable *

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