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

Properties

Name Value
svn:executable *

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