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

Properties

Name Value
svn:executable *

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