/[ascend]/trunk/models/johnpye/tron/asc_tron.c
ViewVC logotype

Contents of /trunk/models/johnpye/tron/asc_tron.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1503 - (show annotations) (download) (as text)
Sat Jun 23 12:44:10 2007 UTC (17 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 39970 byte(s)
Commenting re matrices
1 /* ASCEND modelling environment
2 Copyright (C) 2007 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//**
19 @file
20 Connection of the TRON solver into ASCEND.
21
22 THIS IS STILL VERY MUCH UNDER DEVELOPMENT AND INCOMPLETE. I'VE ACTUALLY
23 ONLY JUST STARTED WRITING IT by starting with asc_conopt.c and modifying.
24
25 The TRON solver is documented in SIAM J Optimiz., vol. 9, no. 4, by
26 Lin and Mor��. http://dx.doi.org/10.1137/S1052623498345075
27 *//*
28 ASCEND wrapper for TRON originally by John Pye, Jun 2007.
29 */
30
31 #include <math.h>
32 #include <ctype.h>
33
34 #include <solver/solver.h>
35 #include <system/slv_stdcalls.h>
36 #include <system/relman.h>
37
38 #include <utilities/ascPanic.h>
39 #include <utilities/ascMalloc.h>
40 #include <utilities/ascDynaLoad.h>
41 #include <utilities/mem.h>
42 #include <utilities/ascEnvVar.h>
43 #include <general/tm_time.h>
44 #include <general/env.h>
45
46 /*
47 'tron.h' declares what we expect to find in the DLL/SO that we will dlopen
48 when this solver is loaded.
49 */
50
51 #include "tron.h"
52
53 ASC_DLLSPEC SolverRegisterFn tron_register;
54
55 /*
56 Output in user defined TRON subroutines
57 */
58 #define TRON_DEBUG
59
60 #ifdef TRON_DEBUG
61 # define TRON_CONSOLE_DEBUG(...) CONSOLE_DEBUG(__VA_ARGS__)
62 #else
63 # define TRON_CONSOLE_DEBUG(...) (void)0
64 #endif
65
66 /*
67 makes lots of extra spew
68 */
69 #define DEBUG FALSE
70
71 #define TRONSYS(s) ((TronSystem *)(s))
72
73 /* for integrity checks */
74 #define TRON_CREATED ((int32)813029392)
75 #define TRON_DESTROYED ((int32)103289182)
76
77 /*------------------------------------------------------------------------------
78 DATA STRUCTURES
79 */
80
81 enum{
82 TRON_PARAM_ITERMAX
83 ,TRON_PARAM_MAXFEV
84 ,TRON_PARAM_FATOL
85 ,TRON_PARAM_FRTOL
86 ,TRON_PARAM_FMIN
87 ,TRON_PARAM_CGTOL
88 ,TRON_PARAM_GTOL
89 ,TRON_PARAM_TIMELIMIT
90 ,TRON_PARAM_SAFECALC
91 ,TRON_PARAMS
92 };
93
94 /**
95 TRON's matrix storage format, for a square n-by-n matrix
96 */
97 typedef struct{
98 double *ltp; /**< lower triangular part, compressed column storage, dimension nnz + n*p. */
99 double *diag; /**< diagonal elements, n elements */
100 int *col_ptr; /** pointers to the columns, dimension n+1. The nonzeros in column j of L are in the
101 col_ptr(j), ... , col_ptr(j+1) - 1 positions of l. */
102 int *row_ind; /** row indices for the strict lower
103 triangular part of L (ltp) in compressed column storage. */
104 } TronMatrix;
105
106 #define TRON_MATRIX_ARG(SYS,MAT) (SYS)->MAT.ltp,(SYS)->MAT.diag,(SYS)->MAT.col_ptr,(SYS)->MAT.row_ind
107
108 typedef struct{
109
110 char task[60];
111 slv_parameters_t params; /* Parameters */
112 slv_status_t status; /* Status (as of iteration end) */
113 slv_system_t slv; /* slv_system_t back-link */
114 struct slv_parameter pa[TRON_PARAMS]; /* dunno... */
115
116 /* Problem definition */
117 struct rel_relation *obj; /* Objective function: NULL = none */
118 struct rel_relation *old_obj;/* Objective function: NULL = none */
119 struct var_variable **vlist; /* Variable list (NULL terminated) */
120 struct rel_relation **rlist; /* Relation list (NULL terminated) */
121 long vtot; /* length of varlist */
122 long rtot; /* length of rellist */
123 int m, n; /* size of system matrix? */
124 double objective; /* value of objective function */
125 int vused; /* Free and incident variables */
126 int rused; /* Included relations */
127
128 /* Solver information */
129 int created; /* Whether the system been created */
130 int presolved; /* Whether the system been presolved */
131 int optimised;
132
133 double *x;
134 double *xupper; /* upper bounds on each variable */
135 double *xlower; /* lower bounds on each variable */
136 double *g; /* gradients of f with respect to each variable */
137 TronMatrix A; /* hessian matrix, passed to TRON? */
138 TronMatrix B; /* some matrix that we get back from TRON */
139 TronMatrix L; /* some other matrix that we get back from TRON */
140
141 double *xc, *s, *dsave, *wa;
142 int *indfree, *isave, *iwa;
143
144 double clock; /* CPU time */
145
146 } TronSystem;
147
148 /* this stuff was in the above as well */
149
150 #if 0
151 struct update_data update; /* Jacobian frequency counters */
152 int32 cap; /* Order of matrix/vectors */
153 int32 rank; /* Symbolic rank of problem */
154 #endif
155
156 #if 0
157 void *parm_array[TRON_PARAMS];
158 #endif
159
160 #if 0
161 /*
162 TRON DATA
163 */
164 struct tron_data con;
165 #endif
166
167 #if 0
168 /*
169 Calculated data (scaled)
170 */
171 struct jacobian_data J; /* linearized system */
172
173 struct vec_vector nominals; /* Variable nominals */
174 struct vec_vector weights; /* Relation weights */
175 struct vec_vector relnoms; /* Relation nominals */
176 struct vec_vector variables; /* Variable values */
177 struct vec_vector residuals; /* Relation residuals */
178
179 real64 objective; /* Objective function evaluation */
180 #endif
181
182 /*------------------*/
183
184 #if 0
185 /*
186 Auxiliary structures
187 */
188
189 struct update_data {
190 int jacobian; /* Countdown on jacobian updating */
191 int weights; /* Countdown on weights updating */
192 int nominals; /* Countdown on nominals updating */
193 int relnoms; /* Countdown on relnom updating */
194 int iterative; /* Countdown on iterative scale update */
195 };
196
197
198 /*
199 varpivots, relpivots used only in optimizing, if we rewrite calc_pivots
200 without them.
201 */
202 struct jacobian_data {
203 linsolqr_system_t sys; /* Linear system */
204 mtx_matrix_t mtx; /* Transpose gradient of residuals */
205 real64 *rhs; /* RHS from linear system */
206 unsigned *varpivots; /* Pivoted variables */
207 unsigned *relpivots; /* Pivoted relations */
208 unsigned *subregions; /* Set of subregion indeces */
209 dof_t *dofdata; /* dof data pointer from server */
210 mtx_region_t reg; /* Current block region */
211 int32 rank; /* Numerical rank of the jacobian */
212 enum factor_method fm; /* Linear factorization method */
213 boolean accurate; /* ? Recalculate matrix */
214 boolean singular; /* ? Can matrix be inverted */
215 boolean old_partition;/* old value of partition flag */
216 };
217 #endif
218
219 /* forward decls */
220
221 static void tron_initialize(TronSystem * sys);
222 static int check_system(TronSystem *sys);
223 static void iteration_begins(TronSystem *sys);
224 static void iteration_ends( TronSystem *sys);
225 static int32 tron_dof_changed(TronSystem *sys);
226 static void update_status( TronSystem *sys);
227 static boolean calc_objective( TronSystem *sys);
228 static boolean calc_residuals( TronSystem *sys);
229
230 /*------------------------------------------------------------------------------
231 SOLVER PARAMETERS
232 */
233
234 /*
235 these are the parameters that AMPL says it offers for the TRON solver:
236
237 itermax the limit on the number of conjugate gradient iterations.
238 maxfev the limit on the number of function evaluations
239 fatol the absolute error desired in the function
240 frtol the relative error desired in the function
241 fmin a lower bound for the function
242 cgtol convergence criteria for the conjugate gradient method
243 gtol relative 2-norm of projected grdients; convergence criteria for the trust region method
244 */
245
246 static
247 int32 tron_get_default_parameters(slv_system_t server, SlvClientToken asys
248 ,slv_parameters_t *parameters
249 ){
250 TronSystem *sys = NULL;
251 struct slv_parameter *new_parms = NULL;
252 int32 make_macros = 0;
253
254 if(server != NULL && asys != NULL) {
255 sys = TRONSYS(asys);
256 make_macros = 1;
257 }
258
259 if(parameters->parms == NULL) {
260 new_parms = ASC_NEW_ARRAY_OR_NULL(struct slv_parameter,TRON_PARAMS);
261 if(new_parms == NULL) {
262 return -1;
263 }
264 parameters->parms = new_parms;
265 parameters->dynamic_parms = 1;
266 }
267
268 parameters->num_parms = 0;
269 asc_assert(TRON_PARAMS==7);
270
271 slv_param_int(parameters,TRON_PARAM_ITERMAX
272 ,(SlvParameterInitInt){{"itermax"
273 ,"Max. CG iterations",1
274 ,"Maximum number of conjugate gradient iterations"
275 }, 30, 1, 20000}
276 );
277
278 slv_param_int(parameters,TRON_PARAM_MAXFEV
279 ,(SlvParameterInitInt){{"maxfev"
280 ,"Max. function evaluations",1
281 ,"Maximum number of function evaluations"
282 }, 30, 1, 20000}
283 );
284
285 slv_param_real(parameters,TRON_PARAM_FATOL
286 ,(SlvParameterInitReal){{"fatol"
287 ,"Error tolerance (absolute)",1
288 ,"The absolute error desired in the function"
289 }, 1e-8, 1e-13, 100000.5}
290 );
291
292 slv_param_real(parameters,TRON_PARAM_FRTOL
293 ,(SlvParameterInitReal){{"frtol"
294 ,"Error tolerance (relative)",1
295 ,"The relative error desired in the function"
296 }, 1e-8, 1e-13, 100000.5}
297 );
298
299 slv_param_real(parameters,TRON_PARAM_FMIN
300 ,(SlvParameterInitReal){{"fmin"
301 ,"Function lower bound",1
302 ,"A lower bound for the function"
303 }, -1e8, 1e8, -1e8}
304 );
305
306 slv_param_real(parameters,TRON_PARAM_CGTOL
307 ,(SlvParameterInitReal){{"cgtol"
308 ,"Tolerance for CG convergence criteria",1
309 ,"Convergence criteria for the conjugate gradient method"
310 }, 1e-8, 1e-13, 100000.5}
311 );
312
313 slv_param_real(parameters,TRON_PARAM_GTOL
314 ,(SlvParameterInitReal){{"gtol"
315 ,"Trust region tolerance",1
316 ,"relative 2-norm of projected gradients; convergence criteria for the trust region method"
317 }, 1e-8, 1e-13, 100000.5}
318 );
319
320 slv_param_int(parameters,TRON_PARAM_TIMELIMIT
321 ,(SlvParameterInitInt){{"timelimit"
322 ,"time limit",1
323 ,"time limit (seconds)"
324 }, 30, 1, 20000}
325 );
326
327 slv_param_bool(parameters,TRON_PARAM_SAFECALC
328 ,(SlvParameterInitBool){{"safecalc"
329 ,"safe calc",1
330 ,"safe floating-point operations"
331 }, 0}
332 );
333
334 asc_assert(parameters->num_parms==TRON_PARAMS);
335
336 return 1;
337 }
338
339 static void tron_get_parameters(slv_system_t server, SlvClientToken asys
340 , slv_parameters_t *parameters
341 ){
342 TronSystem *sys;
343 UNUSED_PARAMETER(server);
344
345 sys = TRONSYS(asys);
346 if (check_system(sys)) return;
347 mem_copy_cast(&(sys->params),parameters,sizeof(slv_parameters_t));
348 }
349
350 static void tron_set_parameters(slv_system_t server, SlvClientToken asys
351 ,slv_parameters_t *parameters
352 ){
353 TronSystem *sys;
354 UNUSED_PARAMETER(server);
355
356 sys = TRONSYS(asys);
357 if (check_system(sys)) return;
358 mem_copy_cast(parameters,&(sys->params),sizeof(slv_parameters_t));
359 }
360
361 /*-----------------------------------------------------------------------------
362 SET-UP AND DESTROY SOLVER DATA
363 */
364
365 static SlvClientToken tron_create(slv_system_t server, int32*statusindex){
366 TronSystem *sys;
367
368 sys = ASC_NEW_CLEAR(TronSystem);
369 if(sys==NULL){
370 *statusindex = 1;
371 return sys;
372 }
373
374 sys->slv = server;
375 sys->params.parms = sys->pa;
376 sys->params.dynamic_parms = 0;
377 tron_get_default_parameters(server,(SlvClientToken)sys,&(sys->params));
378
379 sys->created = TRON_CREATED;
380 sys->presolved = 0;
381
382 sys->params.whose = (*statusindex);
383
384 sys->status.ok = TRUE;
385 sys->status.calc_ok = TRUE;
386 sys->status.costsize = 0;
387 sys->status.cost = NULL; /*redundant, but sanity-preserving */
388 sys->vlist = slv_get_solvers_var_list(server);
389 sys->rlist = slv_get_solvers_rel_list(server);
390 sys->rtot = slv_get_num_solvers_rels(server);
391 sys->vtot = slv_get_num_solvers_vars(server);
392
393 sys->obj = slv_get_obj_relation(server);
394
395 if (sys->vlist == NULL) {
396 ascfree(sys);
397 ERROR_REPORTER_HERE(ASC_PROG_ERR,"TRON called with no variables.");
398 *statusindex = -2;
399 return NULL; /* prolly leak here */
400 }
401 if (sys->rlist == NULL && sys->obj == NULL) {
402 ascfree(sys);
403 ERROR_REPORTER_HERE(ASC_PROG_ERR,"TRON called with no relations or objective.");
404 *statusindex = -1;
405 return NULL; /* prolly leak here */
406 }
407
408 /*
409 apparently we don't give a damn about the objective list or the pars or
410 bounds or extrels or any of the other crap.
411 */
412
413 slv_check_var_initialization(server);
414
415 sys->xc = ASC_NEW_ARRAY(double,sys->n);
416 sys->s = ASC_NEW_ARRAY(double,sys->n);
417 sys->indfree = ASC_NEW_ARRAY(int,sys->n);
418 sys->isave = ASC_NEW_ARRAY(int,sys->n*3);
419 sys->dsave = ASC_NEW_ARRAY(double,sys->n*3);
420 sys->wa = ASC_NEW_ARRAY(double,sys->n*7);
421 sys->iwa = ASC_NEW_ARRAY(int,sys->n*3);
422
423 *statusindex = 0; /* <-- what is that? */
424
425 return((SlvClientToken)sys);
426 }
427
428
429 /**
430 @TODO document this
431 */
432 static int32 tron_destroy(slv_system_t server, SlvClientToken asys){
433 TronSystem *sys;
434 UNUSED_PARAMETER(server);
435
436 sys = TRONSYS(asys);
437 if (check_system(sys)) return 1;
438 /* destroy_vectors(sys); */
439 /* destroy_matrices(sys); */
440 slv_destroy_parms(&(sys->params));
441 sys->created = TRON_DESTROYED;
442 if(sys->status.cost){
443 ASC_FREE(sys->status.cost);
444 }
445
446 /* clear all the work arrays used by TRON */
447 #define F(AA) if(sys->AA)ASC_FREE(sys->AA)
448 F(xc); F(s); F(indfree); F(isave); F(dsave); F(wa); F(iwa);
449 #undef F
450
451 ascfree( (POINTER)asys );
452 return 0;
453 }
454
455 #if 0
456 static void create_matrices(slv_system_t server, TronSystem *sys){
457 sys->J.mtx = mtx_create();
458 mtx_set_order(sys->J.mtx,sys->cap);
459 structural_analysis(server,sys);
460 }
461
462
463 static void create_vectors(TronSystem *sys){
464 sys->nominals.vec = create_array(sys->cap,real64);
465 sys->nominals.rng = &(sys->J.reg.col);
466 sys->weights.vec = create_array(sys->cap,real64);
467 sys->weights.rng = &(sys->J.reg.row);
468 sys->relnoms.vec = create_array(sys->cap,real64);
469 sys->relnoms.rng = &(sys->J.reg.row);
470 sys->variables.vec = create_array(sys->cap,real64);
471 sys->variables.rng = &(sys->J.reg.col);
472 sys->residuals.vec = create_array(sys->cap,real64);
473 sys->residuals.rng = &(sys->J.reg.row);
474 }
475 #endif
476
477 #if 0
478 static void destroy_matrices( TronSystem *sys){
479 if( sys->J.mtx ) {
480 mtx_destroy(sys->J.mtx);
481 }
482 }
483
484 static void destroy_vectors( TronSystem *sys){
485 destroy_array(sys->nominals.vec);
486 destroy_array(sys->weights.vec);
487 destroy_array(sys->relnoms.vec);
488 destroy_array(sys->variables.vec);
489 destroy_array(sys->residuals.vec);
490 }
491 #endif
492
493 /*------------------------------------------------------------------------------
494 HIGH-LEVEL SOLVE ROUTINES
495 */
496
497 static int tron_presolve(slv_system_t server, SlvClientToken asys){
498 struct rel_relation **rp;
499 int32 ind;
500 int32 matrix_creation_needed = 1;
501 TronSystem *sys;
502
503 #if 0
504 int32 cap;
505 struct var_variable **vp;
506 int *cntvect, temp;
507 #endif
508
509 TRON_CONSOLE_DEBUG("PRESOLVE");
510
511 sys = TRONSYS(asys);
512 iteration_begins(sys);
513 check_system(sys);
514 if( sys->vlist == NULL ) {
515 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Variable list was not set.");
516 return -1;
517 }
518 if( sys->rlist == NULL && sys->obj == NULL ) {
519 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Relation list and objective were not set.");
520 return -2;
521 }
522
523 sys->obj = slv_get_obj_relation(server); /*may have changed objective*/
524
525 if(!sys->obj){
526 ERROR_REPORTER_HERE(ASC_PROG_ERR,"No objective function was specified");
527 return -3;
528 }
529
530 if(sys->presolved > 0) { /* system has been presolved before */
531 if(!tron_dof_changed(sys) /*no changes in fixed or included flags*/
532 /* && sys->params.partition == sys->J.old_partition
533 && sys->obj == sys->old_obj */
534 ){
535 matrix_creation_needed = 0;
536 TRON_CONSOLE_DEBUG("YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION");
537 }
538 }
539
540 rp=sys->rlist;
541 for( ind = 0; ind < sys->rtot; ++ind ) {
542 rel_set_satisfied(rp[ind],FALSE);
543 }
544 if( matrix_creation_needed ) {
545
546 #if 0
547 cap = slv_get_num_solvers_rels(sys->slv);
548 sys->cap = slv_get_num_solvers_vars(sys->slv);
549 sys->cap = MAX(sys->cap,cap);
550 vp=sys->vlist;
551 for( ind = 0; ind < sys->vtot; ++ind ) {
552 var_set_in_block(vp[ind],FALSE);
553 }
554 rp=sys->rlist;
555 for( ind = 0; ind < sys->rtot; ++ind ) {
556 rel_set_in_block(rp[ind],FALSE);
557 rel_set_satisfied(rp[ind],FALSE);
558 }
559 #endif
560
561 sys->presolved = 1; /* full presolve recognized here */
562
563 #if 0
564 sys->J.old_partition = sys->params.partition;
565 #endif
566 sys->old_obj = sys->obj;
567
568 slv_sort_rels_and_vars(server,&(sys->m),&(sys->n));
569 TRON_CONSOLE_DEBUG("FOUND %d CONSTRAINTS AND %d VARS",sys->m,sys->n);
570 if (sys->obj != NULL) {
571 TRON_CONSOLE_DEBUG("ADDING OBJECT AS A ROW");
572 sys->m++; /* treat objective as a row */
573 }
574
575 /* destroy_vectors(sys); */
576 /* destroy_matrices(sys); */
577 /* create_matrices(server,sys); */
578 /* create_vectors(sys); */
579
580 sys->status.block.current_reordered_block = -2;
581 }
582
583 /* Reset status */
584 sys->optimised = 0;
585 sys->status.iteration = 0;
586 sys->status.cpu_elapsed = 0.0;
587 sys->status.converged = sys->status.diverged = sys->status.inconsistent = FALSE;
588 sys->status.block.previous_total_size = 0;
589 sys->status.costsize = 1+sys->status.block.number_of;
590
591 if( matrix_creation_needed ) {
592 #if 0
593 destroy_array(sys->status.cost);
594 sys->status.cost = create_zero_array(sys->status.costsize,struct slv_block_cost);
595 for( ind = 0; ind < sys->status.costsize; ++ind ) {
596 sys->status.cost[ind].reorder_method = -1;
597 }
598 } else {
599 reset_cost(sys->status.cost,sys->status.costsize);
600 #endif
601 }
602
603 /* set to go to first unconverged block */
604 sys->status.block.current_block = -1;
605 sys->status.block.current_size = 0;
606 sys->status.calc_ok = TRUE;
607 sys->status.block.iteration = 0;
608 sys->objective = MAXDOUBLE/2000.0;
609
610 update_status(sys);
611 iteration_ends(sys);
612 sys->status.cost[sys->status.block.number_of].time=sys->status.cpu_elapsed;
613
614 return 0;
615 }
616
617 const char TRON_START[] = "START";
618 const char TRON_GH[] = "GH";
619 const char TRON_CONV[] = "CONV";
620
621 /**
622 @TODO document this
623 */
624 static int tron_iterate(slv_system_t server, SlvClientToken asys){
625 TronSystem *sys;
626 sys = TRONSYS(asys);
627
628 if(server == NULL || sys==NULL) return -1;
629 if(check_system(sys)) return -2;
630 if(!sys->status.ready_to_solve){
631 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Not ready to solve.");
632 return 1;
633 }
634
635 if(sys->status.block.current_block==-1) {
636 tron_initialize(sys);
637 sys->status.converged = sys->optimised;
638 update_status(sys);
639 /* do something with scaling here? calc_relnoms(sys)? */
640 }
641
642 iteration_begins(sys);
643
644 /* may have changed objective (how does that happen, BTW?) */
645 sys->obj = slv_get_obj_relation(server);
646
647 if(sys->task[0]=='F' || strncmp(sys->task,TRON_START,5)==0){
648
649 // Evaluate the function at x and store in f.
650
651 }
652 if(strncmp(sys->task,TRON_GH,2)==0 || strncmp(sys->task,TRON_START,5)==0){
653
654 // Evaluate the gradient at x and store in g.
655
656 // Evaluate the Hessian at x and store in compressed
657 // column storage in (a,adiag,acol_ptr,arow_ind)
658 }
659
660 double frtol = SLV_PARAM_REAL(&(sys->params),TRON_PARAM_FRTOL);
661 double fatol = SLV_PARAM_REAL(&(sys->params),TRON_PARAM_FATOL);
662 double cgtol = SLV_PARAM_REAL(&(sys->params),TRON_PARAM_CGTOL);
663 int itermax = SLV_PARAM_INT(&(sys->params),TRON_PARAM_ITERMAX);
664 double fmin = SLV_PARAM_REAL(&(sys->params),TRON_PARAM_FMIN);
665 double delta = SLV_PARAM_REAL(&(sys->params),TRON_PARAM_GTOL);
666 /** @TODO fmin should be taken from the model declaration somehow, not a solar parameter. */
667
668 DTRON(&(sys->n),sys->x,sys->xlower,sys->xupper,&(sys->objective),sys->g
669 ,TRON_MATRIX_ARG(sys,A)
670 ,&frtol,&fatol,&fmin,&cgtol,&itermax,&delta,sys->task
671 ,TRON_MATRIX_ARG(sys,B)
672 ,TRON_MATRIX_ARG(sys,L)
673 ,sys->xc,sys->s,sys->indfree
674 ,sys->isave,sys->dsave,sys->wa,sys->iwa
675 );
676
677 if(strncmp(sys->task,TRON_CONV,4)==0){
678 sys->status.converged = 1;
679 TRON_CONSOLE_DEBUG("System has converged");
680 }
681
682 TRON_CONSOLE_DEBUG("DTRON status is '%s'",sys->task);
683
684 calc_objective(sys);
685 #if 0
686 calc_objectives(sys);
687 sys->residuals.accurate = FALSE;
688 #endif
689 calc_residuals(sys);
690 #if 0
691 update_cost(sys);
692 #endif
693 iteration_ends(sys);
694 update_status(sys);
695
696 return 0;
697 }
698
699 /**
700 @TODO document this
701 */
702 static int tron_solve(slv_system_t server, SlvClientToken asys){
703 TronSystem *sys;
704 int err = 0;
705 sys = TRONSYS(asys);
706 if(server == NULL || sys==NULL) return -1;
707 if(check_system(sys)) return -2;
708 while(sys->status.ready_to_solve)err = err | tron_iterate(server,sys);
709 return err;
710 }
711
712 /*------------------------------------------------------------------------------
713 INTEGRITY CHECKS
714 */
715
716 /**
717 Checks sys for NULL and for integrity.
718 */
719 static int check_system(TronSystem *sys){
720
721 if( sys == NULL ) {
722 ERROR_REPORTER_HERE(ASC_PROG_ERR,"NULL system handle");
723 return 1;
724 }
725
726 switch( sys->created ) {
727 case TRON_CREATED:
728 return 0;
729 case TRON_DESTROYED:
730 ERROR_REPORTER_HERE(ASC_PROG_ERR,"System was recently destroyed.");
731 return 1;
732 default:
733 ERROR_REPORTER_HERE(ASC_PROG_ERR,"System reused or never allocated.");
734 return 1;
735 }
736 }
737
738 /*------------------------------------------------------------------------------
739 CALCULATION ROUTINES
740
741 ok = calc_objective(sys)
742 ok = calc_residuals(sys)
743 ok = calc_J(sys)
744 calc_nominals(sys)
745 calc_weights(sys)
746 scale_J(sys)
747 scale_variables(sys)
748 scale_residuals(sys)
749 */
750
751 #if 0
752 /**
753 Count jacobian elements and set max to the number of elements
754 in the densest row
755 */
756 static int32 num_jacobian_nonzeros(TronSystem *sys, int32 *max){
757 int32 row, len, licn,c,count,row_max;
758 struct rel_relation *rel;
759 rel_filter_t rf;
760 var_filter_t vf;
761 const struct var_variable **list;
762
763 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
764 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
765 vf.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
766 vf.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
767
768 licn = 0;
769 *max = 0;
770 row_max = sys->m;
771 if (sys->obj != NULL) {
772 row_max--;
773 }
774 /* replace at leisure with
775 * relman_jacobian_count(sys->rlist,row_max,&vfilter,&rfilter,max);
776 */
777 for( row = 0; row < row_max; row++ ) {
778 rel = sys->rlist[row];
779 if (rel_apply_filter(rel,&rf)) { /* shouldn't be needed */
780 len = rel_n_incidences(rel);
781 list = rel_incidence_list(rel);
782 count = 0;
783 for (c=0; c < len; c++) {
784 if( var_apply_filter(list[c],&vf) ) {
785 licn++;
786 count++;
787 }
788 }
789 *max = MAX(*max,count);
790 }
791 }
792 if (sys->obj != NULL) {
793 rel = sys->obj;
794 len = rel_n_incidences(rel);
795 list = rel_incidence_list(rel);
796 count = 0;
797 for (c=0; c < len; c++) {
798 if( var_apply_filter(list[c],&vf) ) {
799 licn++;
800 count++;
801 }
802 }
803 *max = MAX(*max,count);
804 }
805 return licn;
806 }
807 #endif
808
809 /**
810 Evaluate the objective function.
811 */
812 static boolean calc_objective(TronSystem *sys){
813 boolean calc_ok = TRUE;
814 asc_assert(sys->obj!=NULL);
815 if(sys->obj){
816 sys->objective = relman_eval(sys->obj,&calc_ok,SLV_PARAM_BOOL(&(sys->params),TRON_PARAM_SAFECALC));
817 }else{
818 sys->objective = 0.0;
819 }
820 return calc_ok;
821 }
822
823 #if 0
824 /**
825 Evaluate all objectives.
826 */
827 static boolean calc_objectives( TronSystem *sys){
828 int32 len,i;
829 static rel_filter_t rfilter;
830 struct rel_relation **rlist=NULL;
831 rfilter.matchbits = (REL_INCLUDED);
832 rfilter.matchvalue =(REL_INCLUDED);
833 rlist = slv_get_solvers_obj_list(sys->slv);
834 len = slv_get_num_solvers_objs(sys->slv);
835 calc_ok = TRUE;
836 for (i = 0; i < len; i++) {
837 if (rel_apply_filter(rlist[i],&rfilter)) {
838 asc_assert(rlist[i]!=NULL);
839 relman_eval(rlist[i],&calc_ok,SAFE_CALC);
840 if DEBUG
841 if (calc_ok == FALSE) {
842 ERROR_REPORTER_HERE(ASC_PROG_ERR,"error in calc_objectives");
843 calc_ok = TRUE;
844 }
845 endif /* DEBUG */
846 }
847 }
848 return calc_ok;
849 }
850 #endif
851
852 #if 0
853 /**
854 Calculate all of the residuals in the current block and compute
855 the residual norm for block status.
856
857 @return true iff calculations preceded without error.
858 */
859 static boolean calc_residuals( TronSystem *sys){
860 int32 row;
861 struct rel_relation *rel;
862 double time0;
863
864 if( sys->residuals.accurate ) return TRUE;
865
866 calc_ok = TRUE;
867 row = sys->residuals.rng->low;
868 time0=tm_cpu_time();
869 for( ; row <= sys->residuals.rng->high; row++ ) {
870 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
871 if DEBUG
872 if (!rel) {
873 int32r;
874 r=mtx_row_to_org(sys->J.mtx,row);
875 ERROR_REPORTER_HERE(ASC_PROG_ERR
876 ,"NULL relation found at row %d rel %d in calc_residuals!",(int)row,r
877 );
878 }
879 endif /* DEBUG */
880 asc_assert(rel!=NULL);
881 sys->residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
882
883 relman_calc_satisfied(rel,sys->params.tolerance.feasible);
884 }
885 sys->status.block.functime += (tm_cpu_time() -time0);
886 sys->status.block.funcs++;
887 square_norm( &(sys->residuals) );
888 sys->status.block.residual = calc_sqrt_D0(sys->residuals.norm2);
889 if(!calc_ok){
890 TRON_CONSOLE_DEBUG("ERROR IN EVALUATION");
891 }
892 return(calc_ok);
893 }
894 #endif
895
896 #if 0
897 /**
898 Calculate the current block of the jacobian.
899 It is initially unscaled.
900 */
901 static boolean calc_J( TronSystem *sys){
902 int32 row;
903 var_filter_t vfilter;
904 double time0;
905 real64 resid;
906
907 calc_ok = TRUE;
908 vfilter.matchbits = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
909 vfilter.matchvalue = (VAR_ACTIVE | VAR_INCIDENT | VAR_SVAR);
910 time0=tm_cpu_time();
911 mtx_clear_region(sys->J.mtx,&(sys->J.reg));
912 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
913 struct rel_relation *rel;
914 rel = sys->rlist[row];
915 relman_diffs(rel,&vfilter,sys->J.mtx,&resid,SAFE_CALC);
916 }
917 sys->status.block.jactime += (tm_cpu_time() - time0);
918 sys->status.block.jacs++;
919
920 if( --(sys->update.nominals) <= 0 ) sys->nominals.accurate = FALSE;
921 if( --(sys->update.weights) <= 0 ) sys->weights.accurate = FALSE;
922
923 return(calc_ok);
924 }
925 #endif
926
927 #if 0
928 /**
929 Retrieve the nominal values of all of the block variables,
930 and ensure that they are all strictly positive.
931 */
932 static void calc_nominals( TronSystem *sys){
933 int32 col;
934 if( sys->nominals.accurate ) return;
935 col = sys->nominals.rng->low;
936 if(strcmp(SCALEOPT,"NONE") == 0 ||
937 strcmp(SCALEOPT,"ITERATIVE") == 0){
938 for( ; col <= sys->nominals.rng->high; col++ ) {
939 sys->nominals.vec[col] = 1;
940 }
941 } else {
942 for( ; col <= sys->nominals.rng->high; col++ ) {
943 struct var_variable *var;
944 real64 n;
945 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
946 n = var_nominal(var);
947 if( n <= 0.0 ) {
948 if( n == 0.0 ) {
949 n = TOO_SMALL;
950
951 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
952 FPRINTF(ASCERR,"Variable ");
953 print_var_name(ASCERR,sys,var);
954 FPRINTF(ASCERR," has nominal value of zero. Resetting to %g.\n",n);
955 error_reporter_end_flush();
956
957 var_set_nominal(var,n);
958 } else {
959 n = -n;
960
961 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
962 FPRINTF(ASCERR,"Variable ");
963 print_var_name(ASCERR,sys,var);
964 FPRINTF(ASCERR," has negative nominal value. Resetting to %g.\n",n);
965 error_reporter_end_flush();
966
967 var_set_nominal(var,n);
968 }
969 }
970 if DEBUG
971 ERROR_REPORTER_START_HERE(ASC_USER_ERROR);
972 FPRINTF(ASCERR,"Column %d is");
973 print_var_name(ASCERR,sys,var);
974 FPRINTF(ASCERR,"\nScaling of column %d is %g\n",col,n);
975 error_reporter_end_flush();
976 endif /* DEBUG */
977 sys->nominals.vec[col] = n;
978 }
979 }
980 square_norm( &(sys->nominals) );
981 sys->update.nominals = UPDATE_NOMINALS;
982 sys->nominals.accurate = TRUE;
983 }
984 #endif
985
986 #if 0
987 /**
988 Calculate the weights of all of the block relations
989 to scale the rows of the Jacobian.
990 */
991 static void calc_weights( TronSystem *sys)
992 {
993 mtx_coord_t nz;
994 real64 sum;
995
996 if( sys->weights.accurate )
997 return;
998
999 nz.row = sys->weights.rng->low;
1000 if(strcmp(SCALEOPT,"NONE") == 0 ||
1001 strcmp(SCALEOPT,"ITERATIVE") == 0) {
1002 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
1003 sys->weights.vec[nz.row] = 1;
1004 }
1005 } else if (strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
1006 strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0) {
1007 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
1008 sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
1009 sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
1010 }
1011 } else if (strcmp(SCALEOPT,"RELNOM") == 0 ||
1012 strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0) {
1013 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
1014 sys->weights.vec[nz.row] =
1015 1.0/rel_nominal(sys->rlist[mtx_row_to_org(sys->J.mtx,nz.row)]);
1016 }
1017 }
1018 square_norm( &(sys->weights) );
1019 sys->update.weights = UPDATE_WEIGHTS;
1020 sys->residuals.accurate = FALSE;
1021 sys->weights.accurate = TRUE;
1022 }
1023 #endif
1024
1025
1026 /**
1027 Reset all flags to setup a new solve.
1028 Should set sys->status.block.current_block = -1
1029 before calling.
1030
1031 @TODO This is currently a HACK! Not sure if should call when done.
1032 */
1033 void tron_initialize( TronSystem *sys){
1034
1035 sys->status.block.current_block++;
1036 /*
1037 * Next line was added to create the aray cost, whis is used by
1038 * the interface to display residuals and number of iterations
1039 */
1040 sys->status.costsize = 1+sys->status.block.number_of;
1041
1042 if( sys->status.block.current_block < sys->status.block.number_of ) {
1043 boolean ok;
1044
1045 sys->status.block.iteration = 0;
1046 sys->status.block.cpu_elapsed = 0.0;
1047 sys->status.block.functime = 0.0;
1048 sys->status.block.jactime = 0.0;
1049 sys->status.block.funcs = 0;
1050 sys->status.block.jacs = 0;
1051
1052 sys->status.calc_ok = TRUE;
1053
1054 if( !(ok = calc_objective(sys)) ) {
1055 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Objective calculation errors detected.");
1056 }
1057
1058 #if 0
1059 sys->status.calc_ok = sys->status.calc_ok && ok;
1060
1061 if(!(sys->params.ignore_bounds) ) {
1062 slv_ensure_bounds(
1063 sys->slv, sys->J.reg.col.low,
1064 sys->J.reg.col.high,MIF(sys)
1065 );
1066 }
1067
1068 sys->residuals.accurate = FALSE;
1069 if( !(ok = calc_residuals(sys)) ) {
1070 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Residual calculation errors detected.");
1071 }
1072
1073 sys->status.calc_ok = sys->status.calc_ok && ok;
1074
1075 /* Must be updated as soon as required */
1076 sys->J.accurate = FALSE;
1077 sys->update.weights = 0;
1078 sys->update.nominals = 0;
1079 sys->update.relnoms = 0;
1080 sys->update.iterative = 0;
1081 sys->variables.accurate = FALSE;
1082 #endif
1083 }
1084 }
1085
1086
1087 /*------------------------------------------------------------------------------
1088 ITERATION BEGIN/END ROUTINES
1089 */
1090
1091 /**
1092 Prepare sys for entering an iteration, increasing the iteration counts
1093 and starting the clock.
1094 */
1095 static void iteration_begins(TronSystem *sys){
1096 sys->clock = tm_cpu_time();
1097 ++(sys->status.block.iteration);
1098 ++(sys->status.iteration);
1099 }
1100
1101
1102 /*
1103 Prepare sys for exiting an iteration, stopping the clock and recording
1104 the cpu time.
1105 */
1106 static void iteration_ends( TronSystem *sys){
1107 double cpu_elapsed; /* elapsed this iteration */
1108
1109 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
1110 sys->status.block.cpu_elapsed += cpu_elapsed;
1111 sys->status.cpu_elapsed += cpu_elapsed;
1112 }
1113
1114
1115 /**
1116 Update the solver status.
1117 */
1118 static void update_status( TronSystem *sys){
1119 boolean unsuccessful;
1120
1121 #if 0
1122 if( !sys->status.converged ) {
1123 sys->status.time_limit_exceeded =
1124 (sys->status.block.cpu_elapsed >= TIME_LIMIT);
1125 sys->status.iteration_limit_exceeded =
1126 (sys->status.block.iteration >= ITER_LIMIT);
1127 }
1128 #endif
1129
1130 unsuccessful = sys->status.diverged || sys->status.inconsistent ||
1131 sys->status.iteration_limit_exceeded || sys->status.time_limit_exceeded;
1132
1133 sys->status.ready_to_solve = !unsuccessful && !sys->status.converged;
1134 sys->status.ok = !unsuccessful && sys->status.calc_ok && !sys->status.struct_singular;
1135 }
1136
1137 /*-----------------------------------------------------------------------------
1138 EXTERNAL ROUTINES (see slv_client.h)
1139 */
1140
1141 static int32 tron_eligible_solver(slv_system_t server){
1142 struct rel_relation **rp;
1143 for( rp=slv_get_solvers_rel_list(server); *rp != NULL ; ++rp ) {
1144 if( rel_less(*rp) || rel_greater(*rp) ){
1145 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"less-than and greater-than relations are not permitted with TRON");
1146 return(FALSE);
1147 }
1148 }
1149 return(TRUE);
1150 }
1151
1152
1153 static int tron_get_status(slv_system_t server, SlvClientToken asys
1154 ,slv_status_t *status
1155 ){
1156 TronSystem *sys;
1157 UNUSED_PARAMETER(server);
1158
1159 sys = TRONSYS(asys);
1160 if (check_system(sys)) return 1;
1161 mem_copy_cast(&(sys->status),status,sizeof(slv_status_t));
1162 return 0;
1163 }
1164
1165 #if 0
1166 /**
1167 Perform structural analysis on the system, setting the flags in
1168 status.
1169
1170 The problem must be set up, the relation/variable list
1171 must be non-NULL. The jacobian (linear) system must be created
1172 and have the correct order (stored in sys->cap). Everything else
1173 will be determined here.
1174
1175 On entry there isn't yet a correspondence between var_sindex and
1176 jacobian column. Here we establish that.
1177
1178 @NOTE this function has been striped of its guts for TRON and may go away
1179 */
1180 static void structural_analysis(slv_system_t server, TronSystem *sys){
1181
1182 var_filter_t vfilter;
1183 rel_filter_t rfilter;
1184
1185 /*
1186 * The server has marked incidence flags already.
1187 */
1188 /* count included equalities */
1189 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1190 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
1191 sys->rused = slv_count_solvers_rels(server,&rfilter);
1192
1193 /* count free and incident vars */
1194 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
1195 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
1196 sys->vused = slv_count_solvers_vars(server,&vfilter);
1197
1198 /* Symbolic analysis */
1199 sys->rtot = slv_get_num_solvers_rels(server);
1200 sys->vtot = slv_get_num_solvers_vars(server);
1201
1202 /*
1203 * The next few lines are used to calculate the rank of the nonlinear
1204 * system. We need it to evaluate if the system is structurally
1205 * singular or not. Calculating this number will keep TRON from
1206 * displaying a "structurally singular" error message
1207 */
1208 if (sys->rtot) {
1209 slv_block_partition(server);
1210 }
1211 sys->J.dofdata = slv_get_dofdata(server);
1212 sys->rank = sys->J.dofdata->structural_rank;
1213 /*
1214 * Unify the partitions since we feed TRON with a single block.
1215 */
1216 slv_block_unify(server);
1217
1218
1219 sys->J.reg.row.low = sys->J.reg.col.low = 0;
1220 sys->J.reg.row.high = sys->m - 1;
1221 if (sys->obj != NULL) sys->J.reg.row.high--;
1222 sys->J.reg.col.high = sys->n - 1;
1223
1224 if(slv_check_bounds(sys->slv,sys->vused,-1,"fixed ")){
1225 sys->status.inconsistent = 1;
1226 }
1227
1228 /* Initialize Status */
1229 sys->status.over_defined = (sys->rused > sys->vused);
1230 sys->status.under_defined = (sys->rused < sys->vused);
1231 sys->status.struct_singular = (sys->rank < sys->rused);
1232 sys->status.block.number_of = (slv_get_solvers_blocks(sys->slv))->nblocks;
1233
1234 }
1235 #endif
1236
1237 /**
1238 Check if any fixed or included flags have
1239 changed since the last presolve.
1240 */
1241 static int32 tron_dof_changed(TronSystem *sys){
1242 int32 ind, result = 0;
1243 /* Currently we have two copies of the fixed and included flags
1244 which must be kept in sync. The var_fixed and rel_included
1245 functions perform the syncronization and hence must be called
1246 over the whole var list and rel list respectively. When we move
1247 to using only one set of flags (bit flags) this function can
1248 be changed to return 1 at the first indication of a change
1249 in the dof. */
1250
1251 /* search for vars that were fixed and are now free */
1252 for( ind = sys->vused; ind < sys->vtot; ++ind ) {
1253 if( !var_fixed(sys->vlist[ind]) && var_active(sys->vlist[ind]) ) {
1254 ++result;
1255 }
1256 }
1257 /* search for rels that were unincluded and are now included */
1258 for( ind = sys->rused; ind < sys->rtot; ++ind ) {
1259 if( rel_included(sys->rlist[ind]) && rel_active(sys->rlist[ind])) {
1260 ++result;
1261 }
1262 }
1263 /* search for vars that were free and are now fixed */
1264 for( ind = sys->vused -1; ind >= 0; --ind ) {
1265 if( var_fixed(sys->vlist[ind]) || !var_active(sys->vlist[ind])) {
1266 ++result;
1267 }
1268 }
1269 /* search for rels that were included and are now unincluded */
1270 for( ind = sys->rused -1; ind >= 0; --ind ) {
1271 if( !rel_included(sys->rlist[ind]) || !rel_active(sys->rlist[ind]) ) {
1272 ++result;
1273 }
1274 }
1275 return result;
1276 }
1277
1278 #if 0
1279 static void reset_cost(struct slv_block_cost *cost,int32 costsize){
1280 int32 ind;
1281 for( ind = 0; ind < costsize; ++ind ) {
1282 cost[ind].size = 0;
1283 cost[ind].iterations = 0;
1284 cost[ind].funcs = 0;
1285 cost[ind].jacs = 0;
1286 cost[ind].functime = 0;
1287 cost[ind].jactime = 0;
1288 cost[ind].time = 0;
1289 cost[ind].resid = 0;
1290 }
1291 }
1292 #endif
1293
1294 #if 0
1295 /**
1296 Update the values of the array cost, which is used by the interface
1297 to display residual and number of iterations. For use after running CONOPT
1298 */
1299 static void update_cost(TronSystem *sys)
1300 {
1301 int32 ci;
1302 if (sys->status.cost == NULL) {
1303 sys->status.cost = create_zero_array(sys->status.costsize,struct slv_block_cost);
1304 } else {
1305 reset_cost(sys->status.cost,sys->status.costsize);
1306 }
1307 ci=sys->status.block.current_block;
1308 sys->status.cost[ci].size = sys->status.block.current_size;
1309 sys->status.cost[ci].iterations = sys->status.block.iteration;
1310 sys->status.cost[ci].resid = sys->status.block.residual;
1311 }
1312 #endif
1313
1314 /*------------------------------------------------------------------------------
1315 LOADING TRON FROM SHARED LIBRARY, IF AVAILABLE
1316 */
1317
1318 typedef struct{
1319 dtron_fn_t *dtron_ptr;
1320 } DTronFns;
1321
1322 DTronFns tron_fptrs;
1323
1324 int tron_dlopen(){
1325 static int loaded=0;
1326 char *libpath;
1327 int status;
1328 char fnsymbol[400], *c;
1329 const char *libname=ASC_TRON_LIB;
1330 const char *envvar;
1331
1332 if(loaded) {
1333 return 0; /* already loaded */
1334 }
1335
1336 CONSOLE_DEBUG("LOADING TRON...");
1337
1338 envvar = ASC_TRON_ENVVAR;
1339
1340 /* need to import this variable into the ascend 'environment' */
1341 if(-1!=env_import(ASC_TRON_ENVVAR,getenv,Asc_PutEnv)){
1342 CONSOLE_DEBUG("Searching in path '%s' (from env var '%s')",getenv(envvar),envvar);
1343 }
1344 libpath = SearchArchiveLibraryPath(libname, ASC_TRON_DLPATH, envvar);
1345
1346 if(libpath==NULL){
1347 ERROR_REPORTER_NOLINE(ASC_PROG_ERR
1348 , "Library '%s' could not be located (check value of env var '%s' and/or default path '%s')"
1349 , libname, envvar, ASC_TRON_DLPATH
1350 );
1351 return 1;
1352 }
1353
1354 status = Asc_DynamicLoad(libpath, NULL);
1355 if (status != 0) {
1356 ASC_FREE(libpath);
1357 return 1; /* failed to load */
1358 }
1359
1360
1361 # if defined(FNAME_UCASE_NODECOR) || defined(FNAME_UCASE_DECOR) || defined(FNAME_UCASE_PREDECOR)
1362 # define FNCASE(C) C=toupper(C)
1363 # elif defined(FNAME_LCASE_NODECOR) || defined(FNAME_LCASE_DECOR)
1364 # define FNCASE(C) C=tolower(C)
1365 # else
1366 # error "Need platform-specific information (asc_tron.c)"
1367 # endif
1368
1369 # if defined(FNAME_UCASE_DECOR) || defined(FNAME_LCASE_DECOR)
1370 # define FNDECOR(S,L) strcat(S,"_")
1371 # elif defined(FNAME_UCASE_PREDECOR) /* on windows, precede with _ and append @L (integer value of L) */
1372 # define FNDECOR(S,L) strcat(S,L);for(c=S+strlen(S)+1;c>S;--c){*c=*(c-1);} *S='_';
1373 # else
1374 # define FNDECOR(S,L) (void)0
1375 # endif
1376
1377 # define FN_PTR_GET(T,A,V,L) \
1378 sprintf(fnsymbol,"%s",#T); \
1379 for(c=fnsymbol;*c!='\0';++c){ \
1380 FNCASE(*c); \
1381 } \
1382 FNDECOR(fnsymbol,L); \
1383 tron_fptrs.T##_ptr = (T##_fn_t *)Asc_DynamicFunction(libpath,fnsymbol); \
1384 if(tron_fptrs.T##_ptr==NULL)status+=1;
1385
1386 FN_PTR_GET(dtron,TRON_DTRON_ARGS,TRON_DTRON_VALS,"");
1387
1388 # undef FN_PTR_GET
1389 # undef FNDECOR
1390 # undef FNCASE
1391
1392 ASC_FREE(libpath);
1393
1394 if(status!=0){
1395 return 1; /* failed to resolve all symbols */
1396 }
1397
1398 loaded = 1; /* save result for next time, as we will never unload it */
1399 return 0;
1400 }
1401
1402 /*------------------------------------------------------------------------------
1403 REGISTERING 'TRON' WITH ASCEND
1404 */
1405
1406 static const SlvFunctionsT tron_internals = {
1407 10
1408 ,"TRON"
1409 ,tron_create
1410 ,tron_destroy
1411 ,tron_eligible_solver
1412 ,tron_get_default_parameters
1413 ,tron_get_parameters
1414 ,tron_set_parameters
1415 ,tron_get_status
1416 ,tron_solve
1417 ,tron_presolve
1418 ,tron_iterate
1419 ,NULL
1420 ,NULL
1421 ,NULL
1422 ,NULL
1423 };
1424
1425 int tron_register(void){
1426 #ifndef ASC_LINKED_TRON
1427 if(tron_dlopen()){
1428 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Failed to load TRON");
1429 return 1;
1430 }
1431 #endif
1432 return solver_register(&tron_internals);
1433 }

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