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

Properties

Name Value
svn:executable *

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