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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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