/[ascend]/trunk/solvers/conopt/asc_conopt.c
ViewVC logotype

Contents of /trunk/solvers/conopt/asc_conopt.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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