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

Properties

Name Value
svn:executable *

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