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