/[ascend]/trunk/base/generic/solver/slv8.c
ViewVC logotype

Contents of /trunk/base/generic/solver/slv8.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:executable *

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