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

Contents of /trunk/ascend4/solver/slv8.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: 87188 byte(s)
Setting up web subdirectory in repository
1 /*
2 * Incorporation of the nonlinear solver CONOPT to ASCEND
3 * by Ken Tyner
4 * Created: 6/97
5 * Version: $Revision: 1.31 $
6 * Version control file: $RCSfile: slv8.c,v $
7 * Date last modified: $Date: 2000/01/25 02:27:50 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1997 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 "utilities/ascConfig.h"
33 #include "utilities/ascMalloc.h"
34 #include "utilities/set.h"
35 #include "general/tm_time.h"
36 #include "utilities/mem.h"
37 #include "general/list.h"
38 #include "compiler/fractions.h"
39 #include "compiler/dimen.h"
40 #include "compiler/functype.h"
41 #include "compiler/func.h"
42 #include "solver/mtx.h"
43 #include "solver/linsol.h"
44 #include "solver/linsolqr.h"
45 #include "solver/slv_types.h"
46 #include "solver/var.h"
47 #include "solver/rel.h"
48 #include "solver/discrete.h"
49 #include "solver/conditional.h"
50 #include "solver/logrel.h"
51 #include "solver/bnd.h"
52 #include "solver/calc.h"
53 #include "solver/relman.h"
54 #include "solver/slv_common.h"
55 #include "solver/slv_client.h"
56 #include "solver/conopt.h"
57 #include "solver/slv8.h"
58 #include "solver/slv_stdcalls.h"
59 #include "solver/conoptdll.h"
60
61 #define slv8_register_conopt_function register_conopt_function
62 #define slv8_coicsm coicsm
63 #define slv8_coimem coimem
64
65 #if !defined(STATIC_CONOPT) && !defined(DYNAMIC_CONOPT)
66
67 int slv8_register(SlvFunctionsT *f)
68 {
69 (void)f; /* stop gcc whine about unused parameter */
70
71 FPRINTF(stderr,"CONOPT not compiled in this ASCEND IV.\n");
72 return 1;
73 }
74
75 #else /* either STATIC_CONOPT or DYNAMIC_CONOPT is defined */
76
77 /* #ifdef DYNAMIC_CONOPT */
78 /* do dynamic loading stuff. yeah, right */
79
80 /* #else *//* following is used if STATIC_CONOPT is defined */
81
82 /*
83 * Output in user defined CONOPT subroutines
84 */
85 #define CONDBG 0
86 #define NONBASIC_DEBUG FALSE
87
88 /*
89 * makes lots of extra spew
90 */
91 #define DEBUG FALSE
92
93 #define SLV8(s) ((slv8_system_t)(s))
94 #define MI8F(s) slv_get_output_file( SLV8(s)->p.output.more_important )
95 #define SERVER (sys->slv)
96 #define slv8_PA_SIZE 33
97 #define SAFE_CALC_PTR (sys->parm_array[0])
98 #define SAFE_CALC ((*(int *)SAFE_CALC_PTR))
99 #define SCALEOPT_PTR (sys->parm_array[1])
100 #define SCALEOPT ((*(char **)SCALEOPT_PTR))
101 #define TOO_SMALL_PTR (sys->parm_array[2])
102 #define TOO_SMALL ((*(real64 *)TOO_SMALL_PTR))
103 #define UPDATE_NOMINALS_PTR (sys->parm_array[3])
104 #define UPDATE_NOMINALS ((*(int *)UPDATE_NOMINALS_PTR))
105 #define UPDATE_RELNOMS_PTR (sys->parm_array[4])
106 #define UPDATE_RELNOMS ((*(int *)UPDATE_RELNOMS_PTR))
107 #define UPDATE_WEIGHTS_PTR (sys->parm_array[5])
108 #define UPDATE_WEIGHTS ((*(int *)UPDATE_WEIGHTS_PTR))
109 #define DUMPCNORM_PTR (sys->parm_array[6])
110 #define DUMPCNORM ((*(int *)DUMPCNORM_PTR))
111 #define CNLOW_PTR (sys->parm_array[7])
112 #define CNLOW ((*(real64 *)CNLOW_PTR))
113 #define CNHIGH_PTR (sys->parm_array[8])
114 #define CNHIGH ((*(real64 *)CNHIGH_PTR))
115 #define UPDATE_JACOBIAN_PTR (sys->parm_array[9])
116 #define UPDATE_JACOBIAN ((*(int *)UPDATE_JACOBIAN_PTR))
117 #define ITSCALELIM_PTR (sys->parm_array[10])
118 #define ITSCALELIM ((*(int *)ITSCALELIM_PTR))
119 #define ITSCALETOL_PTR (sys->parm_array[11])
120 #define ITSCALETOL ((*(real64 *)ITSCALETOL_PTR))
121 #define LIFDS_PTR (sys->parm_array[12])
122 #define LIFDS ((*(int32 *)LIFDS_PTR))
123 #define REORDER_OPTION_PTR (sys->parm_array[13])
124 #define REORDER_OPTION ((*(char **)REORDER_OPTION_PTR))
125 #define CUTOFF_PTR (sys->parm_array[14])
126 #define CUTOFF ((*(int32 *)CUTOFF_PTR))
127 #define RELNOMSCALE_PTR (sys->parm_array[15])
128 #define RELNOMSCALE ((*(int *)RELNOMSCALE_PTR))
129 #define ITER_LIMIT_PTR (sys->parm_array[16])
130 #define ITER_LIMIT ((*(int32 *)ITER_LIMIT_PTR))
131 #define TIME_LIMIT_PTR (sys->parm_array[17])
132 #define TIME_LIMIT ((*(int32 *)TIME_LIMIT_PTR))
133 #define DOMLIM_PTR (sys->parm_array[18])
134 #define DOMLIM ((*(int32 *)DOMLIM_PTR))
135 #define RTMAXJ_PTR (sys->parm_array[19])
136 #define RTMAXJ ((*(real64 *)RTMAXJ_PTR))
137
138 /*
139 * Auxiliar structures
140 */
141
142 struct update_data {
143 int jacobian; /* Countdown on jacobian updating */
144 int weights; /* Countdown on weights updating */
145 int nominals; /* Countdown on nominals updating */
146 int relnoms; /* Countdown on relnom updating */
147 int iterative; /* Countdown on iterative scale update */
148 };
149
150
151 /*
152 * varpivots, relpivots used only in optimizing, if we rewrite calc_pivots
153 * without them.
154 */
155 struct jacobian_data {
156 linsolqr_system_t sys; /* Linear system */
157 mtx_matrix_t mtx; /* Transpose gradient of residuals */
158 real64 *rhs; /* RHS from linear system */
159 unsigned *varpivots; /* Pivoted variables */
160 unsigned *relpivots; /* Pivoted relations */
161 unsigned *subregions; /* Set of subregion indeces */
162 dof_t *dofdata; /* dof data pointer from server */
163 mtx_region_t reg; /* Current block region */
164 int32 rank; /* Numerical rank of the jacobian */
165 enum factor_method fm; /* Linear factorization method */
166 boolean accurate; /* ? Recalculate matrix */
167 boolean singular; /* ? Can matrix be inverted */
168 boolean old_partition; /* old value of partition flag */
169 };
170
171 struct slv8_system_structure {
172
173 /*
174 * Problem definition
175 */
176 slv_system_t slv; /* slv_system_t back-link */
177 struct rel_relation *obj; /* Objective function: NULL = none */
178 struct rel_relation *old_obj;/* Objective function: NULL = none */
179 struct var_variable **vlist; /* Variable list (NULL terminated) */
180 struct rel_relation **rlist; /* Relation list (NULL terminated) */
181
182 /*
183 * Solver information
184 */
185 int integrity; /* ? Has the system been created */
186 int32 presolved; /* ? Has the system been presolved */
187 int32 resolve; /* ? Has the system been resolved */
188 slv_parameters_t p; /* Parameters */
189 slv_status_t s; /* Status (as of iteration end) */
190 struct update_data update; /* Jacobian frequency counters */
191 int32 cap; /* Order of matrix/vectors */
192 int32 rank; /* Symbolic rank of problem */
193 int32 vused; /* Free and incident variables */
194 int32 vtot; /* length of varlist */
195 int32 rused; /* Included relations */
196 int32 rtot; /* length of rellist */
197 double clock; /* CPU time */
198
199 void *parm_array[slv8_PA_SIZE];
200 struct slv_parameter pa[slv8_PA_SIZE];
201
202 /*
203 * CONOPT DATA
204 */
205 struct conopt_data con;
206
207 /*
208 * Calculated data (scaled)
209 */
210 struct jacobian_data J; /* linearized system */
211
212 struct vector_data nominals; /* Variable nominals */
213 struct vector_data weights; /* Relation weights */
214 struct vector_data relnoms; /* Relation nominals */
215 struct vector_data variables; /* Variable values */
216 struct vector_data residuals; /* Relation residuals */
217
218 real64 objective; /* Objective function evaluation */
219 };
220
221
222 /*
223 * Integrity checks
224 * ----------------
225 * check_system(sys)
226 */
227
228 #define OK ((int32)813029392)
229 #define DESTROYED ((int32)103289182)
230
231 /*
232 * Checks sys for NULL and for integrity.
233 */
234 static int check_system(slv8_system_t sys)
235 {
236 if( sys == NULL ) {
237 FPRINTF(stderr,"ERROR: (slv8) check_system\n");
238 FPRINTF(stderr," NULL system handle.\n");
239 return 1;
240 }
241
242 switch( sys->integrity ) {
243 case OK:
244 return 0;
245 case DESTROYED:
246 FPRINTF(stderr,"ERROR: (slv8) check_system\n");
247 FPRINTF(stderr," System was recently destroyed.\n");
248 return 1;
249 default:
250 FPRINTF(stderr,"ERROR: (slv8) check_system\n");
251 FPRINTF(stderr," System reused or never allocated.\n");
252 return 1;
253 }
254 }
255
256 /*
257 * General input/output routines
258 * -----------------------------
259 * print_var_name(out,sys,var)
260 * print_rel_name(out,sys,rel)
261 */
262
263 #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
264 #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
265
266 /*
267 * Debug output routines
268 * ---------------------
269 * debug_delimiter(fp)
270 * debug_out_vector(fp,vec)
271 * debug_out_var_values(fp,sys)
272 * debug_out_rel_residuals(fp,sys)
273 * debug_out_jacobian(fp,sys)
274 * debug_out_hessian(fp,sys)
275 */
276
277 /*
278 * Outputs a hyphenated line.
279 */
280 static void debug_delimiter( FILE *fp)
281 {
282 int32 i;
283 for( i=0; i<60; i++ ) PUTC('-',fp);
284 PUTC('\n',fp);
285 }
286
287 #if DEBUG
288
289 /*
290 * Outputs a vector.
291 */
292 static void debug_out_vector( FILE *fp, slv8_system_t sys,
293 struct vector_data *vec)
294 {
295 int32 ndx;
296 FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n",
297 calc_sqrt_D0(vec->norm2), vec->accurate?"TRUE":"FALSE",
298 vec->rng->low,vec->rng->high);
299 FPRINTF(fp,"Vector --> ");
300 for( ndx=vec->rng->low ; ndx<=vec->rng->high ; ++ndx )
301 FPRINTF(fp, "%g ", vec->vec[ndx]);
302 PUTC('\n',fp);
303 }
304
305 /*
306 * Outputs all variable values in current block.
307 */
308 static void debug_out_var_values( FILE *fp, slv8_system_t sys)
309 {
310 int32 col;
311 struct var_variable *var;
312
313 FPRINTF(fp,"Var values --> \n");
314 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high ; col++ ) {
315 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
316 print_var_name(fp,sys,var);
317 FPRINTF(fp, "\nI Lb Value Ub Scale Col INom\n");
318 FPRINTF(fp,"%d\t%.4g\t%.4g\t%.4g\t%.4g\t%d\t%.4g\n",
319 var_sindex(var),var_lower_bound(var),var_value(var),
320 var_upper_bound(var),var_nominal(var),
321 col,sys->nominals.vec[col]);
322 }
323 }
324
325 /*
326 * Outputs all relation residuals in current block.
327 */
328 static void debug_out_rel_residuals( FILE *fp, slv8_system_t sys)
329 {
330 int32 row;
331
332 FPRINTF(fp,"Rel residuals --> \n");
333 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high ; row++ ) {
334 struct rel_relation *rel;
335 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
336 FPRINTF(fp," %g : ",rel_residual(rel));
337 print_rel_name(fp,sys,rel);
338 PUTC('\n',fp);
339 }
340 PUTC('\n',fp);
341 }
342
343
344 /*
345 * Outputs permutation and values of the nonzero elements in the
346 * the jacobian matrix.
347 */
348 static void debug_out_jacobian( FILE *fp, slv8_system_t sys)
349 {
350 mtx_coord_t nz;
351 real64 value;
352
353 nz.row = sys->J.reg.row.low;
354 for( ; nz.row <= sys->J.reg.row.high; ++(nz.row) ) {
355 FPRINTF(fp," Row %d (rel %d)\n", nz.row,
356 mtx_row_to_org(sys->J.mtx,nz.row));
357 nz.col = mtx_FIRST;
358 while( value = mtx_next_in_row(sys->J.mtx,&nz,&(sys->J.reg.col)),
359 nz.col != mtx_LAST ) {
360 FPRINTF(fp," Col %d (var %d) has value %g\n", nz.col,
361 mtx_col_to_org(sys->J.mtx,nz.col), value);
362 }
363 }
364 }
365
366 /*
367 * Outputs permutation and values of the nonzero elements in the
368 * reduced hessian matrix.
369 */
370 static void debug_out_hessian( FILE *fp, slv8_system_t sys)
371 {
372 mtx_coord_t nz;
373
374 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
375 nz.col = nz.row + sys->J.reg.col.high + 1 - sys->ZBZ.order;
376 FPRINTF(fp," ZBZ[%d (var %d)] = ",
377 nz.row, mtx_col_to_org(sys->J.mtx,nz.col));
378 for( nz.col = 0; nz.col < sys->ZBZ.order; nz.col++ ) {
379 FPRINTF(fp,"%10g ",sys->ZBZ.mtx[nz.row][nz.col]);
380 }
381 PUTC('\n',fp);
382 }
383 }
384
385 #endif /* DEBUG */
386
387
388 /*
389 * Array/vector operations
390 * ----------------------------
391 * destroy_array(p)
392 * create_array(len,type)
393 *
394 * zero_vector(vec)
395 * copy_vector(vec1,vec2)
396 * prod = inner_product(vec1,vec2)
397 * norm2 = square_norm(vec)
398 * matrix_product(mtx,vec,prod,scale,transpose)
399 */
400
401 #define destroy_array(p) \
402 if( (p) != NULL ) ascfree((p))
403 #define create_array(len,type) \
404 ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
405 #define create_zero_array(len,type) \
406 ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
407
408 #define zero_vector(v) slv_zero_vector(v)
409 #define copy_vector(v,t) slv_copy_vector((v),(t))
410 #define inner_product(v,u) slv_inner_product((v),(u))
411 #define square_norm(v) slv_square_norm(v)
412 #define matrix_product(m,v,p,s,t) slv_matrix_product((m),(v),(p),(s),(t))
413
414
415 /*
416 * Calculation routines
417 * --------------------
418 * ok = calc_objective(sys)
419 * ok = calc_residuals(sys)
420 * ok = calc_J(sys)
421 * calc_nominals(sys)
422 * calc_weights(sys)
423 * scale_J(sys)
424 * scale_variables(sys)
425 * scale_residuals(sys)
426 */
427
428 /*
429 * counts jacobian elements and sets max to the number of elements
430 * in the densest row
431 */
432 static int32 num_jacobian_nonzeros(slv8_system_t sys, int32 *max)
433 {
434 int32 row, len, licn,c,count,row_max;
435 struct rel_relation *rel;
436 rel_filter_t rf;
437 var_filter_t vf;
438 const struct var_variable **list;
439
440 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
441 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
442 vf.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
443 vf.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
444
445 licn = 0;
446 *max = 0;
447 row_max = sys->con.m;
448 if (sys->obj != NULL) {
449 row_max--;
450 }
451 /* replace at leisure with
452 * relman_jacobian_count(sys->rlist,row_max,&vfilter,&rfilter,max);
453 */
454 for( row = 0; row < row_max; row++ ) {
455 rel = sys->rlist[row];
456 if (rel_apply_filter(rel,&rf)) { /* shouldn't be needed */
457 len = rel_n_incidences(rel);
458 list = rel_incidence_list(rel);
459 count = 0;
460 for (c=0; c < len; c++) {
461 if( var_apply_filter(list[c],&vf) ) {
462 licn++;
463 count++;
464 }
465 }
466 *max = MAX(*max,count);
467 }
468 }
469 if (sys->obj != NULL) {
470 rel = sys->obj;
471 len = rel_n_incidences(rel);
472 list = rel_incidence_list(rel);
473 count = 0;
474 for (c=0; c < len; c++) {
475 if( var_apply_filter(list[c],&vf) ) {
476 licn++;
477 count++;
478 }
479 }
480 *max = MAX(*max,count);
481 }
482 return licn;
483 }
484
485
486 /*
487 * Evaluate the objective function.
488 */
489 static boolean calc_objective( slv8_system_t sys)
490 {
491 calc_ok = TRUE;
492 sys->objective =
493 (sys->obj ? relman_eval(sys->obj,&calc_ok,SAFE_CALC) : 0.0);
494 return calc_ok;
495 }
496
497
498 /*
499 * Evaluate all objectives.
500 */
501 static boolean calc_objectives( slv8_system_t sys)
502 {
503 int32 len,i;
504 static rel_filter_t rfilter;
505 struct rel_relation **rlist=NULL;
506 rfilter.matchbits = (REL_INCLUDED);
507 rfilter.matchvalue =(REL_INCLUDED);
508 rlist = slv_get_solvers_obj_list(SERVER);
509 len = slv_get_num_solvers_objs(SERVER);
510 calc_ok = TRUE;
511 for (i = 0; i < len; i++) {
512 if (rel_apply_filter(rlist[i],&rfilter)) {
513 relman_eval(rlist[i],&calc_ok,SAFE_CALC);
514 #if DEBUG
515 if (calc_ok == FALSE) {
516 FPRINTF(stderr,"error in calc_objectives\n");
517 calc_ok = TRUE;
518 }
519 #endif /* DEBUG */
520 }
521 }
522 return calc_ok;
523 }
524
525 /*
526 * Calculates all of the residuals in the current block and computes
527 * the residual norm for block status. Returns true iff calculations
528 * preceded without error.
529 */
530 static boolean calc_residuals( slv8_system_t sys)
531 {
532 int32 row;
533 struct rel_relation *rel;
534 double time0;
535
536 if( sys->residuals.accurate ) return TRUE;
537
538 calc_ok = TRUE;
539 row = sys->residuals.rng->low;
540 time0=tm_cpu_time();
541 for( ; row <= sys->residuals.rng->high; row++ ) {
542 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
543 #if DEBUG
544 if (!rel) {
545 int32r;
546 r=mtx_row_to_org(sys->J.mtx,row);
547 FPRINTF(MIF(sys),"NULL relation found !!\n");
548 FPRINTF(MIF(sys),"at row %d rel %d in calc_residuals\n",(int)row,r);
549 FFLUSH(MIF(sys));
550 }
551 #endif /* DEBUG */
552 sys->residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
553
554 relman_calc_satisfied(rel,sys->p.tolerance.feasible);
555 }
556 sys->s.block.functime += (tm_cpu_time() -time0);
557 sys->s.block.funcs++;
558 square_norm( &(sys->residuals) );
559 sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
560 return(calc_ok);
561 }
562
563
564 /*
565 * Calculates the current block of the jacobian.
566 * It is initially unscaled.
567 */
568 static boolean calc_J( slv8_system_t sys)
569 {
570 int32 row;
571 var_filter_t vfilter;
572 double time0;
573 real64 resid;
574
575 calc_ok = TRUE;
576 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
577 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
578 time0=tm_cpu_time();
579 mtx_clear_region(sys->J.mtx,&(sys->J.reg));
580 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
581 struct rel_relation *rel;
582 rel = sys->rlist[row];
583 relman_diffs(rel,&vfilter,sys->J.mtx,&resid,SAFE_CALC);
584 }
585 sys->s.block.jactime += (tm_cpu_time() - time0);
586 sys->s.block.jacs++;
587
588 if( --(sys->update.nominals) <= 0 ) sys->nominals.accurate = FALSE;
589 if( --(sys->update.weights) <= 0 ) sys->weights.accurate = FALSE;
590
591 return(calc_ok);
592 }
593
594 /*
595 * Retrieves the nominal values of all of the block variables,
596 * insuring that they are all strictly positive.
597 */
598 static void calc_nominals( slv8_system_t sys)
599 {
600 int32 col;
601 FILE *fp = MIF(sys);
602 if( sys->nominals.accurate ) return;
603 fp = MIF(sys);
604 col = sys->nominals.rng->low;
605 if(strcmp(SCALEOPT,"NONE") == 0 ||
606 strcmp(SCALEOPT,"ITERATIVE") == 0){
607 for( ; col <= sys->nominals.rng->high; col++ ) {
608 sys->nominals.vec[col] = 1;
609 }
610 } else {
611 for( ; col <= sys->nominals.rng->high; col++ ) {
612 struct var_variable *var;
613 real64 n;
614 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
615 n = var_nominal(var);
616 if( n <= 0.0 ) {
617 if( n == 0.0 ) {
618 n = TOO_SMALL;
619 FPRINTF(fp,"ERROR: (slv8) calc_nominals\n");
620 FPRINTF(fp," Variable ");
621 print_var_name(fp,sys,var);
622 FPRINTF(fp," \nhas nominal value of zero.\n");
623 FPRINTF(fp," Resetting to %g.\n",n);
624 var_set_nominal(var,n);
625 } else {
626 n = -n;
627 FPRINTF(fp,"ERROR: (slv8) calc_nominals\n");
628 FPRINTF(fp," Variable ");
629 print_var_name(fp,sys,var);
630 FPRINTF(fp," \nhas negative nominal value.\n");
631 FPRINTF(fp," Resetting to %g.\n",n);
632 var_set_nominal(var,n);
633 }
634 }
635 #if DEBUG
636 FPRINTF(fp,"Column %d is");
637 print_var_name(fp,sys,var);
638 FPRINTF(fp,"\nScaling of column %d is %g\n",col,n);
639 #endif /* DEBUG */
640 sys->nominals.vec[col] = n;
641 }
642 }
643 square_norm( &(sys->nominals) );
644 sys->update.nominals = UPDATE_NOMINALS;
645 sys->nominals.accurate = TRUE;
646 }
647
648
649 /*
650 * Calculates the weights of all of the block relations
651 * to scale the rows of the Jacobian.
652 */
653 static void calc_weights( slv8_system_t sys)
654 {
655 mtx_coord_t nz;
656 real64 sum;
657
658 if( sys->weights.accurate )
659 return;
660
661 nz.row = sys->weights.rng->low;
662 if(strcmp(SCALEOPT,"NONE") == 0 ||
663 strcmp(SCALEOPT,"ITERATIVE") == 0) {
664 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
665 sys->weights.vec[nz.row] = 1;
666 }
667 } else if (strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
668 strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0) {
669 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
670 sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
671 sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
672 }
673 } else if (strcmp(SCALEOPT,"RELNOM") == 0 ||
674 strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0) {
675 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
676 sys->weights.vec[nz.row] =
677 1.0/rel_nominal(sys->rlist[mtx_row_to_org(sys->J.mtx,nz.row)]);
678 }
679 }
680 square_norm( &(sys->weights) );
681 sys->update.weights = UPDATE_WEIGHTS;
682 sys->residuals.accurate = FALSE;
683 sys->weights.accurate = TRUE;
684 }
685
686
687 /*
688 * Scales the jacobian.
689 */
690 static void scale_J( slv8_system_t sys)
691 {
692 int32 row;
693 int32 col;
694
695 if( sys->J.accurate ) return;
696
697 calc_nominals(sys);
698 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ )
699 mtx_mult_col(sys->J.mtx,col,sys->nominals.vec[col],&(sys->J.reg.row));
700
701 calc_weights(sys);
702 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ )
703 mtx_mult_row(sys->J.mtx,row,sys->weights.vec[row],&(sys->J.reg.col));
704 }
705
706 static void jacobian_scaled(slv8_system_t sys)
707 {
708 int32 col;
709 if (DUMPCNORM) {
710 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
711 real64 cnorm;
712 cnorm =
713 calc_sqrt_D0(mtx_sum_sqrs_in_col(sys->J.mtx,col,&(sys->J.reg.row)));
714 if (cnorm >CNHIGH || cnorm <CNLOW) {
715 FPRINTF(MIF(sys),"[col %d org %d] %g\n", col,
716 mtx_col_to_org(sys->J.mtx,col), cnorm);
717 }
718 }
719 }
720
721 sys->update.jacobian = UPDATE_JACOBIAN;
722 sys->J.accurate = TRUE;
723 sys->J.singular = FALSE; /* yet to be determined */
724 #if DEBUG
725 FPRINTF(LIF(sys),"\nJacobian: \n");
726 debug_out_jacobian(LIF(sys),sys);
727 #endif /* DEBUG */
728 }
729
730
731
732 static void scale_variables( slv8_system_t sys)
733 {
734 int32 col;
735
736 if( sys->variables.accurate ) return;
737
738 col = sys->variables.rng->low;
739 for( ; col <= sys->variables.rng->high; col++ ) {
740 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
741 sys->variables.vec[col] = var_value(var)/sys->nominals.vec[col];
742 }
743 square_norm( &(sys->variables) );
744 sys->variables.accurate = TRUE;
745 #if DEBUG
746 FPRINTF(LIF(sys),"Variables: ");
747 debug_out_vector(LIF(sys),sys,&(sys->variables));
748 #endif /* DEBUG */
749 }
750
751
752 /*
753 * Scales the previously calculated residuals.
754 */
755 static void scale_residuals( slv8_system_t sys)
756 {
757 int32 row;
758
759 if( sys->residuals.accurate ) return;
760
761 row = sys->residuals.rng->low;
762 for( ; row <= sys->residuals.rng->high; row++ ) {
763 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
764 sys->residuals.vec[row] = rel_residual(rel)*sys->weights.vec[row];
765 }
766 square_norm( &(sys->residuals) );
767 sys->residuals.accurate = TRUE;
768 #if DEBUG
769 FPRINTF(LIF(sys),"Residuals: ");
770 debug_out_vector(LIF(sys),sys,&(sys->residuals));
771 #endif /* DEBUG */
772 }
773
774
775 /*
776 * Calculates relnoms for all relations in sys
777 * using variable nominals.
778 */
779 static void calc_relnoms(slv8_system_t sys)
780 {
781 int32 row, col;
782 struct var_variable *var;
783 struct rel_relation *rel;
784 real64 *var_list;
785 var_list = create_array(sys->cap,real64);
786 col = 0;
787 var = sys->vlist[col];
788 /* store current variable values and
789 set variable value to nominal value */
790 while(var != NULL){
791 var_list[col] = var_value(var);
792 var_set_value(var, var_nominal(var));
793 col++;
794 var = sys->vlist[col];
795 }
796 row = 0;
797 rel = sys->rlist[row];
798 /* calculate relation nominal */
799 while(rel != NULL){
800 relman_scale(rel);
801 row++;
802 rel = sys->rlist[row];
803 }
804 col = 0;
805 var = sys->vlist[col];
806 /* restore variable values */
807 while(var != NULL){
808 var_set_value(var, var_list[col]);
809 col++;
810 var = sys->vlist[col];
811 }
812 destroy_array(var_list);
813 }
814
815 /*
816 * Returns the maximum ratio of magnitudes of any two nonzero
817 * elements in the same column of mtx. Only considers elements
818 * in region reg.
819 */
820 static real64 col_max_ratio(mtx_matrix_t *mtx,
821 mtx_region_t *reg)
822 {
823 real64 ratio;
824 real64 max_ratio;
825 real64 num, denom, dummy;
826 mtx_coord_t coord;
827 max_ratio = 0;
828 for(coord.col = reg->col.low;
829 coord.col <= reg->col.high; coord.col++) {
830 ratio = 0;
831 num = mtx_col_max(*mtx,&(coord),&(reg->row),&(dummy));
832 denom = mtx_col_min(*mtx,&(coord),&(reg->row),&(dummy),1e-7);
833 if(denom >0){
834 ratio = num/denom;
835 }
836 if(ratio > 10000000){
837 /* FPRINTF(stderr,"HELPME\n");*/
838 }
839 if(ratio > max_ratio){
840 max_ratio = ratio;
841 }
842 }
843 if(max_ratio == 0){
844 max_ratio = 1;
845 }
846 return max_ratio;
847 }
848
849
850 /*
851 * Returns the maximum ratio of magnitudes of any two nonzero
852 * elements in the same row of mtx. Only considers elements
853 * in region reg.
854 */
855 static real64 row_max_ratio(mtx_matrix_t *mtx,
856 mtx_region_t *reg)
857 {
858 real64 ratio;
859 real64 max_ratio;
860 real64 num, denom, dummy;
861 mtx_coord_t coord;
862 max_ratio = 0;
863 for(coord.row = reg->row.low;
864 coord.row <= reg->row.high; coord.row++) {
865 ratio = 0;
866 num = mtx_row_max(*mtx,&(coord),&(reg->col),&(dummy));
867 denom = mtx_row_min(*mtx,&(coord),&(reg->col),&(dummy),1e-7);
868 if(denom >0){
869 ratio = num/denom;
870 }
871 if(ratio > 10000000){
872 /* FPRINTF(stderr,"HELPME\n");*/
873 }
874 if(ratio > max_ratio){
875 max_ratio = ratio;
876 }
877 }
878 if(max_ratio == 0){
879 max_ratio = 1;
880 }
881 return max_ratio;
882 }
883
884 /*
885 * Calculates scaling factor suggested by Fourer.
886 * For option = 0, returns scaling factor for
887 * row number loc.
888 * For option = 1, returns scaling factor for
889 * col number loc.
890 */
891 static real64 calc_fourer_scale(mtx_matrix_t mtx,
892 mtx_region_t reg,
893 int32 loc,
894 int32 option)
895 {
896 mtx_coord_t coord;
897 real64 min, max, dummy;
898 real64 scale;
899 if(option == 0){
900 if((loc < reg.row.low) || (loc > reg.row.high)){
901 return 1;
902 }
903 coord.row = loc;
904 min = mtx_row_min(mtx,&(coord),&(reg.col),&(dummy),1e-7);
905 max = mtx_row_max(mtx,&(coord),&(reg.col),&(dummy));
906 scale = min*max;
907 if(scale > 0){
908 scale = sqrt(scale);
909 } else {
910 scale = 1;
911 }
912 return scale;
913 } else {
914 if(loc < reg.col.low || loc > reg.col.high){
915 return 1;
916 }
917 coord.col = loc;
918 min = mtx_col_min(mtx,&(coord),&(reg.row),&(dummy),1e-7);
919 max = mtx_col_max(mtx,&(coord),&(reg.row),&(dummy));
920 scale = min*max;
921 if(scale > 0){
922 scale = sqrt(scale);
923 } else {
924 scale = 1;
925 }
926 return scale;
927 }
928 }
929
930 /*
931 * This funcion is an implementation of the scaling
932 * routine by Fourer on p304 of Mathematical Programing
933 * vol 23, (1982).
934 * This function will scale the Jacobian and store the scaling
935 * factors in sys->nominals and sys->weights.
936 * If the Jacobian has been previously scaled
937 * by another method (during this iteration) then these vectors
938 * should contain the scale factors used in that scaling.
939 */
940 static void scale_J_iterative(slv8_system_t sys)
941 {
942 real64 rho_col_old, rho_col_new;
943 real64 rho_row_old, rho_row_new;
944 int32 k;
945 int32 done;
946 int32 row, col;
947 real64 *colvec = sys->nominals.vec;
948 real64 *rowvec = sys->weights.vec;
949 real64 rowscale, colscale;
950 rho_col_old = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
951 rho_row_old = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
952 k = 0;
953 done = 0;
954 while(done == 0){
955 k++;
956 for(row = sys->J.reg.row.low;
957 row <= sys->J.reg.row.high; row++){
958 rowscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,row,0);
959 mtx_mult_row(sys->J.mtx,row,rowscale,&(sys->J.reg.col));
960 rowvec[row] *= rowscale;
961 }
962 for(col = sys->J.reg.col.low;
963 col <= sys->J.reg.col.high; col++){
964 colscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,col,1);
965 mtx_mult_col(sys->J.mtx,col,colscale,&(sys->J.reg.row));
966 colvec[col] *= colscale;
967 }
968 rho_col_new = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
969 rho_row_new = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
970 if((rho_col_new >= ITSCALETOL*rho_col_old &&
971 rho_row_new >= ITSCALETOL*rho_row_old)
972 || k >= ITSCALELIM){
973 done = 1;
974 /* FPRINTF(MIF(sys),"%d ITERATIVE SCALING ITERATIONS\n",k);*/
975 } else {
976 rho_row_old = rho_row_new;
977 rho_col_old = rho_col_new;
978 }
979 }
980 square_norm( &(sys->nominals) );
981 sys->update.nominals = UPDATE_NOMINALS;
982 sys->nominals.accurate = TRUE;
983
984 square_norm( &(sys->weights) );
985 sys->update.weights = UPDATE_WEIGHTS;
986 sys->residuals.accurate = FALSE;
987 sys->weights.accurate = TRUE;
988 }
989
990
991 /*
992 * Scale system dependent on interface parameters
993 */
994 static void scale_system( slv8_system_t sys )
995 {
996 if(strcmp(SCALEOPT,"NONE") == 0){
997 if(sys->J.accurate == FALSE){
998 calc_nominals(sys);
999 calc_weights(sys);
1000 jacobian_scaled(sys);
1001 }
1002 scale_variables(sys);
1003 scale_residuals(sys);
1004 return;
1005 }
1006 if(strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
1007 strcmp(SCALEOPT,"RELNOM") == 0){
1008 if(sys->J.accurate == FALSE){
1009 scale_J(sys);
1010 jacobian_scaled(sys);
1011 }
1012 scale_variables(sys);
1013 scale_residuals(sys);
1014 return;
1015 }
1016 if(strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0 ||
1017 strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0){
1018 if(sys->J.accurate == FALSE){
1019 --sys->update.iterative;
1020 if(sys->update.iterative <= 0) {
1021 scale_J(sys);
1022 scale_J_iterative(sys);
1023 sys->update.iterative =
1024 UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1025 } else {
1026 sys->weights.accurate = TRUE;
1027 sys->nominals.accurate = TRUE;
1028 scale_J(sys); /* will use current scaling vectors */
1029 }
1030 jacobian_scaled(sys);
1031 }
1032 scale_variables(sys);
1033 scale_residuals(sys);
1034 return;
1035 }
1036 if(strcmp(SCALEOPT,"ITERATIVE") == 0){
1037 if(sys->J.accurate == FALSE){
1038 --sys->update.iterative;
1039 if(sys->update.iterative <= 0) {
1040 calc_nominals(sys);
1041 calc_weights(sys);
1042 scale_J_iterative(sys);
1043 sys->update.iterative =
1044 UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1045 } else {
1046 sys->weights.accurate = TRUE;
1047 sys->nominals.accurate = TRUE;
1048 scale_J(sys); /* will use current scaling vectors */
1049 }
1050 jacobian_scaled(sys);
1051 }
1052 scale_variables(sys);
1053 scale_residuals(sys);
1054 }
1055 return;
1056 }
1057
1058
1059 /*
1060 * Resets all flags to setup a new solve.
1061 * Should set sys->s.block.current_block = -1
1062 * before calling.
1063 * This is currently a HACK!
1064 * not sure if should call when done.
1065 */
1066 static void conopt_initialize( slv8_system_t sys)
1067 {
1068
1069 sys->s.block.current_block++;
1070 /*
1071 * Next line was added to create the aray cost, whis is used by
1072 * the interface to display residuals and number of iterations
1073 */
1074 sys->s.costsize = 1+sys->s.block.number_of;
1075
1076 if( sys->s.block.current_block < sys->s.block.number_of ) {
1077 boolean ok;
1078
1079 sys->s.block.iteration = 0;
1080 sys->s.block.cpu_elapsed = 0.0;
1081 sys->s.block.functime = 0.0;
1082 sys->s.block.jactime = 0.0;
1083 sys->s.block.funcs = 0;
1084 sys->s.block.jacs = 0;
1085
1086 sys->s.calc_ok = TRUE;
1087
1088 if(sys->p.output.less_important && (LIFDS ||
1089 sys->s.block.current_size > 1)) {
1090 debug_delimiter(LIF(sys));
1091 debug_delimiter(LIF(sys));
1092 }
1093 if(sys->p.output.less_important && LIFDS) {
1094 FPRINTF(LIF(sys),"\n%-40s ---> %d in [%d..%d]\n",
1095 "Current block number", sys->s.block.current_block,
1096 0, sys->s.block.number_of-1);
1097 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Current block size",
1098 sys->s.block.current_size);
1099 }
1100 if( !(ok = calc_objective(sys)) ) {
1101 FPRINTF(MIF(sys),"Objective calculation errors detected.\n");
1102 }
1103 if(sys->p.output.less_important && sys->obj) {
1104 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Objective", sys->objective);
1105 }
1106 sys->s.calc_ok = sys->s.calc_ok && ok;
1107
1108 if (!(sys->p.ignore_bounds) ) {
1109 slv_insure_bounds(SERVER, sys->J.reg.col.low,
1110 sys->J.reg.col.high,MIF(sys));
1111 }
1112
1113 sys->residuals.accurate = FALSE;
1114 if( !(ok = calc_residuals(sys)) ) {
1115 FPRINTF(MIF(sys),"Residual calculation errors detected.\n");
1116 }
1117 if( sys->p.output.less_important &&
1118 (sys->s.block.current_size >1 ||
1119 LIFDS) ) {
1120 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
1121 sys->s.block.residual);
1122 }
1123 sys->s.calc_ok = sys->s.calc_ok && ok;
1124
1125 /* Must be updated as soon as required */
1126 sys->J.accurate = FALSE;
1127 sys->update.weights = 0;
1128 sys->update.nominals = 0;
1129 sys->update.relnoms = 0;
1130 sys->update.iterative = 0;
1131 sys->variables.accurate = FALSE;
1132 }
1133 }
1134
1135
1136 /*
1137 * Iteration begin/end routines
1138 * ----------------------------
1139 * iteration_begins(sys)
1140 * iteration_ends(sys)
1141 */
1142
1143 /*
1144 * Prepares sys for entering an iteration, increasing the iteration counts
1145 * and starting the clock.
1146 */
1147 static void iteration_begins( slv8_system_t sys)
1148 {
1149 sys->clock = tm_cpu_time();
1150 ++(sys->s.block.iteration);
1151 ++(sys->s.iteration);
1152 if(sys->p.output.less_important && LIFDS) {
1153 FPRINTF(LIF(sys),"\n%-40s ---> %d\n",
1154 "Iteration", sys->s.block.iteration);
1155 FPRINTF(LIF(sys),"%-40s ---> %d\n",
1156 "Total iteration", sys->s.iteration);
1157 }
1158 }
1159
1160
1161 /*
1162 * Prepares sys for exiting an iteration, stopping the clock and recording
1163 * the cpu time.
1164 */
1165 static void iteration_ends( slv8_system_t sys)
1166 {
1167 double cpu_elapsed; /* elapsed this iteration */
1168
1169 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
1170 sys->s.block.cpu_elapsed += cpu_elapsed;
1171 sys->s.cpu_elapsed += cpu_elapsed;
1172 if(sys->p.output.less_important && LIFDS) {
1173 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1174 "Elapsed time", sys->s.block.cpu_elapsed);
1175 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1176 "Total elapsed time", sys->s.cpu_elapsed);
1177 }
1178 }
1179
1180
1181 /*
1182 * Updates the solver status.
1183 */
1184 static void update_status( slv8_system_t sys)
1185 {
1186 boolean unsuccessful;
1187
1188 if( !sys->s.converged ) {
1189 sys->s.time_limit_exceeded =
1190 (sys->s.block.cpu_elapsed >= TIME_LIMIT);
1191 sys->s.iteration_limit_exceeded =
1192 (sys->s.block.iteration >= ITER_LIMIT);
1193 }
1194
1195 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
1196 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
1197
1198 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
1199 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
1200 }
1201
1202
1203 static
1204 int32 slv8_get_default_parameters(slv_system_t server, SlvClientToken asys,
1205 slv_parameters_t *parameters)
1206 {
1207 slv8_system_t sys;
1208 union parm_arg lo,hi,val;
1209 struct slv_parameter *new_parms = NULL;
1210 static char *reorder_names[] = {
1211 "SPK1","TEAR_DROP","OVER_TEAR"
1212 };
1213 static char *scaling_names[] = {
1214 "NONE","ROW_2NORM","RELNOM" /*,"2NORM+ITERATIVE",
1215 "RELNOM+ITERATIVE","ITERATIVE" */
1216 };
1217 int32 make_macros = 0;
1218 if (server != NULL && asys != NULL) {
1219 sys = SLV8(asys);
1220 make_macros = 1;
1221 }
1222
1223 #if DEBUG /* keep purify from whining on UMR */
1224 lo.argr = hi.argr = val.argr = 0.0;
1225 #endif /* DEBUG */
1226
1227 if (parameters->parms == NULL) {
1228 /* an external client wants our parameter list.
1229 * an instance of slv8_system_structure has this pointer
1230 * already set in slv8_create
1231 */
1232 new_parms = (struct slv_parameter *)
1233 ascmalloc((slv8_PA_SIZE)*sizeof(struct slv_parameter));
1234 if (new_parms == NULL) {
1235 return -1;
1236 }
1237 parameters->parms = new_parms;
1238 parameters->dynamic_parms = 1;
1239 }
1240 parameters->num_parms = 0;
1241
1242 /* begin defining parameters */
1243
1244 slv_define_parm(parameters, real_parm,
1245 "infinity","RTMAXV","Internal value of infinity",
1246 U_p_real(val,10e20),U_p_real(lo,10),U_p_real(hi,MAX_REAL),2);
1247
1248 slv_define_parm(parameters, real_parm,
1249 "maxjac","RTMAXJ","Maximum Jacobian Element",
1250 U_p_real(val,2e8),U_p_real(lo,10),U_p_real(hi,MAX_REAL),2);
1251 SLV_RPARM_MACRO(RTMAXJ_PTR,parameters);
1252
1253 slv_define_parm(parameters, real_parm,
1254 "hessian_ub","RTMXJ2","upper bound on 2nd derivatives",
1255 U_p_real(val,1e4),U_p_real(lo,0),U_p_real(hi,MAX_REAL),2);
1256
1257 slv_define_parm(parameters, real_parm,
1258 "maxfeastol", "RTNWMA",
1259 "max residual considered feasible (may be scaled)",
1260 U_p_real(val, 1e-3),U_p_real(lo, 1e-13),U_p_real(hi,10e10),2);
1261
1262 slv_define_parm(parameters, real_parm,
1263 "minfeastol", "RTNWMI",
1264 "residuals below this always considered feasible",
1265 U_p_real(val, 4e-10),U_p_real(lo, 1e-20),U_p_real(hi,10e10),2);
1266
1267 slv_define_parm(parameters, real_parm,
1268 "oneDsearch","RTONED","accuracy of one dimensional search",
1269 U_p_real(val,0.2),U_p_real(lo,0.1),U_p_real(hi,0.7),2);
1270
1271 slv_define_parm(parameters, real_parm,
1272 "stepmult","RVSTLM","steplength multiplier",
1273 U_p_real(val,4),U_p_real(lo,0),U_p_real(hi,MAX_REAL),2);
1274
1275 slv_define_parm(parameters, real_parm,
1276 "objtol","RTOBJR","relative objective tolerance",
1277 U_p_real(val,3e-13),U_p_real(lo,0),U_p_real(hi,1),2);
1278
1279 slv_define_parm(parameters, real_parm,
1280 "pivottol","RTPIVA","absolute pivot tolerance",
1281 U_p_real(val,1e-7),U_p_real(lo,1e-15),U_p_real(hi,1),2);
1282
1283 slv_define_parm(parameters, real_parm,
1284 "pivtolrel","RTPIVR","relative pivot tolerance",
1285 U_p_real(val,0.05),U_p_real(lo,0),U_p_real(hi,1),2);
1286
1287 slv_define_parm(parameters, real_parm,
1288 "opttol","RTREDG","optimality tolerance",
1289 U_p_real(val,2e-5),U_p_real(lo,0),U_p_real(hi,MAX_REAL),2);
1290
1291 slv_define_parm(parameters, int_parm,
1292 "log_freq","LFILOG","Log Frequency",
1293 U_p_int(val,10),U_p_int(lo,1),U_p_int(hi,MAX_INT),1);
1294
1295 slv_define_parm(parameters, int_parm,
1296 "iterationlimit", "LFITER", "maximum number of iterations",
1297 U_p_int(val, 1000),U_p_int(lo, 1),U_p_int(hi,MAX_INT),1);
1298 SLV_IPARM_MACRO(ITER_LIMIT_PTR,parameters);
1299
1300 slv_define_parm(parameters, int_parm,
1301 "slowproglim","LFNICR","limit for slow progress",
1302 U_p_int(val,5),U_p_int(lo,1),U_p_int(hi,MAX_INT),1);
1303
1304 slv_define_parm(parameters, int_parm,
1305 "maxhessdim","LFNICR","maximum Hessian dimension",
1306 U_p_int(val,500),U_p_int(lo,1),U_p_int(hi,MAX_INT),1);
1307
1308 slv_define_parm(parameters, int_parm,
1309 "supbasiclim","LFMXNS","limit on new superbasics",
1310 U_p_int(val,5),U_p_int(lo,0),U_p_int(hi,MAX_INT),1);
1311
1312 slv_define_parm(parameters, int_parm,
1313 "errlim","max # func errs",
1314 "limit on function evaluation errors",
1315 U_p_int(val,500),U_p_int(lo,0),U_p_int(hi,MAX_INT),1);
1316 SLV_IPARM_MACRO(DOMLIM_PTR,parameters);
1317
1318 slv_define_parm(parameters, char_parm,
1319 "scaleopt", "scaling option", "scaling option",
1320 U_p_string(val,scaling_names[1]),
1321 U_p_strings(lo,scaling_names),
1322 U_p_int(hi,sizeof(scaling_names)/sizeof(char *)),3);
1323 SLV_CPARM_MACRO(SCALEOPT_PTR,parameters);
1324
1325 slv_define_parm(parameters, bool_parm,
1326 "lifds", "show singletons details", "show singletons details",
1327 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 3);
1328 SLV_BPARM_MACRO(LIFDS_PTR,parameters);
1329
1330 slv_define_parm(parameters, bool_parm,
1331 "safe_calc", "safe calculations", "safe calculations",
1332 U_p_bool(val, 1),U_p_bool(lo,0),U_p_bool(hi,1), 3);
1333 SLV_BPARM_MACRO(SAFE_CALC_PTR,parameters);
1334
1335 slv_define_parm(parameters, real_parm,
1336 "toosmall", "default for zero nominal",
1337 "default for zero nominal",
1338 U_p_real(val, 1e-8),U_p_real(lo, 1e-12),U_p_real(hi,1.0), 3);
1339 SLV_RPARM_MACRO(TOO_SMALL_PTR,parameters);
1340
1341 slv_define_parm(parameters, int_parm,
1342 "upwts", "Row scaling update frequency",
1343 "Row scaling update frequency",
1344 U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
1345 SLV_IPARM_MACRO(UPDATE_WEIGHTS_PTR,parameters);
1346
1347 slv_define_parm(parameters, int_parm,
1348 "upnom", "Column scaling update frequency",
1349 "Column scaling update frequency",
1350 U_p_int(val, 1000),U_p_int(lo,0),U_p_int(hi,20000), 3);
1351 SLV_IPARM_MACRO(UPDATE_NOMINALS_PTR,parameters);
1352
1353 slv_define_parm(parameters, bool_parm,
1354 "cncols", "Check poorly scaled columns",
1355 "Check poorly scaled columns",
1356 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 3);
1357 SLV_BPARM_MACRO(DUMPCNORM_PTR,parameters);
1358
1359 slv_define_parm(parameters, real_parm,
1360 "cnlow", "smallest allowable column norm",
1361 "smallest allowable column norm",
1362 U_p_real(val, 0.01),U_p_real(lo, 0),U_p_real(hi,10e10), 3);
1363 SLV_RPARM_MACRO(CNLOW_PTR,parameters);
1364
1365 slv_define_parm(parameters, real_parm,
1366 "cnhigh", "largest allowable column norm",
1367 "largest allowable column norm",
1368 U_p_real(val, 100.0),U_p_real(lo,0),U_p_real(hi,10e10), 3);
1369 SLV_RPARM_MACRO(CNHIGH_PTR,parameters);
1370
1371 slv_define_parm(parameters, int_parm,
1372 "upjac", "Jacobian update frequency",
1373 "Jacobian update frequency",
1374 U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
1375 SLV_IPARM_MACRO(UPDATE_JACOBIAN_PTR,parameters);
1376
1377 slv_define_parm(parameters, int_parm,
1378 "itscalelim", "Iteration lim for iterative scale",
1379 "Iteration lim for iterative scale",
1380 U_p_int(val, 10),U_p_int(lo,0),U_p_int(hi,20000), 3);
1381 SLV_IPARM_MACRO(ITSCALELIM_PTR,parameters);
1382
1383 slv_define_parm(parameters, real_parm,
1384 "itscaletol", "Iterative scaling tolerance",
1385 "Iterative scaling tolerance",
1386 U_p_real(val, 0.99999),U_p_real(lo,0),U_p_real(hi,1.0), 3);
1387 SLV_RPARM_MACRO(ITSCALETOL_PTR,parameters);
1388
1389 slv_define_parm(parameters, char_parm,
1390 "reorder", "reorder method", "reorder method",
1391 U_p_string(val,reorder_names[0]),
1392 U_p_strings(lo,reorder_names),
1393 U_p_int(hi,sizeof(reorder_names)/sizeof(char *)),3);
1394 SLV_CPARM_MACRO(REORDER_OPTION_PTR,parameters);
1395
1396 slv_define_parm(parameters, int_parm,
1397 "cutoff", "block size cutoff (MODEL-based)",
1398 "block size cutoff (MODEL-based)",
1399 U_p_int(val, 200),U_p_int(lo,0),U_p_int(hi,20000), 3);
1400 SLV_IPARM_MACRO(CUTOFF_PTR,parameters);
1401
1402 slv_define_parm(parameters, bool_parm,
1403 "relnomscale", "calc rel nominals", "calc rel nominals",
1404 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 3);
1405 SLV_BPARM_MACRO(RELNOMSCALE_PTR,parameters);
1406
1407 slv_define_parm(parameters, int_parm,
1408 "timelimit", "time limit (CPU sec/block)",
1409 "time limit (CPU sec/block)",
1410 U_p_int(val,1500),U_p_int(lo, 1),U_p_int(hi,20000),1);
1411 SLV_IPARM_MACRO(TIME_LIMIT_PTR,parameters);
1412
1413 return 1;
1414 }
1415
1416
1417 /*
1418 * External routines
1419 * -----------------
1420 * See slv_client.h
1421 */
1422
1423 static SlvClientToken slv8_create(slv_system_t server, int32*statusindex)
1424 {
1425 slv8_system_t sys;
1426
1427 sys = (slv8_system_t)asccalloc(1, sizeof(struct slv8_system_structure) );
1428 if (sys==NULL) {
1429 *statusindex = 1;
1430 return sys;
1431 }
1432 SERVER = server;
1433 sys->p.parms = sys->pa;
1434 sys->p.dynamic_parms = 0;
1435 slv8_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
1436 sys->integrity = OK;
1437 sys->presolved = 0;
1438 sys->resolve = 0;
1439 sys->p.output.more_important = stdout;
1440 sys->p.output.less_important = stdout;
1441
1442 sys->p.whose = (*statusindex);
1443
1444 sys->s.ok = TRUE;
1445 sys->s.calc_ok = TRUE;
1446 sys->s.costsize = 0;
1447 sys->s.cost = NULL; /*redundant, but sanity preserving */
1448 sys->vlist = slv_get_solvers_var_list(server);
1449 sys->rlist = slv_get_solvers_rel_list(server);
1450 sys->obj = slv_get_obj_relation(server);
1451 if (sys->vlist == NULL) {
1452 ascfree(sys);
1453 FPRINTF(MIF(sys),"CONOPT called with no variables.\n");
1454 *statusindex = -2;
1455 return NULL; /* prolly leak here */
1456 }
1457 if (sys->rlist == NULL && sys->obj == NULL) {
1458 ascfree(sys);
1459 FPRINTF(MIF(sys),"CONOPT called with no relations or objective.\n");
1460 *statusindex = -1;
1461 return NULL; /* prolly leak here */
1462 }
1463 /* we don't give a damn about the objective list or the pars or
1464 * bounds or extrels or any of the other crap.
1465 */
1466 slv_check_var_initialization(server);
1467 *statusindex = 0;
1468 return((SlvClientToken)sys);
1469 }
1470
1471 static void destroy_matrices( slv8_system_t sys)
1472 {
1473 if( sys->J.mtx ) {
1474 mtx_destroy(sys->J.mtx);
1475 }
1476 }
1477
1478 static void destroy_vectors( slv8_system_t sys)
1479 {
1480 destroy_array(sys->nominals.vec);
1481 destroy_array(sys->weights.vec);
1482 destroy_array(sys->relnoms.vec);
1483 destroy_array(sys->variables.vec);
1484 destroy_array(sys->residuals.vec);
1485 }
1486
1487
1488 static int32 slv8_eligible_solver(slv_system_t server)
1489 {
1490 struct rel_relation **rp;
1491 for( rp=slv_get_solvers_rel_list(server); *rp != NULL ; ++rp ) {
1492 if( rel_less(*rp) || rel_greater(*rp) ) return(FALSE);
1493 }
1494 return(TRUE);
1495 }
1496
1497
1498 static void slv8_get_parameters(slv_system_t server, SlvClientToken asys,
1499 slv_parameters_t *parameters)
1500 {
1501 slv8_system_t sys;
1502 (void)server; /* stop gcc whine about unused parameter */
1503
1504 sys = SLV8(asys);
1505 if (check_system(sys)) return;
1506 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
1507 }
1508
1509
1510 static void slv8_set_parameters(slv_system_t server, SlvClientToken asys,
1511 slv_parameters_t *parameters)
1512 {
1513 slv8_system_t sys;
1514 (void)server; /* stop gcc whine about unused parameter */
1515
1516 sys = SLV8(asys);
1517 if (check_system(sys)) return;
1518 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
1519 }
1520
1521
1522 static void slv8_get_status(slv_system_t server, SlvClientToken asys,
1523 slv_status_t *status)
1524 {
1525 slv8_system_t sys;
1526 (void)server; /* stop gcc whine about unused parameter */
1527
1528 sys = SLV8(asys);
1529 if (check_system(sys)) return;
1530 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
1531 }
1532
1533
1534 static linsolqr_system_t slv8_get_linsolqr_sys(slv_system_t server,
1535 SlvClientToken asys)
1536 {
1537 slv8_system_t sys;
1538 (void)server; /* stop gcc whine about unused parameter */
1539
1540 sys = SLV8(asys);
1541 if (check_system(sys)) return NULL;
1542 return(sys->J.sys);
1543 }
1544
1545
1546 static linsol_system_t slv8_get_linsol_sys(slv_system_t server,
1547 SlvClientToken asys)
1548 {
1549 (void)server; /* stop gcc whine about unused parameter */
1550 (void)asys; /* stop gcc whine about unused parameter */
1551 return( NULL );
1552 }
1553
1554
1555 /*
1556 * Performs structural analysis on the system, setting the flags in
1557 * status. The problem must be set up, the relation/variable list
1558 * must be non-NULL. The
1559 * jacobian (linear) system must be created and have the correct order
1560 * (stored in sys->cap). Everything else will be determined here.
1561 * On entry there isn't yet a correspondence between var_sindex and
1562 * jacobian column. Here we establish that.
1563 */
1564 static void structural_analysis(slv_system_t server, slv8_system_t sys)
1565 {
1566 /* this function has been striped of its guts for conopt
1567 and may go away */
1568
1569 var_filter_t vfilter;
1570 rel_filter_t rfilter;
1571
1572 /*
1573 * The server has marked incidence flags already.
1574 */
1575 /* count included equalities */
1576 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1577 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1578 sys->rused = slv_count_solvers_rels(server,&rfilter);
1579
1580 /* count free and incident vars */
1581 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
1582 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
1583 sys->vused = slv_count_solvers_vars(server,&vfilter);
1584
1585 /* Symbolic analysis */
1586 sys->rtot = slv_get_num_solvers_rels(server);
1587 sys->vtot = slv_get_num_solvers_vars(server);
1588
1589 /*
1590 * The next few lines are used to calculate the rank of the nonlinear
1591 * system. We need it to evaluate if the system is structurally
1592 * singular or not. Calculating this number will keep CONOPT from
1593 * displaying a "structurally singular" error message
1594 */
1595 if (sys->rtot) {
1596 slv_block_partition(server);
1597 }
1598 sys->J.dofdata = slv_get_dofdata(server);
1599 sys->rank = sys->J.dofdata->structural_rank;
1600 /*
1601 * Unify the partitions since we feed CONOPT with a single block.
1602 */
1603 slv_block_unify(server);
1604
1605
1606 sys->J.reg.row.low = sys->J.reg.col.low = 0;
1607 sys->J.reg.row.high = sys->con.m -1;
1608 if (sys->obj != NULL) sys->J.reg.row.high--;
1609 sys->J.reg.col.high = sys->con.n -1;
1610 slv_check_bounds(SERVER,sys->vused,sys->vtot-1,MIF(sys),"fixed");
1611
1612 /* Initialize Status */
1613 sys->s.over_defined = (sys->rused > sys->vused);
1614 sys->s.under_defined = (sys->rused < sys->vused);
1615 sys->s.struct_singular = (sys->rank < sys->rused);
1616 sys->s.block.number_of = (slv_get_solvers_blocks(SERVER))->nblocks;
1617
1618 }
1619
1620
1621 static void create_matrices(slv_system_t server, slv8_system_t sys)
1622 {
1623 sys->J.mtx = mtx_create();
1624 mtx_set_order(sys->J.mtx,sys->cap);
1625 structural_analysis(server,sys);
1626 }
1627
1628
1629 static void create_vectors(sys)
1630 slv8_system_t sys;
1631 {
1632 sys->nominals.vec = create_array(sys->cap,real64);
1633 sys->nominals.rng = &(sys->J.reg.col);
1634 sys->weights.vec = create_array(sys->cap,real64);
1635 sys->weights.rng = &(sys->J.reg.row);
1636 sys->relnoms.vec = create_array(sys->cap,real64);
1637 sys->relnoms.rng = &(sys->J.reg.row);
1638 sys->variables.vec = create_array(sys->cap,real64);
1639 sys->variables.rng = &(sys->J.reg.col);
1640 sys->residuals.vec = create_array(sys->cap,real64);
1641 sys->residuals.rng = &(sys->J.reg.row);
1642 }
1643
1644
1645 static void slv8_dump_internals(slv_system_t server,
1646 SlvClientToken sys,int32 level)
1647 {
1648 (void)server; /* stop gcc whine about unused parameter */
1649
1650 check_system(sys);
1651 if (level > 0) {
1652 FPRINTF(MI8F(sys),"ERROR: (slv8) slv8_dump_internals\n");
1653 FPRINTF(MI8F(sys)," slv8 does not dump its internals.\n");
1654 }
1655 }
1656
1657
1658 /*
1659 * Here we will check if any fixed or included flags have
1660 * changed since the last presolve.
1661 */
1662 static int32 slv8_dof_changed(slv8_system_t sys)
1663 {
1664 int32 ind, result = 0;
1665 /* Currently we have two copies of the fixed and included flags
1666 which must be kept in sync. The var_fixed and rel_included
1667 functions perform the syncronization and hence must be called
1668 over the whole var list and rel list respectively. When we move
1669 to using only one set of flags (bit flags) this function can
1670 be changed to return 1 at the first indication of a change
1671 in the dof. */
1672
1673 /* search for vars that were fixed and are now free */
1674 for( ind = sys->vused; ind < sys->vtot; ++ind ) {
1675 if( !var_fixed(sys->vlist[ind]) && var_active(sys->vlist[ind]) ) {
1676 ++result;
1677 }
1678 }
1679 /* search for rels that were unincluded and are now included */
1680 for( ind = sys->rused; ind < sys->rtot; ++ind ) {
1681 if( rel_included(sys->rlist[ind]) && rel_active(sys->rlist[ind])) {
1682 ++result;
1683 }
1684 }
1685 /* search for vars that were free and are now fixed */
1686 for( ind = sys->vused -1; ind >= 0; --ind ) {
1687 if( var_fixed(sys->vlist[ind]) || !var_active(sys->vlist[ind])) {
1688 ++result;
1689 }
1690 }
1691 /* search for rels that were included and are now unincluded */
1692 for( ind = sys->rused -1; ind >= 0; --ind ) {
1693 if( !rel_included(sys->rlist[ind]) || !rel_active(sys->rlist[ind]) ) {
1694 ++result;
1695 }
1696 }
1697 return result;
1698 }
1699
1700
1701 static void reset_cost(struct slv_block_cost *cost,int32 costsize)
1702 {
1703 int32 ind;
1704 for( ind = 0; ind < costsize; ++ind ) {
1705 cost[ind].size = 0;
1706 cost[ind].iterations = 0;
1707 cost[ind].funcs = 0;
1708 cost[ind].jacs = 0;
1709 cost[ind].functime = 0;
1710 cost[ind].jactime = 0;
1711 cost[ind].time = 0;
1712 cost[ind].resid = 0;
1713 }
1714 }
1715
1716
1717 /*
1718 * After running CONOPT, we need to update the values of the array
1719 * cost, which is used by the interface to display residual and number
1720 * of iterations
1721 */
1722 static void update_cost(slv8_system_t sys)
1723 {
1724 int32 ci;
1725 if (sys->s.cost == NULL) {
1726 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
1727 } else {
1728 reset_cost(sys->s.cost,sys->s.costsize);
1729 }
1730 ci=sys->s.block.current_block;
1731 sys->s.cost[ci].size = sys->s.block.current_size;
1732 sys->s.cost[ci].iterations = sys->s.block.iteration;
1733 sys->s.cost[ci].resid = sys->s.block.residual;
1734 }
1735
1736
1737 static void slv8_presolve(slv_system_t server, SlvClientToken asys)
1738 {
1739 struct var_variable **vp;
1740 struct rel_relation **rp;
1741 int32 cap, ind;
1742 int32 matrix_creation_needed = 1;
1743 slv8_system_t sys;
1744
1745 sys = SLV8(asys);
1746 iteration_begins(sys);
1747 check_system(sys);
1748 if( sys->vlist == NULL ) {
1749 FPRINTF(MIF(sys),"ERROR: (slv8) slv8_presolve\n");
1750 FPRINTF(MIF(sys)," Variable list was never set.\n");
1751 return;
1752 }
1753 if( sys->rlist == NULL && sys->obj == NULL ) {
1754 FPRINTF(MIF(sys),"ERROR: (slv8) slv8_presolve\n");
1755 FPRINTF(MIF(sys)," Relation list and objective never set.\n");
1756 return;
1757 }
1758
1759 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
1760
1761 if(sys->presolved > 0) { /* system has been presolved before */
1762 if(!slv8_dof_changed(sys) /*no changes in fixed or included flags*/
1763 && sys->p.partition == sys->J.old_partition
1764 && sys->obj == sys->old_obj) {
1765 #if DEBUG
1766 FPRINTF(MIF(sys),"YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION\n");
1767 #endif /* DEBUG */
1768 matrix_creation_needed = 0;
1769 }
1770 }
1771
1772 rp=sys->rlist;
1773 for( ind = 0; ind < sys->rtot; ++ind ) {
1774 rel_set_satisfied(rp[ind],FALSE);
1775 }
1776 if( matrix_creation_needed ) {
1777
1778 cap = slv_get_num_solvers_rels(SERVER);
1779 sys->cap = slv_get_num_solvers_vars(SERVER);
1780 sys->cap = MAX(sys->cap,cap);
1781 vp=sys->vlist;
1782 for( ind = 0; ind < sys->vtot; ++ind ) {
1783 var_set_in_block(vp[ind],FALSE);
1784 }
1785 rp=sys->rlist;
1786 for( ind = 0; ind < sys->rtot; ++ind ) {
1787 rel_set_in_block(rp[ind],FALSE);
1788 rel_set_satisfied(rp[ind],FALSE);
1789 }
1790
1791 sys->presolved = 1; /* full presolve recognized here */
1792 sys->resolve = 0; /* initialize resolve flag */
1793 sys->J.old_partition = sys->p.partition;
1794 sys->old_obj = sys->obj;
1795
1796 slv_sort_rels_and_vars(server,&(sys->con.m),&(sys->con.n));
1797 if (sys->obj != NULL) {
1798 sys->con.m++; /* treat objective as a row */
1799 }
1800
1801 sys->con.ipsz[0] = sys->con.n;
1802 sys->con.ipsz[1] = sys->con.m;
1803 sys->con.nz = num_jacobian_nonzeros(sys, &(sys->con.maxrow));
1804 sys->con.ipsz[2] = sys->con.nz;
1805 sys->con.nintgr = NINTGR;
1806 /*
1807 * Memory estimation by calling the CONOPT subroutine coimem
1808 * The use of conopt_estimate_memory is a hack to avoid
1809 * unresolved external during the linking of the CONOPT library.
1810 * See conopt.h
1811 */
1812 conopt_estimate_memory(&(sys->con.nintgr),&(sys->con.ipsz[0]),
1813 &(sys->con.minmem),&(sys->con.estmem));
1814
1815 if (sys->con.work != NULL) {
1816 sys->con.work = (real64 *)ascrealloc(sys->con.work,
1817 sys->con.estmem*sizeof(real64));
1818 } else {
1819 /* calloc here doesn't help the crash */
1820 sys->con.work = (real64 *)ascmalloc(sys->con.estmem*sizeof(real64));
1821 }
1822 sys->con.lwork = sys->con.estmem;
1823 destroy_vectors(sys);
1824 destroy_matrices(sys);
1825 create_matrices(server,sys);
1826 create_vectors(sys);
1827
1828 sys->s.block.current_reordered_block = -2;
1829 }
1830
1831 /* Reset status */
1832 sys->con.optimized = 0;
1833 sys->s.iteration = 0;
1834 sys->s.cpu_elapsed = 0.0;
1835 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
1836 sys->s.block.previous_total_size = 0;
1837 sys->s.costsize = 1+sys->s.block.number_of;
1838
1839 if( matrix_creation_needed ) {
1840 destroy_array(sys->s.cost);
1841 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
1842 for( ind = 0; ind < sys->s.costsize; ++ind ) {
1843 sys->s.cost[ind].reorder_method = -1;
1844 }
1845 } else {
1846 reset_cost(sys->s.cost,sys->s.costsize);
1847 }
1848
1849 /* set to go to first unconverged block */
1850 sys->s.block.current_block = -1;
1851 sys->s.block.current_size = 0;
1852 sys->s.calc_ok = TRUE;
1853 sys->s.block.iteration = 0;
1854 sys->objective = MAXDOUBLE/2000.0;
1855
1856 update_status(sys);
1857 iteration_ends(sys);
1858 sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;
1859 }
1860
1861
1862 static void slv8_resolve(slv_system_t server, SlvClientToken asys)
1863 {
1864 struct var_variable **vp;
1865 struct rel_relation **rp;
1866 slv8_system_t sys;
1867 (void)server; /* stop gcc whine about unused parameter */
1868
1869 sys = SLV8(asys);
1870
1871 check_system(sys);
1872 for( vp = sys->vlist ; *vp != NULL ; ++vp ) {
1873 var_set_in_block(*vp,FALSE);
1874 }
1875 for( rp = sys->rlist ; *rp != NULL ; ++rp ) {
1876 rel_set_in_block(*rp,FALSE);
1877 rel_set_satisfied(*rp,FALSE);
1878 }
1879
1880 sys->resolve = 1; /* resolved recognized here */
1881
1882 /* Reset status */
1883 sys->s.iteration = 0;
1884 sys->s.cpu_elapsed = 0.0;
1885 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
1886 sys->s.block.previous_total_size = 0;
1887
1888 /* go to first unconverged block */
1889 sys->s.block.current_block = -1;
1890 sys->s.block.current_size = 0;
1891 sys->s.calc_ok = TRUE;
1892 sys->s.block.iteration = 0;
1893 sys->objective = MAXDOUBLE/2000.0;
1894
1895 update_status(sys);
1896 }
1897
1898
1899 /*
1900 * CONOPT ROUTINES
1901 */
1902
1903 /*
1904 * WORK TO DO:
1905 * Fix interface so that solvers define status messages. We
1906 * should not be stuck with one standard set that all solvers
1907 * must deal with.
1908 * Reimplement old code to detect linear coefficients and use
1909 * in conopt hookup.
1910 * Implement better communication routines.
1911 * Give more contol to interface (ex. turn off error counting,
1912 * switch from we alocate to conopt allocates, etc.).
1913 * Record marginal values and make available to user.
1914 * Set up interface such that any variable can be selected and
1915 * either maximized or minimized.
1916 * Get rid of our Jacobian calculation routine and stuff conopt's
1917 * workspace directly (in unsorted order). This will require
1918 * rewriting the scaling functions. (This code really should
1919 * not be in the solver)
1920 * Fix up restart code...we don't keep track of which values
1921 * change so must update entire Jacobian and residual vector
1922 * but may save some time by not having to redo column pointers etc.
1923 * Implement a way to bailout...check for ^C and tell conopt to
1924 * return as soon as possible.
1925 * Currently will not take problem like MIN x^2...it will complain
1926 * about empty rows. Must formulate as y=x^2; MIN y; until we
1927 * fix the way we handle objectives.
1928 */
1929
1930
1931 /*
1932 * COIRMS Based on the information provided in Coispz, CONOPT will
1933 * allocate the number of vectors into which the user can define
1934 * the details of the model. The details of the model are defined
1935 * here.
1936 *
1937 * COIRMS(lower, curr, upper, vsta, type,rhs, fv, esta, colsta,
1938 * rowno, value, nlflag, n, m, n1, nz, usrmem)
1939 *
1940 * lower - lower bounds on the variables
1941 * curr - intial values of the variables
1942 * upper - upper bounds on the variables
1943 * vsta - initial status of the variable(o nonbasic, 1 basic)
1944 * type - types of equations (0 equality,1 greater,2 less)
1945 * rhs - values of the right hand sides
1946 * fv - sum of the nonlinear terms in the initial point
1947 * esta - initial status of the slack in the constraint (nonbasic,basic)
1948 * colsta- start of column pointers
1949 * rowno - row or equation numbers of the nonzeros
1950 * value - values of the jacobian elements
1951 * nlflag- nonlinearity flags(0 nonzero constant,1 varying)
1952 * n - number of variables
1953 * m - number of constraints
1954 * n1 - n+1
1955 * nz - number of jacobian elements
1956 * usrmem- user memory defined by conopt
1957 */
1958 static void slv8_coirms(real64 *lower, real64 *curr, real64 *upper,
1959 int32 *vsta, int32 *type, real64 *rhs,
1960 real64 *fv, int32 *esta, int32 *colsta,
1961 int32 *rowno, real64 *value, int32 *nlflag,
1962 int32 *n, int32 *m, int32 *n1, int32 *nz,
1963 real64 *usrmem)
1964 {
1965 int32 col,row,count,count_old,len,c,r,offset, obj_count;
1966 real64 nominal, up, low;
1967 struct var_variable *var;
1968 struct rel_relation *rel;
1969 const struct rel_relation **rlist=NULL;
1970 static rel_filter_t rfilter;
1971 static var_filter_t vfilter;
1972 real64 *derivatives;
1973 int32 *variables;
1974 mtx_coord_t coord;
1975 slv8_system_t sys;
1976
1977 /*
1978 * stop gcc whining about unused parameter
1979 */
1980 (void)vsta; (void)rhs; (void)esta; (void)n;
1981
1982 sys = (slv8_system_t)usrmem;
1983 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1984 rfilter.matchvalue =(REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1985
1986 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
1987 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
1988
1989 calc_J(sys);
1990 calc_residuals(sys);
1991 scale_system(sys);
1992
1993 for (offset = col = sys->J.reg.col.low;
1994 col <= sys->J.reg.col.high; col++) {
1995 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1996 nominal = sys->nominals.vec[col];
1997 low = var_lower_bound(var)/nominal;
1998 up = var_upper_bound(var)/nominal;
1999 /* KHACK: get rid of hard coded numbers */
2000 lower[col-offset] = low > -1e20 ? low : -1e20;
2001 upper[col-offset] = up < 1e20 ? up : 1e20;
2002 curr[col-offset] = sys->variables.vec[col]; /* already scaled */
2003 vsta[col-offset] = !var_nonbasic(var);
2004 }
2005 for (offset = row = sys->J.reg.row.low;
2006 row <= sys->J.reg.row.high; row++) {
2007 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2008 nominal = sys->weights.vec[row];
2009 fv[row-offset] = sys->residuals.vec[row]; /* already scaled */
2010 }
2011 for (row = 0; row < *m; row++) {
2012 type[row] = 0;
2013 }
2014 if (sys->obj != NULL) {
2015 type[*m-1] = 3;
2016 }
2017
2018 if (sys->obj != NULL) {
2019 len = rel_n_incidences(sys->obj);
2020 variables = (int32 *)ascmalloc(len*sizeof(int32));
2021 derivatives = (real64 *)ascmalloc(len*sizeof(real64));
2022 relman_diff2(sys->obj,&vfilter,derivatives,variables,
2023 &(obj_count),SAFE_CALC);
2024 }
2025
2026 count = count_old = 0;
2027
2028 colsta[0] = 1;
2029
2030 for (offset = col = sys->J.reg.col.low;
2031 col <= sys->J.reg.col.high; col++) {
2032 coord.col = col;
2033 var = sys->vlist[col];
2034 #if CONDBG
2035 if (!var_apply_filter(var,&vfilter) ) {
2036 FPRINTF(stderr,"var doesn't pass filter\n");
2037 }
2038 #endif /* CONDBG */
2039 len = var_n_incidences(var);
2040 rlist = var_incidence_list(var);
2041 count_old = count;
2042 for (c=0; c < len; c++) {
2043 /* assuming obj on list... check this */
2044 if (rel_apply_filter(rlist[c],&rfilter)) {
2045 coord.row = rel_sindex(rlist[c]);
2046 rowno[count] = C2F(rel_sindex(rlist[c]) - offset);
2047 value[count] = mtx_value(sys->J.mtx,&coord);
2048 nlflag[count] = 1; /* fix this later */
2049 if(rlist[c] == sys->obj) {
2050 #if CONDBG
2051 FPRINTF(stderr,"found objective in unexpected location\n");
2052 #endif /* CONDBG */
2053 }
2054 if (fabs(value[count]) > RTMAXJ) {
2055 #if CONDBG
2056 FPRINTF(stderr,"Large Jacobian value being set to RTMAXJ\n");
2057 #endif /* CONDBG */
2058 if (value[count] > 0) {
2059 value[count] = RTMAXJ-1;
2060 } else {
2061 value[count] = -RTMAXJ+1;
2062 }
2063 }
2064 count++;
2065 }
2066 if(rlist[c] == sys->obj) {
2067 for (r = 0; r < obj_count; r++) {
2068 if ( variables[r] == var_sindex(var) ) {
2069 rowno[count] = *m;
2070 value[count] = derivatives[r];
2071 nlflag[count] = 1; /* fix this later */
2072 if (fabs(value[count]) > RTMAXJ) {
2073 if (value[count] > 0) {
2074 value[count] = RTMAXJ-1;
2075 } else {
2076 value[count] = -RTMAXJ+1;
2077 }
2078 }
2079 count++;
2080 }
2081 }
2082 }
2083 }
2084 if (count_old != count) {
2085 colsta[col - offset] = C2F(count_old);
2086 }
2087 }
2088 colsta[F2C(*n1)] = *nz + 1;
2089 if (sys->obj != NULL) {
2090 ascfree(variables);
2091 ascfree(derivatives);
2092 }
2093 }
2094
2095 #if 0 /* not compatible with our version */
2096 /*
2097 * COIFBL Defines the nonlinearities of the model by returning
2098 * numerical values. It works on a block of rows during each call.
2099 * COIFBL( x, g, otn, nto, from, to, jac, stcl, rnum, cnum, nl, strw,
2100 * llen, indx, mode, errcnt, n, m, n1, m1, nz, usrmem)
2101 *
2102 * x - punt of evaluation provided by conopt
2103 * g - vector of function values
2104 * otn - old to new permutation vector
2105 * nto - new to old permutation vector
2106 * from - range in permutation
2107 * to - range in permutation
2108 * jac - vector of jacobian values.
2109 * The following are vectors defining the jacobian structure
2110 * stcl - start of column pointers
2111 * rnum - row numbers
2112 * cnum - column numbers
2113 * nl - nonlinearity flags
2114 * strw - start row pointers
2115 * llen - count of linear jacobian elements
2116 * indx - pointers from the row-wise representation
2117 * mode - indicator of mode of evaluation
2118 * errcnt- number of function evaluation errors
2119 * n - umber of variables
2120 * m - number of constraints
2121 * n1 - n+1
2122 * m1 - m+1
2123 * nz - number of jacobian elements
2124 * usrmem- user memory defined by conopt
2125 */
2126 static void slv8_coifbl(real64 *x, real64 *g, int32 *otn, int32 *nto,
2127 int32 *from, int32 *to, real64 *jac, int32 *stcl,
2128 int32 *rnum, int32 *cnum, int32 *nl, int32 *strw,
2129 int32 *llen, int32 *indx, int32 *mode, int32 *errcnt,
2130 int32 *n, int32 *m, int32 *n1, int32 *m1, int32 *nz,
2131 real64 *usrmem)
2132 {
2133 int32 offset, col, row, j, jj, len, c;
2134 real64 nominal, value;
2135 struct var_variable *var;
2136 struct rel_relation *rel;
2137 int32 *jac_row;
2138 int32 *variables;
2139 real64 *derivatives;
2140 static var_filter_t vfilter;
2141 slv8_system_t sys;
2142
2143 /*
2144 * stop gcc whining about unused parameter
2145 */
2146 (void)nto; (void)stcl; (void)rnum; (void)nl;
2147 (void)errcnt; (void)n1; (void)m1; (void)nz;
2148
2149 sys = (slv8_system_t)usrmem;
2150
2151 for (offset = col = sys->J.reg.col.low;
2152 col <= sys->J.reg.col.high; col++) {
2153 var = sys->vlist[col];
2154 nominal = sys->nominals.vec[col];
2155 value = x[col-offset] * nominal;
2156 var_set_value(var, value);
2157 }
2158 /* NOTE: could be more efficient when mode = 3 */
2159 if (*mode == 1 || *mode == 3) {
2160 for (offset = row = sys->J.reg.row.low;
2161 row <= sys->J.reg.row.high; row++) {
2162 if (F2C(*to) <= otn[row-offset] && otn[row-offset] <= F2C(*from)) {
2163 rel = sys->rlist[row];
2164 g[row-offset] = relman_eval(rel,&calc_ok,SAFE_CALC)
2165 * sys->weights.vec[row];
2166 }
2167 }
2168 if (F2C(*to) <= otn[F2C(*m)] && otn[F2C(*m)] <= F2C(*from)) {
2169 if(calc_objective(sys)){
2170 g[F2C(*m)] = sys->objective;
2171 } else {
2172 FPRINTF(MIF(sys),"slv8_coifbl: ERROR IN OBJECTIVE CALCULATION\n");
2173 }
2174 }
2175 }
2176 jac_row = (int32 *)ascmalloc((*n)*sizeof(int32));
2177 if (*mode == 2 || *mode == 3) {
2178 len = sys->con.maxrow;
2179 variables = (int32 *)ascmalloc(len*sizeof(int32));
2180 derivatives = (real64 *)ascmalloc(len*sizeof(real64));
2181 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
2182 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
2183 for (offset = row = sys->J.reg.row.low;
2184 row <= sys->J.reg.row.high; row++) {
2185 if (F2C(*to) <= otn[row-offset] && otn[row-offset] <= F2C(*from)) {
2186 rel = sys->rlist[row];
2187 relman_diff2(rel,&vfilter,derivatives,variables,
2188 &(len),SAFE_CALC);
2189 for (c = 0; c < len; c++) {
2190 jac_row[variables[c]] = derivatives[c];
2191 }
2192 for (j = strw[row-offset] + llen[row-offset];
2193 j < strw[row-offset + 1]; j++) {
2194 jj = indx[F2C(j)];
2195 jac[F2C(jj)] = jac_row[F2C(cnum[F2C(jj)])]
2196 * sys->weights.vec[row]
2197 * sys->nominals.vec[F2C(cnum[F2C(jj)])];
2198 if(fabs(jac[F2C(jj)]) > RTMAXJ) {
2199 if (jac[F2C(jj)] < 0) {
2200 jac[F2C(jj)] = -RTMAXJ+1;
2201 } else {
2202 jac[F2C(jj)] = RTMAXJ-1;
2203 }
2204 #if CONDBG
2205 FPRINTF(stderr,"large jac element\n");
2206 #endif /* CONDBG */
2207 }
2208 }
2209 }
2210 }
2211 }
2212 }
2213 #endif /* 0 */
2214
2215 /*
2216 * COIFDE Defines the nonlinearities of the model by returning
2217 * numerical values. It works on one row or equation at a time
2218 * COIFDE(x, g, jac, rowno, jcnm, mode, errcnt, newpt, n, nj, usrmem)
2219 *
2220 * x - punt of evaluation provided by conopt
2221 * g - function value
2222 * jac - jacobian values
2223 * rowno - number of the row for which nonlinearities will be eval
2224 * jcnm - list of column number fon the NL nonzeros
2225 * mode - indicator of mode of evaluation
2226 * errcnt - sum of number of func evaluation errors thus far
2227 * newpt - new point indicator
2228 * nj - number of nonlinear nonzero jacobian elements
2229 * n - number of variables
2230 * usrmem - user memory
2231 */
2232 static void slv8_coifde(real64 *x, real64 *g, real64 *jac, int32 *rowno,
2233 int32 *jcnm, int32 *mode, int32 *errcnt,
2234 int32 *newpt, int32 *n, int32 *nj, real64 *usrmem)
2235 {
2236 int32 offset, col, row, len, c;
2237 real64 nominal, value;
2238 struct var_variable *var;
2239 struct rel_relation *rel;
2240 int32 *variables;
2241 real64 *derivatives;
2242 static var_filter_t vfilter;
2243 slv8_system_t sys;
2244
2245 /*
2246 * stop gcc whining about unused parameter
2247 */
2248 (void)jcnm; (void)n; (void)nj;
2249
2250 sys = (slv8_system_t)usrmem;
2251 if (*newpt == 1) {
2252 for (offset = col = sys->J.reg.col.low;
2253 col <= sys->J.reg.col.high; col++) {
2254 var = sys->vlist[col];
2255 nominal = sys->nominals.vec[col];
2256 value = x[col-offset] * nominal;
2257 var_set_value(var, value);
2258 }
2259 }
2260 /* NOTE: could be more efficient when mode = 3
2261 * (with future versions of CONOPT)
2262 */
2263 if (*mode == 1 || *mode == 3) {
2264 offset = sys->J.reg.row.low;
2265 row = F2C(*rowno + offset);
2266 if ((*rowno == sys->con.m) && (sys->obj != NULL)){
2267 if(calc_objective(sys)){
2268 *g = sys->objective;
2269 } else {
2270 FPRINTF(MIF(sys),"slv8_coifde: ERROR IN OBJECTIVE CALCULATION\n");
2271 }
2272 } else {
2273 rel = sys->rlist[row];
2274 *g = relman_eval(rel,&calc_ok,SAFE_CALC)
2275 * sys->weights.vec[row];
2276 if (!calc_ok) {
2277 (*errcnt)++;
2278 }
2279 }
2280 }
2281 if (*mode == 2 || *mode == 3) {
2282 len = sys->con.maxrow;
2283 variables = (int32 *)ascmalloc(len*sizeof(int32));
2284 derivatives = (real64 *)ascmalloc(len*sizeof(real64));
2285 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
2286 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
2287
2288 offset = sys->J.reg.row.low;
2289 row = F2C(*rowno + offset);
2290 if ((*rowno == sys->con.m) && (sys->obj != NULL)){
2291 rel = sys->obj;
2292 calc_ok = relman_diff2(rel,&vfilter,derivatives,variables,
2293 &(len),SAFE_CALC);
2294 for (c = 0; c < len; c++) {
2295 jac[variables[c]] = derivatives[c]
2296 * sys->nominals.vec[variables[c]];
2297 }
2298 if (!calc_ok) {
2299 (*errcnt)++;
2300 }
2301 } else {
2302 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2303 calc_ok = relman_diff2(rel,&vfilter,derivatives,variables,
2304 &(len),SAFE_CALC);
2305 for (c = 0; c < len; c++) {
2306 jac[variables[c]] = derivatives[c]
2307 * sys->weights.vec[row] * sys->nominals.vec[variables[c]];
2308 }
2309 if (!calc_ok) {
2310 (*errcnt)++;
2311 }
2312 }
2313 for (c = 0; c < len; c++) {
2314 if(fabs(jac[variables[c]]) > RTMAXJ) {
2315 #if CONDBG
2316 FPRINTF(stderr,"large jac element\n");
2317 #endif /* CONDBG */
2318 if (jac[variables[c]] < 0) {
2319 jac[variables[c]] = -RTMAXJ+1;
2320 } else {
2321 jac[variables[c]] = RTMAXJ-1;
2322 }
2323 }
2324 }
2325 ascfree(variables);
2326 ascfree(derivatives);
2327 }
2328 }
2329
2330
2331 /*
2332 * COISTA Pass the solution from CONOPT to the modeler. It returns
2333 * completion status
2334 * COISTA(modsta, solsts, iter, objval, usrmem)
2335 *
2336 * modsta - model status
2337 * solsta - solver status
2338 * iter - number of iterations
2339 * objval - objective value
2340 * usrmem - user memory
2341 */
2342 static void slv8_coista(int32 *modsta, int32 *solsta, int32 *iter,
2343 real64 *objval, real64 *usrmem)
2344 {
2345 slv8_system_t sys;
2346 sys = (slv8_system_t)usrmem;
2347 sys->con.modsta = *modsta;
2348 sys->con.solsta = *solsta;
2349
2350 #if NONBASIC_DEBUG
2351 FPRINTF(ASCERR," Model Status = %d\n",*modsta);
2352 FPRINTF(ASCERR," Solver Status = %d\n",*solsta);
2353 #endif /* NONBASIC_DEBUG */
2354
2355 sys->con.iter = *iter;
2356 sys->con.obj = sys->objective = *objval;
2357 }
2358
2359
2360 /*
2361 * COIRS passes the solution from CONOPT to the modeler. It returns
2362 * solution values
2363 * COIRS(val, xmar, xbas, xsta, yval, ymar, ybas, ysta, n, m, usrmem)
2364 *
2365 * xval - the solution values of the variables
2366 * xmar - corresponding marginal values
2367 * xbas - basis indicator for column (at bound, basic, nonbasic)
2368 * xsta - status of column (normal, nonoptimal, infeasible,unbounded)
2369 * yval - values of the left hand side in all the rows
2370 * ymar - corresponding marginal values
2371 * ybas - basis indicator for row
2372 * ysta - status of row
2373 * n - number of variables
2374 * m - number of constraints
2375 * usrmem - user memory
2376 */
2377 static void slv8_coirs(real64 *xval, real64 *xmar, int32 *xbas, int32 *xsta,
2378 real64 *yval, real64 *ymar, int32 *ybas, int32 * ysta,
2379 int32 *n, int32 *m, real64 *usrmem)
2380 {
2381 int32 offset, col, c;
2382 real64 nominal, value;
2383 struct var_variable *var;
2384 slv8_system_t sys;
2385
2386 /*
2387 * stop gcc whining about unused parameter
2388 */
2389 (void)xmar; (void)xsta; (void)yval;
2390 (void)ymar; (void)ysta; (void)ybas; (void)m;
2391
2392
2393 sys = (slv8_system_t)usrmem;
2394
2395 offset = sys->J.reg.col.low;
2396 for (c = 0; c < *n; c++) {
2397 col = c + offset;
2398 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2399 nominal = sys->nominals.vec[col];
2400 /*
2401 * value of the variable
2402 */
2403 value = xval[c] * nominal;
2404 var_set_value(var, value);
2405 /*
2406 * status (basic, nonbasic) of the variable
2407 */
2408 if (xbas[c] != 2) {
2409 var_set_nonbasic(var,TRUE);
2410 #if NONBASIC_DEBUG
2411 FPRINTF(ASCERR," c = %d\n",c);
2412 FPRINTF(ASCERR," status = %d\n",xbas[c]);
2413 FPRINTF(ASCERR," nonbasic \n");
2414 #endif /* NONBASIC_DEBUG */
2415 } else {
2416 var_set_nonbasic(var,FALSE);
2417 }
2418 }
2419
2420 /* should pull out additional info here */
2421 }
2422
2423
2424 /*
2425 * COIUSZ communicates and update of an existing model to CONOPT
2426 * COIUSZ(nintg, ipsz, nreal, rpsz, usrmem)
2427 *
2428 * nintg - number of positions in ipsz
2429 * ipsz - array describing problem size and options
2430 * nreal - number of positions in rpsz
2431 * rpsz - array of reals describing problem size and options
2432 * usrmem- user memory
2433 */
2434 static void slv8_coiusz(int32 *nintg, int32 *ipsz, int32 *nreal,
2435 real64 *rpsz, real64 *usrmem)
2436 {
2437 /*
2438 * "zero changes" subroutines. the default for all the values is a
2439 * no change situation
2440 */
2441
2442 /*
2443 * stop gcc whining about unused parameter
2444 */
2445 (void)nintg; (void)ipsz; (void)nreal; (void)rpsz;
2446 (void)usrmem;
2447
2448 #if 0
2449 slv8_system_t sys;
2450
2451 /*
2452 * To Ken: This was in the subroutine before. But all the values
2453 * are the same as in coipsz. So, no redefintion is necessary since
2454 * the defaults contain the same information
2455 */
2456
2457 /*
2458 * stop gcc whining about unused parameter
2459 */
2460 (void)nintg; (void)nreal;
2461
2462 sys = (slv8_system_t)usrmem;
2463
2464 ipsz[F2C(1)] = sys->con.n;
2465 ipsz[F2C(2)] = sys->con.m;
2466 ipsz[F2C(3)] = sys->con.nz;
2467 ipsz[F2C(4)] = 0; /* FIX THESE AT A LATER DATE!!!! */
2468 ipsz[F2C(5)] = sys->con.nz; /* ASSUMING ALL NONLINEAR FOR NOW */
2469 if (sys->obj != NULL) {
2470 ipsz[F2C(6)] = relman_obj_direction(sys->obj);
2471 ipsz[F2C(7)] = sys->con.m; /* objective will be last row */
2472 } else {
2473 ipsz[F2C(7)] = 0;
2474 }
2475 ipsz[F2C(8)] = ITER_LIMIT;
2476 ipsz[F2C(12)] = 1; /* NON DEFAULT VALUE */
2477 ipsz[F2C(13)] = 1; /* NON DEFAULT VALUE */
2478
2479 rpsz[F2C(7)] = TIME_LIMIT;
2480
2481 #endif
2482
2483 return;
2484 }
2485
2486
2487 /*
2488 * COIOPT communicates non-default option values to CONOPT
2489 * COIOPT(name, rval, ival, lval, usrmem)
2490 * name - the name of a CONOPT CR-cell defined by the modeler
2491 * rval - the value to be assigned to name if the cells contains a real
2492 * ival - the value to be assigned to name if the cells contains an int
2493 * lval - the value to be assigned to name if the cells contains a log value
2494 * usrmem - user memory
2495 */
2496 static void slv8_coiopt(char *name, real64 *rval, int32 *ival,
2497 int32 *logical, real64 *usrmem)
2498 {
2499 slv8_system_t sys;
2500
2501 /*
2502 * stop gcc whining about unused parameter
2503 */
2504 (void)logical;
2505
2506 sys = (slv8_system_t)usrmem;
2507 name = memset(name,' ',8);
2508 while (sys->con.opt_count < slv8_PA_SIZE) {
2509 if (strlen(sys->p.parms[sys->con.opt_count].interface_label) == 6){
2510 if (strncmp(sys->p.parms[sys->con.opt_count].interface_label,
2511 "R",1) == 0) {
2512 name =
2513 strncpy(name, sys->p.parms[sys->con.opt_count].interface_label,6);
2514 *rval = sys->p.parms[sys->con.opt_count].info.r.value;
2515 sys->con.opt_count++;
2516 return;
2517 } else if (strncmp(sys->p.parms[sys->con.opt_count].interface_label,
2518 "L",1) == 0) {
2519 name =
2520 strncpy(name,sys->p.parms[sys->con.opt_count].interface_label,6);
2521 *ival = sys->p.parms[sys->con.opt_count].info.i.value;
2522 sys->con.opt_count++;
2523 return;
2524 }
2525 }
2526 sys->con.opt_count++;
2527 }
2528
2529 /* sending blank to quit iterative calling */
2530 name = memset(name,' ',8);
2531 }
2532
2533
2534 /*
2535 * COIPSZ communicates the model size and structure to CONOPT
2536 * COIPSZ(nintg, ipsz, nreal, rpsz, usrmem)
2537 *
2538 * ningt - number of positions in ipsz
2539 * ipsz - array describing problem size and options
2540 * nreal - number of positions in rpsz
2541 * rpsz - array of reals describing problem size and options
2542 * usrmem - user memory
2543 */
2544 static void slv8_coipsz(int32 *nintg, int32 *ipsz, int32 *nreal,
2545 real64 *rpsz, real64 *usrmem)
2546 {
2547 slv8_system_t sys;
2548
2549 /*
2550 * stop gcc whining about unused parameter
2551 */
2552 (void)nintg; (void)nreal;
2553
2554 sys = (slv8_system_t)usrmem;
2555
2556 ipsz[F2C(1)] = sys->con.n;
2557 ipsz[F2C(2)] = sys->con.m;
2558 ipsz[F2C(3)] = sys->con.nz;
2559 ipsz[F2C(4)] = 0; /* FIX THESE AT A LATER DATE!!!! */
2560 ipsz[F2C(5)] = sys->con.nz; /* ASSUMING ALL NONLINEAR FOR NOW */
2561 if (sys->obj != NULL) {
2562 ipsz[F2C(6)] = relman_obj_direction(sys->obj);
2563 ipsz[F2C(7)] = sys->con.m; /* objective will be last row */
2564 } else {
2565 ipsz[F2C(7)] = 0;
2566 }
2567 ipsz[F2C(8)] = ITER_LIMIT;
2568 ipsz[F2C(9)] = DOMLIM;
2569 ipsz[F2C(10)] = 1; /* OUTPUT TO SUBROUTINE */
2570 ipsz[F2C(11)] = 0; /* NO OUTPUT TO SCREEN */
2571 ipsz[F2C(12)] = 1; /* NON DEFAULT VALUE */
2572 ipsz[F2C(13)] = 1; /* NON DEFAULT VALUE */
2573 ipsz[F2C(14)] = 1; /* NON DEFAULT VALUE */
2574 ipsz[F2C(15)] = 1; /* NON DEFAULT VALUE */
2575 ipsz[F2C(16)] = 1; /* NON DEFAULT VALUE */
2576 ipsz[F2C(17)] = 0;
2577 ipsz[F2C(18)] = 0;
2578 ipsz[F2C(19)] = 0;
2579 ipsz[F2C(20)] = 0;
2580 ipsz[F2C(21)] = 0;
2581 ipsz[F2C(22)] = 1; /* NON DEFAULT VALUE */
2582 /*skipping remainder of ipsz which are fortran io parameters */
2583
2584 rpsz[F2C(1)] = 1e20;
2585 rpsz[F2C(2)] = -1e20;
2586 rpsz[F2C(3)] = 1.2e20;
2587 /*rpsz[F2C(4)] = NA*/
2588 /*rpsz[F2C(5)] = eps*/
2589 rpsz[F2C(6)] = 0;
2590 rpsz[F2C(7)] = TIME_LIMIT;
2591 rpsz[F2C(8)] = 1;
2592
2593 }
2594
2595 /* conopt communication subroutines */
2596 void slv8_coiec COIEC_ARGS {
2597 char *name=NULL;
2598 struct var_variable **vp;
2599 slv8_system_t sys;
2600 sys = (slv8_system_t)usrmem;
2601 vp=slv_get_solvers_var_list(SERVER);
2602 /* assumes cur = org */
2603 vp = vp + (*colno + sys->J.reg.col.low);
2604 name= var_make_name(SERVER,*vp);
2605 FPRINTF(stderr,"ERROR in variable:\n %s\n",name);
2606 FPRINTF(stdout," %.*s\n",*msglen,&msg[0]);
2607 if (name) {
2608 ascfree(name);
2609 }
2610 }
2611 void slv8_coier COIER_ARGS {
2612 char *name=NULL;
2613 struct rel_relation **rp;
2614 slv8_system_t sys;
2615 sys = (slv8_system_t)usrmem;
2616 rp=slv_get_solvers_rel_list(SERVER);
2617 rp = rp + (*rowno + sys->J.reg.row.low);;
2618 name= rel_make_name(SERVER,*rp);
2619 FPRINTF(stderr,"ERROR in relation:\n %s\n",name);
2620 FPRINTF(stdout," %.*s\n",*msglen,&msg[0]);
2621 if (name) {
2622 ascfree(name);
2623 }
2624 }
2625 void slv8_coienz COIENZ_ARGS {
2626 char *relname=NULL;
2627 char *varname=NULL;
2628 struct rel_relation **rp;
2629 struct var_variable **vp;
2630
2631 slv8_system_t sys;
2632 sys = (slv8_system_t)usrmem;
2633
2634 rp=slv_get_solvers_rel_list(SERVER);
2635 vp=slv_get_solvers_var_list(SERVER);
2636 /* assumes cur = org */
2637 rp = rp + (*rowno + sys->J.reg.row.low);;
2638 vp = vp + (*colno + sys->J.reg.col.low);
2639 relname= rel_make_name(SERVER,*rp);
2640 varname= var_make_name(SERVER,*vp);
2641 FPRINTF(stderr,"ERROR in jacobian element:\n");
2642 FPRINTF(stderr," relation: %s\n variable: %s\n",
2643 relname, varname);
2644 FPRINTF(stdout,"%.*s\n",*msglen,&msg[0]);
2645 if (relname) {
2646 ascfree(relname);
2647 }
2648 if (varname) {
2649 ascfree(varname);
2650 }
2651 }
2652
2653 void slv8_coimsg COIMSG_ARGS {
2654 int32 stop, i, len;
2655 char *line[15];
2656 /* should put option to make stop = *smsg or *nmsg
2657 * and option to route output
2658 */
2659 stop = *nmsg;
2660 for (i = 0; i < stop; i++) {
2661 len = llen[i];
2662 line[i] = &msgv[i*80];
2663 FPRINTF(stdout,"%.*s\n",len,line[i]);
2664 }
2665 }
2666
2667 void slv8_coiprg COIPRG_ARGS {
2668 slv8_system_t sys;
2669 sys = (slv8_system_t)usrmem;
2670 if (sys->con.progress_count == 0) {
2671 FPRINTF(stderr,
2672 " iter phase numinf numnop nsuper ");
2673 FPRINTF(stderr,
2674 " suminf objval rgmax\n");
2675 }
2676 FPRINTF(stdout,"%6i ",intrep[0]);
2677 FPRINTF(stdout," %6i ",intrep[1]);
2678 FPRINTF(stdout," %6i ",intrep[2]);
2679 FPRINTF(stdout," %6i ",intrep[3]);
2680 FPRINTF(stdout," %6i ",intrep[4]);
2681 FPRINTF(stdout," %16e ",rl[0]);
2682 FPRINTF(stdout," %16e ",rl[1]);
2683 FPRINTF(stdout," %7.2e ",rl[2]);
2684 FPRINTF(stdout,"\n");
2685 sys->con.progress_count++;
2686 if (sys->con.progress_count == 10) { /* 10 should be iface parm */
2687 sys->con.progress_count = 0;
2688 }
2689
2690 }
2691
2692 void slv8_coiorc COIORC_ARGS {
2693 if (*resid != 0.0) {
2694 char *relname=NULL;
2695 char *varname=NULL;
2696 struct rel_relation **rp;
2697 struct var_variable **vp;
2698 int32 row, col;
2699
2700 slv8_system_t sys;
2701 sys = (slv8_system_t)usrmem;
2702
2703 rp=slv_get_solvers_rel_list(SERVER);
2704 vp=slv_get_solvers_var_list(SERVER);
2705 /* assumes cur = org */
2706 row = F2C(*rowno + sys->J.reg.row.low);
2707 rp = rp + row;
2708 col = F2C(*colno + sys->J.reg.col.low);
2709 vp = vp + col;
2710 relname= rel_make_name(SERVER,*rp);
2711 varname= var_make_name(SERVER,*vp);
2712
2713 FPRINTF(stderr,"ERROR: Infeasible specification discovered at:\n");
2714 FPRINTF(stderr," relation: %s\n residual: %g\n",
2715 relname, (*resid)/sys->weights.vec[row]);
2716 FPRINTF(stderr," variable: %s\n value: %g\n",
2717 varname, (*value)*sys->nominals.vec[col]);
2718
2719 if (relname) {
2720 ascfree(relname);
2721 }
2722 if (varname) {
2723 ascfree(varname);
2724 }
2725 }
2726 }
2727 void slv8_coiscr COISCR_ARGS {
2728 FPRINTF(stdout,"%.*s\n",*len,&msg[0]);
2729 }
2730
2731
2732
2733 /*
2734 * slv_conopt iterate calls conopt_start, which calls coicsm
2735 * to starts CONOPT. The use of conopt_start is a H A C K to avoid
2736 * unresolved external during the linking of the CONOPT library.
2737 * See conopt.h
2738 */
2739 static void slv_conopt_iterate(slv8_system_t sys)
2740 {
2741 real64 **usrmem;
2742 conopt_pointers conopt_ptrs;
2743
2744 conopt_ptrs = (conopt_pointers)asccalloc
2745 (1, sizeof(struct conopt_function_pointers ) );
2746 conopt_ptrs->coirms_ptr = slv8_coirms;
2747 conopt_ptrs->coifbl_ptr = NULL;
2748 conopt_ptrs->coifde_ptr = slv8_coifde;
2749 conopt_ptrs->coirs_ptr = slv8_coirs;
2750 conopt_ptrs->coista_ptr = slv8_coista;
2751 conopt_ptrs->coiusz_ptr = NULL;
2752 conopt_ptrs->coiopt_ptr = slv8_coiopt;
2753 conopt_ptrs->coipsz_ptr = slv8_coipsz;
2754
2755 conopt_ptrs->coimsg_ptr = slv8_coimsg;
2756 conopt_ptrs->coiscr_ptr = slv8_coiscr;
2757 conopt_ptrs->coiec_ptr = slv8_coiec;
2758 conopt_ptrs->coier_ptr = slv8_coier;
2759 conopt_ptrs->coienz_ptr = slv8_coienz;
2760 conopt_ptrs->coiprg_ptr = slv8_coiprg;
2761 conopt_ptrs->coiorc_ptr = slv8_coiorc;
2762
2763 usrmem = (real64 **)(ascmalloc(2*sizeof(real64 *)));
2764
2765 /*
2766 * We pass the pointers to sys and conopt_ptrs instead of a usrmem array.
2767 * Cast the appropriate element of usrmem back to slv9_system_t and
2768 * conopt_pointers to access the information required
2769 */
2770 usrmem[0] = (real64 *)conopt_ptrs;
2771 usrmem[1] = (real64 *)sys;
2772
2773 sys->con.opt_count = 0; /* reset count on slv8_coiopt calls */
2774 sys->con.progress_count = 0; /* reset count on coiprg calls */
2775
2776 sys->con.kept = 1;
2777 conopt_start(&(sys->con.kept), usrmem, &(sys->con.lwork),
2778 sys->con.work, &(sys->con.maxusd), &(sys->con.curusd));
2779 /*
2780 * We need to check conopt's status codes before claim
2781 * optimization complete. For now just see status file
2782 */
2783 sys->con.optimized = 1;
2784
2785 ascfree(conopt_ptrs);
2786 ascfree(usrmem);
2787 }
2788
2789
2790 /*
2791 * Function created to provide the interface with the correct values
2792 * for number of iterations, residuals, solved variables, etc
2793 */
2794 static void update_block_information(slv8_system_t sys)
2795 {
2796 int32 row,col;
2797
2798 row = sys->J.reg.row.high - sys->J.reg.row.low + 1;
2799 col = sys->J.reg.col.high - sys->J.reg.col.low + 1;
2800 sys->s.block.current_size = MAX(row,col);
2801
2802 sys->s.block.iteration = sys->con.iter;
2803 sys->s.iteration = sys->con.iter;
2804
2805 if ( (sys->con.modsta == 1 || sys->con.modsta == 2)
2806 && sys->con.solsta == 1 ) {
2807 sys->s.converged = TRUE;
2808 sys->s.block.previous_total_size += sys->s.block.current_size;
2809 } else {
2810 if (sys->con.solsta == 2 ) {
2811 sys->s.converged = FALSE;
2812 sys->s.inconsistent = FALSE;
2813 } else {
2814 sys->s.converged = FALSE;
2815 sys->s.inconsistent = TRUE;
2816 }
2817 }
2818 }
2819
2820
2821 static void slv8_iterate(slv_system_t server, SlvClientToken asys)
2822 {
2823 slv8_system_t sys;
2824 FILE *mif;
2825 FILE *lif;
2826 sys = SLV8(asys);
2827 mif = MIF(sys);
2828 lif = LIF(sys);
2829 if (server == NULL || sys==NULL) return;
2830 if (check_system(SLV8(sys))) return;
2831 if( !sys->s.ready_to_solve ) {
2832 FPRINTF(mif,"ERROR: (slv8) slv8_iterate\n");
2833 FPRINTF(mif," Not ready to solve.\n");
2834 return;
2835 }
2836
2837 if (sys->s.block.current_block==-1) {
2838 conopt_initialize(sys);
2839 sys->s.converged = sys->con.optimized;
2840 update_status(sys);
2841 if( RELNOMSCALE == 1 || (strcmp(SCALEOPT,"RELNOM") == 0) ||
2842 (strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0) ){
2843 calc_relnoms(sys);
2844 }
2845 }
2846 if (sys->p.output.less_important && (sys->s.block.current_size >1 ||
2847 LIFDS)) {
2848 debug_delimiter(lif);
2849 }
2850 iteration_begins(sys);
2851 if (1 || sys->J.reg.row.high != sys->J.reg.col.high) {
2852 /*may have changed objective*/
2853 sys->obj = slv_get_obj_relation(server);
2854 slv_conopt_iterate(sys);
2855 update_block_information(sys); /* update values of block information */
2856 calc_objective(sys);
2857 calc_objectives(sys);
2858 sys->residuals.accurate = FALSE;
2859 calc_residuals(sys);
2860 update_cost(sys);
2861 iteration_ends(sys);
2862 update_status(sys);
2863 return;
2864 }
2865 }
2866
2867
2868 static void slv8_solve(slv_system_t server, SlvClientToken asys)
2869 {
2870 slv8_system_t sys;
2871 sys = SLV8(asys);
2872 if (server == NULL || sys==NULL) return;
2873 if (check_system(sys)) return;
2874 while( sys->s.ready_to_solve ) slv8_iterate(server,sys);
2875 }
2876
2877 static mtx_matrix_t slv8_get_jacobian(slv_system_t server, SlvClientToken sys)
2878 {
2879 if (server == NULL || sys==NULL) return NULL;
2880 if (check_system(SLV8(sys))) return NULL;
2881 return SLV8(sys)->J.mtx;
2882 }
2883
2884 static int32 slv8_destroy(slv_system_t server, SlvClientToken asys)
2885 {
2886 slv8_system_t sys;
2887
2888 /*
2889 * stop gcc whining about unused parameter
2890 */
2891 (void)server;
2892
2893 sys = SLV8(asys);
2894 if (check_system(sys)) return 1;
2895 destroy_vectors(sys);
2896 destroy_matrices(sys);
2897 slv_destroy_parms(&(sys->p));
2898 sys->integrity = DESTROYED;
2899 if (sys->s.cost) ascfree(sys->s.cost);
2900 if (sys->con.work != NULL) {
2901 ascfree(sys->con.work);
2902 sys->con.work = NULL;
2903 }
2904 ascfree( (POINTER)asys );
2905 return 0;
2906 }
2907
2908 int32 slv8_register(SlvFunctionsT *sft)
2909 {
2910 if (sft==NULL) {
2911 FPRINTF(stderr,"slv8_register called with NULL pointer\n");
2912 return 1;
2913 }
2914 #ifdef DYNAMIC_CONOPT
2915 if (conopt_load() == 1) {
2916 FPRINTF(stderr,"Registration failure: CONOPT dll unavailable\n");
2917 return 1;
2918 }
2919 #endif /* DYNAMIC_CONOPT */
2920
2921 sft->name = "CONOPT";
2922 sft->ccreate = slv8_create;
2923 sft->cdestroy = slv8_destroy;
2924 sft->celigible = slv8_eligible_solver;
2925 sft->getdefparam = slv8_get_default_parameters;
2926 sft->getparam = slv8_get_parameters;
2927 sft->setparam = slv8_set_parameters;
2928 sft->getstatus = slv8_get_status;
2929 sft->solve = slv8_solve;
2930 sft->presolve = slv8_presolve;
2931 sft->iterate = slv8_iterate;
2932 sft->resolve = slv8_resolve;
2933 sft->getlinsol = slv8_get_linsol_sys;
2934 sft->getlinsys = slv8_get_linsolqr_sys;
2935 sft->getsysmtx = slv8_get_jacobian;
2936 sft->dumpinternals = slv8_dump_internals;
2937 return 0;
2938 }
2939
2940 /* #endif */ /* #else clause of DYNAMIC_CONOPT */
2941 #endif /* #else clause of !STATIC_CONOPT && !DYNAMIC_CONOPT */

Properties

Name Value
svn:executable *

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