/[ascend]/trunk/ipslv-temp/slv5.westdir_notruncate.c
ViewVC logotype

Contents of /trunk/ipslv-temp/slv5.westdir_notruncate.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 916 - (show annotations) (download) (as text)
Tue Oct 31 02:39:16 2006 UTC (16 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 129140 byte(s)
Adding Vicente's IPSlv code (needs revising and testing against the new code)\
1 /*
2 * IPSlv ASCEND Interior Point Method Solver
3 * by Vicente Rico-Ramirez based on QRSlv
4 * Created:
5 * Version: $ $
6 * Version control file: $ $
7 * Date last modified: $ $
8 * Last modified by: $Author: $
9 *
10 * This file is part of the SLV solver.
11 *
12 * The SLV solver is free software; you can redistribute
13 * it and/or modify it under the terms of the GNU General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
16 *
17 * The SLV solver is distributed in hope that it will be
18 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with the program; if not, write to the Free Software Foundation,
24 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
25 * COPYING. COPYING is found in ../compiler.
26 *
27 */
28
29 #include <math.h>
30 #include <stdarg.h>
31 #include "utilities/ascConfig.h"
32 #include "utilities/ascSignal.h"
33 #include "utilities/ascMalloc.h"
34 #include "utilities/set.h"
35 #include "general/time.h"
36 #include "utilities/mem.h"
37 #include "utilities/ascPanic.h"
38 #include "general/list.h"
39 #include "compiler/fractions.h"
40 #include "compiler/dimen.h"
41 #include "compiler/functype.h"
42 #include "compiler/func.h"
43 #include "solver/mtx.h"
44 #include "solver/linsol.h"
45 #include "solver/linsolqr.h"
46 #include "solver/slv_types.h"
47 #include "solver/var.h"
48 #include "solver/rel.h"
49 #include "solver/discrete.h"
50 #include "solver/conditional.h"
51 #include "solver/logrel.h"
52 #include "solver/bnd.h"
53 #include "solver/calc.h"
54 #include "solver/relman.h"
55 #include "solver/slv_common.h"
56 #include "solver/slv_client.h"
57 #include "solver/slv5.h"
58 #include "solver/slv_stdcalls.h"
59
60 #if !defined(STATIC_IPSLV) && !defined(DYNAMIC_IPSLV)
61 int slv5_register(SlvFunctionsT *f)
62 {
63 (void)f; /* stop gcc whine about unused parameter */
64
65 FPRINTF(stderr,"IPSlv not compiled in this ASCEND IV.\n");
66 return 1;
67 }
68 #else /* either STATIC_IPSLV or DYNAMIC_IPSLV is defined */
69 #ifdef DYNAMIC_IPSLV
70 /* do dynamic loading stuff. yeah, right */
71 #else /* following is used if STATIC_IPSLV is defined */
72
73 #define DEBUG FALSE
74 #define DEBUG_ITERATION TRUE
75 #define DEBUG_COMPLEMENTARY_VAR TRUE
76 #define DEBUG_CENTERING TRUE
77 #define SLV5(s) ((slv5_system_t)(s))
78 #define SERVER (sys->slv)
79 #define slv5_PA_SIZE 47 /* MUST INCREMENT WHEN ADDING PARAMETERS */
80 #define slv5_RA_SIZE 11
81
82 /* do not delete (or extend) this array definition.
83 */
84 #define IEX(n) slv5_iaexpln[(n)]
85 #define slv5_IA_SIZE 16
86 static char *slv5_iaexpln[slv5_IA_SIZE] = {
87 "If lifds != 0 and showlessimportant is TRUE, show direct solve details",
88 "If savlin != 0, write out matrix data file at each iteration to SlvLinsol.dat",
89 "scale residuals by relation nominals for evaluating progress",
90 "Cutoff is the block size cutoff for MODEL-based reordering of partitions",
91 "Update jacobian every this many major iterations",
92 "Update row scalings every this many major iterations",
93 "Update column scalings every this many major iterations",
94 "Require misunderstood reduction somewhere in the stepping algorithm",
95 "Require residual >= some other number in the stepping algorithm",
96 "Check jacobian for poorly scaled columns and whine if found",
97 "Truncate whole step vector rather than componentwise at variable bound",
98 "Reorder option. 0 = MODEL based, 1 = MODEL based2, 2 = simple spk1",
99 "Use safe calculation routines",
100 "Update relation nominal scalings every this many major iterations",
101 "Max iterations for iterative scaling",
102 "scaleopt = 0: 2norm,= 1: relnom,= 2 2norm + iterative,= 3: relnom + iterative,= 4: iterative"
103 };
104
105 /* change slv5_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
106 #define UPDATE_JACOBIAN_PTR (sys->parm_array[0])
107 #define UPDATE_JACOBIAN ((*(int *)UPDATE_JACOBIAN_PTR))
108 #define UPDATE_WEIGHTS_PTR (sys->parm_array[1])
109 #define UPDATE_WEIGHTS ((*(int *)UPDATE_WEIGHTS_PTR))
110 #define UPDATE_NOMINALS_PTR (sys->parm_array[2])
111 #define UPDATE_NOMINALS ((*(int *)UPDATE_NOMINALS_PTR))
112 #define REDUCE_PTR (sys->parm_array[3])
113 #define REDUCE ((*(int *)REDUCE_PTR))
114 #define EXACT_LINE_SEARCH_PTR (sys->parm_array[4])
115 #define EXACT_LINE_SEARCH ((*(int *)EXACT_LINE_SEARCH_PTR))
116 #define DUMPCNORM_PTR (sys->parm_array[5])
117 #define DUMPCNORM ((*(int *)DUMPCNORM_PTR))
118 #define TRUNCATE_PTR (sys->parm_array[6])
119 #define TRUNCATE ((*(int *)TRUNCATE_PTR))
120 #define SAFE_CALC_PTR (sys->parm_array[7])
121 #define SAFE_CALC ((*(int *)SAFE_CALC_PTR))
122 #define SCALEOPT_PTR (sys->parm_array[8])
123 #define SCALEOPT ((*(char **)SCALEOPT_PTR))
124 #define UPDATE_RELNOMS_PTR (sys->parm_array[9])
125 #define UPDATE_RELNOMS ((*(int *)UPDATE_RELNOMS_PTR))
126 #define ITSCALELIM_PTR (sys->parm_array[10])
127 #define ITSCALELIM ((*(int *)ITSCALELIM_PTR))
128 #define RELNOMSCALE_PTR (sys->parm_array[11])
129 #define RELNOMSCALE ((*(int *)RELNOMSCALE_PTR))
130 #define TOO_SMALL_PTR (sys->parm_array[12])
131 #define TOO_SMALL ((*(real64 *)TOO_SMALL_PTR))
132 #define CNLOW_PTR (sys->parm_array[13])
133 #define CNLOW ((*(real64 *)CNLOW_PTR))
134 #define CNHIGH_PTR (sys->parm_array[14])
135 #define CNHIGH ((*(real64 *)CNHIGH_PTR))
136 #define TOWARD_BOUNDS_PTR (sys->parm_array[15])
137 #define TOWARD_BOUNDS ((*(real64 *)TOWARD_BOUNDS_PTR))
138 #define POSITIVE_DEFINITE_PTR (sys->parm_array[16])
139 #define POSITIVE_DEFINITE ((*(real64 *)POSITIVE_DEFINITE_PTR))
140 #define DETZERO_PTR (sys->parm_array[17])
141 #define DETZERO ((*(real64 *)DETZERO_PTR))
142 #define STEPSIZEERR_MAX_PTR (sys->parm_array[18])
143 #define STEPSIZEERR_MAX ((*(real64 *)STEPSIZEERR_MAX_PTR))
144 #define PARMRNG_MIN_PTR (sys->parm_array[19])
145 #define PARMRNG_MIN ((*(real64 *)PARMRNG_MIN_PTR))
146 #define MIN_COEF_PTR (sys->parm_array[20])
147 #define MIN_COEF ((*(real64 *)MIN_COEF_PTR))
148 #define MAX_COEF_PTR (sys->parm_array[21])
149 #define MAX_COEF ((*(real64 *)MAX_COEF_PTR))
150 #define ITSCALETOL_PTR (sys->parm_array[22])
151 #define ITSCALETOL ((*(real64 *)ITSCALETOL_PTR))
152 #define IGNORE_BOUNDS_PTR (sys->parm_array[23])
153 #define IGNORE_BOUNDS ((*(int32 *)IGNORE_BOUNDS_PTR))
154 #define SHOW_MORE_IMPT_PTR (sys->parm_array[24])
155 #define SHOW_MORE_IMPT ((*(int32 *)SHOW_MORE_IMPT_PTR))
156 #define RHO_PTR (sys->parm_array[25])
157 #define RHO ((*(real64 *)RHO_PTR))
158 #define PARTITION_PTR (sys->parm_array[26])
159 #define PARTITION ((*(int32 *)PARTITION_PTR))
160 #define SHOW_LESS_IMPT_PTR (sys->parm_array[27])
161 #define SHOW_LESS_IMPT ((*(int32 *)SHOW_LESS_IMPT_PTR))
162 #define AUTO_RESOLVE_PTR (sys->parm_array[28])
163 #define AUTO_RESOLVE ((*(int32 *)AUTO_RESOLVE_PTR))
164 #define TIME_LIMIT_PTR (sys->parm_array[29])
165 #define TIME_LIMIT ((*(int32 *)TIME_LIMIT_PTR))
166 #define ITER_LIMIT_PTR (sys->parm_array[30])
167 #define ITER_LIMIT ((*(int32 *)ITER_LIMIT_PTR))
168 #define STAT_TOL_PTR (sys->parm_array[31])
169 #define STAT_TOL ((*(real64 *)STAT_TOL_PTR))
170 #define TERM_TOL_PTR (sys->parm_array[32])
171 #define TERM_TOL ((*(real64 *)TERM_TOL_PTR))
172 #define SING_TOL_PTR (sys->parm_array[33])
173 #define SING_TOL ((*(real64 *)SING_TOL_PTR))
174 #define PIVOT_TOL_PTR (sys->parm_array[34])
175 #define PIVOT_TOL ((*(real64 *)PIVOT_TOL_PTR))
176 #define FEAS_TOL_PTR (sys->parm_array[35])
177 #define FEAS_TOL ((*(real64 *)FEAS_TOL_PTR))
178 #define LIFDS_PTR (sys->parm_array[36])
179 #define LIFDS ((*(int32 *)LIFDS_PTR))
180 #define SAVLIN_PTR (sys->parm_array[37])
181 #define SAVLIN ((*(int32 *)SAVLIN_PTR))
182 #define REORDER_OPTION_PTR (sys->parm_array[38])
183 #define REORDER_OPTION ((*(char **)REORDER_OPTION_PTR))
184 #define CUTOFF_PTR (sys->parm_array[39])
185 #define CUTOFF ((*(int32 *)CUTOFF_PTR))
186 #define FACTOR_OPTION_PTR (sys->parm_array[40])
187 #define FACTOR_OPTION ((*(char **)FACTOR_OPTION_PTR))
188 #define CONVOPT_PTR (sys->parm_array[41])
189 #define CONVOPT ((*(char **)CONVOPT_PTR))
190 #define LINTIME_PTR (sys->parm_array[42])
191 #define LINTIME ((*(int *)LINTIME_PTR))
192 #define SIGMA_CONSTANT_PTR (sys->parm_array[43])
193 #define SIGMA_CONSTANT ((*(int32 *)SIGMA_CONSTANT_PTR))
194 #define SIGMA_PTR (sys->parm_array[44])
195 #define SIGMA ((*(real64 *)SIGMA_PTR))
196 #define MU_MAXIMUM_PTR (sys->parm_array[45])
197 #define MU_MAXIMUM ((*(int32 *)MU_MAXIMUM_PTR))
198 #define MUK_PTR (sys->parm_array[46])
199 #define MUK ((*(real64 *)MUK_PTR))
200 /* change slv5_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
201
202
203 #define REX(n) slv5_raexpln[(n)]
204 static
205 char *slv5_raexpln[slv5_RA_SIZE] = {
206 "Var nominal to use if user specifies 0.0",
207 "Smallest column norm we won't complain about if checking",
208 "Largest column norm we won't complain about if checking",
209 "If bound is in the way, we go this fraction toward it",
210 "Hessian fudge number when optimizing",
211 "Minimum 2x2 determinant of newton/gradient we consider non-parallel",
212 "Step size must be determined this precisely, or prngmin happy",
213 "Parameter range must be this narrow to exit inner loop if step size unhappy",
214 "'Largest' drop in maxstep allowed",
215 "'Smallest' drop in maxstep allowed",
216 "scale termination ratio for iterative method"};
217
218 struct update_data {
219 int jacobian; /* Countdown on jacobian updating */
220 int weights; /* Countdown on weights updating */
221 int nominals; /* Countdown on nominals updating */
222 int relnoms; /* Countdown on relnom updating */
223 int iterative; /* Countdown on iterative scale update */
224 };
225
226
227 struct jacobian_data {
228 linsolqr_system_t sys; /* Linear system */
229 mtx_matrix_t mtx; /* Transpose gradient of residuals */
230 real64 *rhs; /* RHS from linear system */
231 dof_t *dofdata; /* dof data pointer from server */
232 mtx_region_t reg; /* Current block region */
233 int32 rank; /* Numerical rank of the jacobian */
234 enum factor_method fm; /* Linear factorization method */
235 boolean accurate; /* ? Recalculate matrix */
236 boolean singular; /* ? Can matrix be inverted */
237 boolean old_partition;/* old value of partition flag */
238 };
239
240 struct slv5_system_structure {
241
242 /*
243 * Problem definition
244 */
245 slv_system_t slv; /* slv_system_t back-link */
246 struct rel_relation *obj; /* Objective function: NULL = none */
247 struct var_variable **vlist; /* Variable list (NULL terminated) */
248 struct rel_relation **rlist; /* Relation list (NULL terminated) */
249
250 /*
251 * Solver information
252 */
253 int integrity; /* ? Has the system been created */
254 int32 presolved; /* ? Has the system been presolved */
255 slv_parameters_t p; /* Parameters */
256 slv_status_t s; /* Status (as of iteration end) */
257 struct update_data update; /* Jacobian frequency counters */
258 int32 cap; /* Order of matrix/vectors */
259 int32 rank; /* Symbolic rank of problem */
260 int32 vused; /* Free and incident variables */
261 int32 vtot; /* length of varlist */
262 int32 rused; /* Included relations */
263 int32 rtot; /* length of rellist */
264 double clock; /* CPU time */
265 void *parm_array[slv5_PA_SIZE]; /* array of pointers to param values */
266 struct slv_parameter pa[slv5_PA_SIZE];/* &pa[0] => sys->p.parms */
267
268 /*
269 * Calculated data (scaled)
270 */
271 struct jacobian_data J; /* linearized system */
272 struct vector_data nominals; /* Variable nominals */
273 struct vector_data weights; /* Relation weights */
274 struct vector_data relnoms; /* Relation nominals */
275 struct vector_data variables; /* Variable values */
276 struct vector_data residuals; /* Relation residuals */
277 struct vector_data newton_residuals; /* Newton Relation residuals */
278 struct vector_data perturbed_residuals;/* Perturbed Relation residuals */
279 struct vector_data gamma; /* Feasibility steepest descent */
280 struct vector_data Jgamma; /* Product of J and gamma */
281 struct vector_data newton; /* Dependent variables */
282 struct vector_data perturbed_newton; /* Perturbed Newton direction */
283 struct vector_data varstep1; /* 1st order in variables */
284 struct vector_data varstep2; /* 2nd order in variables */
285 struct vector_data varnewstep;/* hypothetical newton step in variables */
286 struct vector_data varstep; /* Step in variables */
287
288 real64 objective; /* Objective function evaluation */
289 real64 phi; /* Unconstrained minimizer */
290 real64 maxstep; /* Maximum step size allowed */
291 real64 progress; /* Steepest directional derivative */
292 real64 sigma; /* penalty parameter */
293 real64 sigma_const; /* penalty parameter */
294 real64 mu; /* Complementary gap */
295 real64 muaff; /* Complementary gap after newton */
296 real64 sigmamu; /* Complementary gap times penalty */
297 real64 comp; /* no. of complementary eqns. */
298
299
300 };
301
302
303 /*
304 * Integrity checks
305 * ----------------
306 * check_system(sys)
307 */
308
309 #define OK ((int)813029392)
310 #define DESTROYED ((int)103289182)
311 static int check_system(slv5_system_t sys)
312 /*
313 * Checks sys for NULL and for integrity.
314 */
315 {
316 if( sys == NULL ) {
317 FPRINTF(stderr,"ERROR: (slv5) check_system\n");
318 FPRINTF(stderr," NULL system handle.\n");
319 return 1;
320 }
321
322 switch( sys->integrity ) {
323 case OK:
324 return 0;
325 case DESTROYED:
326 FPRINTF(stderr,"ERROR: (slv5) check_system\n");
327 FPRINTF(stderr," System was recently destroyed.\n");
328 return 1;
329 default:
330 FPRINTF(stderr,"ERROR: (slv5) check_system\n");
331 FPRINTF(stderr," System reused or never allocated.\n");
332 return 1;
333 }
334 }
335
336 /*
337 * General input/output routines
338 * -----------------------------
339 * print_var_name(out,sys,var)
340 * print_rel_name(out,sys,rel)
341 */
342
343 #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
344 #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
345
346 /*
347 * Debug output routines
348 * ---------------------
349 * debug_delimiter(fp)
350 * debug_out_vector(fp,vec)
351 * debug_out_var_values(fp,sys)
352 * debug_out_rel_residuals(fp,sys)
353 * debug_out_jacobian(fp,sys)
354 * debug_write_array(fp,real64 *,length)
355 * debug_out_parameters(fp)
356 */
357
358 /*
359 * Outputs a hyphenated line.
360 */
361 static void debug_delimiter( FILE *fp)
362 {
363 int i;
364 for( i=0; i<60; i++ ) PUTC('-',fp);
365 PUTC('\n',fp);
366 }
367
368 #if DEBUG
369 /*
370 * Outputs a vector.
371 */
372 static void debug_out_vector( FILE *fp, slv5_system_t sys,
373 struct vector_data *vec)
374 {
375 int32 ndx;
376 FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n",
377 calc_sqrt_D0(vec->norm2), vec->accurate?"TRUE":"FALSE",
378 vec->rng->low,vec->rng->high);
379 FPRINTF(fp,"Vector --> ");
380 for( ndx=vec->rng->low ; ndx<=vec->rng->high ; ++ndx )
381 FPRINTF(fp, "%g ", vec->vec[ndx]);
382 PUTC('\n',fp);
383 }
384
385 /*
386 * Outputs all variable values in current block.
387 */
388 static void debug_out_var_values( FILE *fp, slv5_system_t sys)
389 {
390 int32 col;
391 struct var_variable *var;
392
393 FPRINTF(fp,"Var values --> \n");
394 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high ; col++ ) {
395 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
396 print_var_name(fp,sys,var);
397 FPRINTF(fp, "\nI Lb Value Ub Scale Col INom\n");
398 FPRINTF(fp,"%d\t%.4g\t%.4g\t%.4g\t%.4g\t%d\t%.4g\n",
399 var_sindex(var),var_lower_bound(var),var_value(var),
400 var_upper_bound(var),var_nominal(var),
401 col,sys->nominals.vec[col]);
402 }
403 }
404
405 /*
406 * Outputs all relation residuals in current block.
407 */
408 static void debug_out_rel_residuals( FILE *fp, slv5_system_t sys)
409 {
410 int32 row;
411
412 FPRINTF(fp,"Rel residuals --> \n");
413 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high ; row++ ) {
414 struct rel_relation *rel;
415 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
416 FPRINTF(fp," %g : ",rel_residual(rel));
417 print_rel_name(fp,sys,rel);
418 PUTC('\n',fp);
419 }
420 PUTC('\n',fp);
421 }
422
423 /*
424 * Outputs permutation and values of the nonzero elements in the
425 * the jacobian matrix.
426 */
427 static void debug_out_jacobian( FILE *fp, slv5_system_t sys)
428 {
429 mtx_coord_t nz;
430 real64 value;
431
432 nz.row = sys->J.reg.row.low;
433 for( ; nz.row <= sys->J.reg.row.high; ++(nz.row) ) {
434 FPRINTF(fp," Row %d (rel %d)\n", nz.row,
435 mtx_row_to_org(sys->J.mtx,nz.row));
436 nz.col = mtx_FIRST;
437 while( value = mtx_next_in_row(sys->J.mtx,&nz,&(sys->J.reg.col)),
438 nz.col != mtx_LAST ) {
439 FPRINTF(fp," Col %d (var %d) has value %g\n", nz.col,
440 mtx_col_to_org(sys->J.mtx,nz.col), value);
441 }
442 }
443 }
444
445 #endif
446
447 static void debug_write_array(FILE *fp,real64 *vec, int32 length)
448 {
449 int32 i;
450 for (i=0; i< length;i++)
451 FPRINTF(fp,"%.20g\n",vec[i]);
452 }
453
454 static char savlinfilename[]="SlvLinsol.dat.\0";
455 static char savlinfilebase[]="SlvLinsol.dat.\0";
456 static int savlinnum=0;
457
458 /* The number to postfix to savlinfilebase. increases with file accesses. */
459
460 /*
461 * Array/vector operations
462 * ----------------------------
463 * destroy_array(p)
464 * create_array(len,type)
465 * zero_vector(vec)
466 * copy_vector(vec1,vec2)
467 * prod = inner_product(vec1,vec2)
468 * norm2 = square_norm(vec)
469 * matrix_product(mtx,vec,prod,scale,transpose)
470 */
471
472 #define destroy_array(p) \
473 if( (p) != NULL ) ascfree((p))
474 #define create_array(len,type) \
475 ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
476 #define create_zero_array(len,type) \
477 ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
478
479 #define zero_vector(v) slv_zero_vector(v)
480 #define copy_vector(v,t) slv_copy_vector((v),(t))
481 #define inner_product(v,u) slv_inner_product((v),(u))
482 #define square_norm(v) slv_square_norm(v)
483 #define matrix_product(m,v,p,s,t) slv_matrix_product((m),(v),(p),(s),(t))
484
485 /*
486 * Calculation routines
487 * --------------------
488 * ok = calc_objective(sys)
489 * ok = calc_boundaries(sys)
490 * ok = calc_residuals(sys)
491 * ok = calc_mu(sys)
492 * ok = calc_newton_residuals(sys)
493 * ok = calc_muaff(sys)
494 * ok = calc_sigma(sys)
495 * ok = calc_sigmamu(sys)
496 * ok = calc_perturbed_residuals(sys)
497 * ok = calc_J(sys)
498 * calc_nominals(sys)
499 * calc_weights(sys)
500 * scale_J(sys)
501 * scale_variables(sys)
502 * scale_residuals(sys)
503 * scale_perturbed_residuals(sys)
504 * calc_pivots(sys)
505 * calc_rhs(sys)
506 * calc_newton(sys)
507 * calc_perturbed_newton(sys)
508 * calc_gamma(sys)
509 * calc_Jgamma(sys)
510 * calc_varstep1(sys)
511 * calc_varstep2(sys)
512 * calc_varnewstep(sys)
513 * calc_varstep(sys)
514 * calc_phi(sys)
515 */
516
517 /*
518 * Evaluate the objective function.
519 */
520 static boolean calc_objective( slv5_system_t sys)
521 {
522 calc_ok = TRUE;
523 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
524 sys->objective = (sys->obj ? relman_eval(sys->obj,&calc_ok,SAFE_CALC) : 0.0);
525 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
526 return calc_ok;
527 }
528
529 /*
530 * Evaluate all objectives.
531 */
532 static boolean calc_objectives( slv5_system_t sys)
533 {
534 int32 len,i;
535 static rel_filter_t rfilter;
536 struct rel_relation **rlist=NULL;
537 rfilter.matchbits = (REL_INCLUDED);
538 rfilter.matchvalue =(REL_INCLUDED);
539 rlist = slv_get_solvers_obj_list(SERVER);
540 len = slv_get_num_solvers_objs(SERVER);
541 calc_ok = TRUE;
542 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
543 for (i = 0; i < len; i++) {
544 if (rel_apply_filter(rlist[i],&rfilter)) {
545 relman_eval(rlist[i],&calc_ok,SAFE_CALC);
546 #if DEBUG
547 if (calc_ok == FALSE) {
548 FPRINTF(stderr,"error in calc_objectives\n");
549 calc_ok = TRUE;
550 }
551 #endif
552 }
553 }
554 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
555 return calc_ok;
556 }
557
558
559 /*
560 * Calculates all of the residuals of included inequalities.
561 * Returns true iff (calculations preceded without error and
562 * all inequalities are satisfied.)
563 */
564 static boolean calc_inequalities( slv5_system_t sys)
565 {
566 struct rel_relation **rp;
567 boolean satisfied=TRUE;
568 static rel_filter_t rfilter;
569 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
570 rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
571
572 calc_ok = TRUE;
573 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
574 for (rp=sys->rlist;*rp != NULL; rp++) {
575 if (rel_apply_filter(*rp,&rfilter)) {
576 relman_eval(*rp,&calc_ok,SAFE_CALC);
577 satisfied= satisfied &&
578 relman_calc_satisfied(*rp,FEAS_TOL);
579 }
580 }
581 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
582 return (calc_ok && satisfied);
583 }
584
585 /*
586 * Calculates the penalty of objective function
587 */
588 static void calc_comp( slv5_system_t sys )
589 {
590 int32 row;
591 struct rel_relation *rel;
592 real64 comp;
593
594 comp = 0.0;
595 row = sys->residuals.rng->low;
596 for( ; row <= sys->residuals.rng->high; row++ ) {
597 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
598 #if DEBUG
599 if (!rel) {
600 int r;
601 r=mtx_row_to_org(sys->J.mtx,row);
602 FPRINTF(ASCERR,"NULL relation found !!\n");
603 FPRINTF(ASCERR,"at row %d rel %d in calc_comp\n",(int)row,r);
604 FFLUSH(ASCERR);
605 }
606 #endif
607 if (rel_active(rel) && rel_included (rel) && rel_complementary(rel)) {
608 #if DEBUG
609 FPRINTF(ASCERR,"Complementary equation in slv5 \n");
610 #endif /* DEBUG */
611 comp = comp + 1.0;
612 }
613 }
614
615 #if DEBUG_ITERATION
616 FPRINTF(ASCERR," No. of complementary eqns = %g \n",comp);
617 #endif /* DEBUG_ITERATION */
618 sys->comp = comp;
619 }
620
621
622 /*
623 * Calculates all of the residuals in the current block and computes
624 * the residual norm for block status. Returns true iff calculations
625 * preceded without error.
626 */
627 static boolean calc_residuals( slv5_system_t sys)
628 {
629 int32 row;
630 struct rel_relation *rel;
631 double time0;
632
633 if( sys->residuals.accurate ) return TRUE;
634
635 calc_ok = TRUE;
636 row = sys->residuals.rng->low;
637 time0=tm_cpu_time();
638 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
639 for( ; row <= sys->residuals.rng->high; row++ ) {
640 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
641 #if DEBUG
642 if (!rel) {
643 int r;
644 r=mtx_row_to_org(sys->J.mtx,row);
645 FPRINTF(stderr,"NULL relation found !!\n");
646 FPRINTF(stderr,"at row %d rel %d in calc_residuals\n",(int)row,r);
647 FFLUSH(stderr);
648 }
649 #endif
650 sys->residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
651
652 if (strcmp(CONVOPT,"ABSOLUTE") == 0) {
653 relman_calc_satisfied(rel,FEAS_TOL);
654 } else if (strcmp(CONVOPT,"RELNOM_SCALE") == 0) {
655 relman_calc_satisfied_scaled(rel,FEAS_TOL);
656 }
657 }
658 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
659 sys->s.block.functime += (tm_cpu_time() -time0);
660 sys->s.block.funcs++;
661 square_norm( &(sys->residuals) );
662 sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
663 return(calc_ok);
664 }
665
666
667
668 /*
669 * Calculates the complementary gap in the current point
670 */
671 static boolean calc_mu( slv5_system_t sys)
672 {
673
674 int32 row;
675 struct rel_relation *rel;
676 real64 muk;
677
678 muk = 0.0;
679 row = sys->residuals.rng->low;
680
681 #if DEBUG
682 FPRINTF(ASCERR,"row low is = %d \n",row);
683 FPRINTF(ASCERR,"row high is = %d \n",sys->residuals.rng->high);
684 #endif /* DEBUG */
685
686 for(row ; row <= sys->residuals.rng->high; row++ ) {
687 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
688 #if DEBUG
689 if (!rel) {
690 int r;
691 r = mtx_row_to_org(sys->J.mtx,row);
692 FPRINTF(ASCERR,"NULL relation found !!\n");
693 FPRINTF(ASCERR,"at row %d rel %d in calc_mu \n",(int)row,r);
694 FFLUSH(ASCERR);
695 }
696 #endif
697 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
698 muk = muk + rel_residual(rel);
699 #if DEBUG
700 FPRINTF(ASCERR,"Complementary equation in calc_mu \n");
701 FPRINTF(ASCERR,"row = %d\n",row);
702 FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
703 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
704 FPRINTF(ASCERR,"Partial muk is = %g \n",muk);
705 #endif /* DEBUG */
706 }
707 }
708
709 if (sys->comp > 0.0) {
710 muk = muk / sys->comp;
711 } else {
712 muk = 0.0;
713 }
714
715 sys->mu = muk;
716
717 #if DEBUG_CENTERING
718 FPRINTF(ASCERR,"muk is = %g \n",sys->mu);
719 #endif /* DEBUG_CENTERING */
720
721 return(calc_ok);
722 }
723
724
725 /*
726 * Calculates all of the residuals in the current block and computes
727 * the residual norm for block status after a hypothetical Newton
728 * step. Returns true iff calculations
729 * preceded without error.
730 */
731 static boolean calc_newton_residuals( slv5_system_t sys)
732 {
733 int32 row;
734 struct rel_relation *rel;
735
736 if( sys->newton_residuals.accurate ) return TRUE;
737
738 calc_ok = TRUE;
739 row = sys->newton_residuals.rng->low;
740 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
741 for( ; row <= sys->newton_residuals.rng->high; row++ ) {
742 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
743 #if DEBUG
744 if (!rel) {
745 int r;
746 r=mtx_row_to_org(sys->J.mtx,row);
747 FPRINTF(stderr,"NULL relation found !!\n");
748 FPRINTF(stderr,"at row %d rel %d in calc_newton_residuals\n",(int)row,r);
749 FFLUSH(stderr);
750 }
751 #endif
752 sys->newton_residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
753 }
754 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
755 square_norm( &(sys->newton_residuals) );
756 return(calc_ok);
757 }
758
759 /*
760 * Calculates the complementary gap after a hypothetical Newton Step
761 */
762 static boolean calc_muaff( slv5_system_t sys)
763 {
764 int32 row;
765 struct rel_relation *rel;
766 real64 muaff;
767
768 muaff = 0.0;
769 row = sys->newton_residuals.rng->low;
770 for(row ; row <= sys->newton_residuals.rng->high; row++ ) {
771 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
772 #if DEBUG
773 if (!rel) {
774 int r;
775 r=mtx_row_to_org(sys->J.mtx,row);
776 FPRINTF(ASCERR,"NULL relation found !!\n");
777 FPRINTF(ASCERR,"at row %d rel %d in calc_muaff\n",(int)row,r);
778 FFLUSH(ASCERR);
779 }
780 #endif
781 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
782 muaff = muaff + sys->newton_residuals.vec[row];
783 #if DEBUG
784 FPRINTF(ASCERR,"Complementary equation in calc_muaff \n");
785 FPRINTF(ASCERR,"row = %d\n",row);
786 FPRINTF(ASCERR,"residual vector = %g\n",sys->newton_residuals.vec[row]);
787 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
788 FPRINTF(ASCERR,"Partial muaff is = %g \n",muaff);
789 #endif /* DEBUG */
790 }
791 }
792
793 if (sys->comp > 0.0) {
794 muaff = muaff / sys->comp;
795 } else {
796 muaff = 0.0;
797 }
798
799 sys->muaff = muaff;
800
801 #if DEBUG_CENTERING
802 FPRINTF(ASCERR,"muaff is = %g \n",sys->muaff);
803 #endif /* DEBUG_CENTERING */
804
805 return(calc_ok);
806 }
807
808
809 /*
810 * Calculates the penalty parameter sigma
811 */
812 static boolean calc_sigma( slv5_system_t sys)
813 {
814 real64 sigma, frac;
815
816 if ((sys->mu) > 0.0) {
817 frac = (sys->muaff) / (sys->mu);
818 sigma = (frac) * (frac) * (frac);
819 } else {
820 frac = 0.0;
821 sigma = 0.0;
822 }
823
824 sys->sigma = sigma;
825 #if DEBUG_CENTERING
826 FPRINTF(ASCERR,"sigma is = %g \n",sys->sigma);
827 #endif /* DEBUG_CENTERING */
828 return(calc_ok);
829 }
830
831 /*
832 * Calculates the penalty parameter time the complementary gap
833 */
834 static boolean calc_sigmamu( slv5_system_t sys)
835 {
836 real64 sigmamu;
837
838 if ((sys->mu) > 0.0 && (sys->sigma) > 0.0) {
839 sigmamu = (sys->mu) * (sys->sigma);
840 } else {
841 sigmamu = 0.0;
842 }
843
844 sys->sigmamu = sigmamu;
845 #if DEBUG_CENTERING
846 FPRINTF(ASCERR,"sigmamu is = %g \n",sys->sigmamu);
847 #endif /* DEBUG_CENTERING */
848 return(calc_ok);
849 }
850
851
852 /*
853 * Calculates the perturbed residuals in the current block and computes
854 * the perturbed residual norm. Returns true iff calculations
855 * preceded without error.
856 */
857 static boolean calc_perturbed_residuals( slv5_system_t sys)
858 {
859 int32 row;
860 struct rel_relation *rel;
861
862 if( sys->perturbed_residuals.accurate ) return TRUE;
863
864 calc_ok = TRUE;
865 row = sys->perturbed_residuals.rng->low;
866 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
867 for( ; row <= sys->perturbed_residuals.rng->high; row++ ) {
868 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
869 #if DEBUG
870 if (!rel) {
871 int r;
872 r=mtx_row_to_org(sys->J.mtx,row);
873 FPRINTF(stderr,"NULL relation found !!\n");
874 FPRINTF(stderr,"at row %d rel %d in calc_perturbed _residuals\n",
875 (int)row,r);
876 FFLUSH(stderr);
877 }
878 #endif
879 sys->perturbed_residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
880 if (rel_complementary(rel) && rel_active(rel)) {
881 sys->perturbed_residuals.vec[row] = sys->perturbed_residuals.vec[row] -
882 sys->sigmamu ;
883 }
884 }
885 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
886 square_norm( &(sys->perturbed_residuals) );
887 return(calc_ok);
888 }
889
890 /*
891 * Calculates the current block of the jacobian.
892 * It is initially unscaled.
893 */
894 static boolean calc_J( slv5_system_t sys)
895 {
896 int32 row;
897 var_filter_t vfilter;
898 double time0;
899 real64 resid;
900
901 if( sys->J.accurate )
902 return TRUE;
903
904 calc_ok = TRUE;
905 vfilter.matchbits = (VAR_INBLOCK | VAR_ACTIVE);
906 vfilter.matchvalue = (VAR_INBLOCK | VAR_ACTIVE);
907 time0=tm_cpu_time();
908 mtx_clear_region(sys->J.mtx,&(sys->J.reg));
909 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
910 struct rel_relation *rel;
911 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
912 relman_diffs(rel,&vfilter,sys->J.mtx,&resid,SAFE_CALC);
913 }
914 sys->s.block.jactime += (tm_cpu_time() - time0);
915 sys->s.block.jacs++;
916
917 if( --(sys->update.nominals) <= 0 ) sys->nominals.accurate = FALSE;
918 if( --(sys->update.weights) <= 0 ) sys->weights.accurate = FALSE;
919
920 linsolqr_matrix_was_changed(sys->J.sys);
921 return(calc_ok);
922 }
923
924 /*
925 * Retrieves the nominal values of all of the block variables,
926 * insuring that they are all strictly positive.
927 */
928 static void calc_nominals( slv5_system_t sys)
929 {
930 int32 col;
931 FILE *fp = MIF(sys);
932 if( sys->nominals.accurate ) return;
933 fp = MIF(sys);
934 col = sys->nominals.rng->low;
935 if(strcmp(SCALEOPT,"NONE") == 0 ||
936 strcmp(SCALEOPT,"ITERATIVE") == 0){
937 for( ; col <= sys->nominals.rng->high; col++ ) {
938 sys->nominals.vec[col] = 1;
939 }
940 } else {
941 for( ; col <= sys->nominals.rng->high; col++ ) {
942 struct var_variable *var;
943 real64 n;
944 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
945 n = var_nominal(var);
946 if( n <= 0.0 ) {
947 if( n == 0.0 ) {
948 n = TOO_SMALL;
949 FPRINTF(fp,"ERROR: (slv5) calc_nominals\n");
950 FPRINTF(fp," Variable ");
951 print_var_name(fp,sys,var);
952 FPRINTF(fp," \nhas nominal value of zero.\n");
953 FPRINTF(fp," Resetting to %g.\n",n);
954 var_set_nominal(var,n);
955 } else {
956 n = -n;
957 FPRINTF(fp,"ERROR: (slv5) calc_nominals\n");
958 FPRINTF(fp," Variable ");
959 print_var_name(fp,sys,var);
960 FPRINTF(fp," \nhas negative nominal value.\n");
961 FPRINTF(fp," Resetting to %g.\n",n);
962 var_set_nominal(var,n);
963 }
964 }
965 #if DEBUG
966 FPRINTF(fp,"Column %d is");
967 print_var_name(fp,sys,var);
968 FPRINTF(fp,"\nScaling of column %d is %g\n",col,n);
969 #endif
970 sys->nominals.vec[col] = n;
971 }
972 }
973 square_norm( &(sys->nominals) );
974 sys->update.nominals = UPDATE_NOMINALS;
975 sys->nominals.accurate = TRUE;
976 }
977
978 /*
979 * Calculates the weights of all of the block relations
980 * to scale the rows of the Jacobian.
981 */
982 static void calc_weights( slv5_system_t sys)
983 {
984 mtx_coord_t nz;
985 real64 sum;
986
987 if( sys->weights.accurate )
988 return;
989
990 nz.row = sys->weights.rng->low;
991 if(strcmp(SCALEOPT,"NONE") == 0 ||
992 strcmp(SCALEOPT,"ITERATIVE") == 0) {
993 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
994 sys->weights.vec[nz.row] = 1;
995 }
996 } else if (strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
997 strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0) {
998 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
999 sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
1000 sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
1001 }
1002 } else if (strcmp(SCALEOPT,"RELNOM") == 0 ||
1003 strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0) {
1004 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
1005 sys->weights.vec[nz.row] =
1006 1.0/rel_nominal(sys->rlist[mtx_row_to_org(sys->J.mtx,nz.row)]);
1007 }
1008 }
1009 square_norm( &(sys->weights) );
1010 sys->update.weights = UPDATE_WEIGHTS;
1011 sys->residuals.accurate = FALSE;
1012 sys->weights.accurate = TRUE;
1013 }
1014
1015 /*
1016 * Scales the jacobian.
1017 */
1018 static void scale_J( slv5_system_t sys)
1019 {
1020 int32 row;
1021 int32 col;
1022
1023 if( sys->J.accurate ) return;
1024
1025 calc_nominals(sys);
1026 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ )
1027 mtx_mult_col(sys->J.mtx,col,sys->nominals.vec[col],&(sys->J.reg.row));
1028
1029 calc_weights(sys);
1030 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ )
1031 mtx_mult_row(sys->J.mtx,row,sys->weights.vec[row],&(sys->J.reg.col));
1032 }
1033
1034 static void jacobian_scaled(slv5_system_t sys)
1035 {
1036 int32 col;
1037 if (DUMPCNORM) {
1038 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
1039 real64 cnorm;
1040 cnorm =
1041 calc_sqrt_D0(mtx_sum_sqrs_in_col(sys->J.mtx,col,&(sys->J.reg.row)));
1042 if (cnorm >CNHIGH || cnorm <CNLOW) {
1043 FPRINTF(stderr,"[col %d org %d] %g\n", col,
1044 mtx_col_to_org(sys->J.mtx,col), cnorm);
1045 }
1046 }
1047 }
1048
1049 sys->update.jacobian = UPDATE_JACOBIAN;
1050 sys->J.accurate = TRUE;
1051 sys->J.singular = FALSE; /* yet to be determined */
1052 #if DEBUG
1053 FPRINTF(LIF(sys),"\nJacobian: \n");
1054 debug_out_jacobian(LIF(sys),sys);
1055 #endif
1056 }
1057
1058 static void scale_variables( slv5_system_t sys)
1059 {
1060 int32 col;
1061
1062 if( sys->variables.accurate ) return;
1063
1064 col = sys->variables.rng->low;
1065 for( ; col <= sys->variables.rng->high; col++ ) {
1066 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1067 sys->variables.vec[col] = var_value(var)/sys->nominals.vec[col];
1068 }
1069 square_norm( &(sys->variables) );
1070 sys->variables.accurate = TRUE;
1071 #if DEBUG
1072 FPRINTF(LIF(sys),"Variables: ");
1073 debug_out_vector(LIF(sys),sys,&(sys->variables));
1074 #endif
1075 }
1076
1077 /*
1078 * Scales the previously calculated residuals.
1079 */
1080 static void scale_residuals( slv5_system_t sys)
1081 {
1082 int32 row;
1083
1084 if( sys->residuals.accurate ) return;
1085
1086 row = sys->residuals.rng->low;
1087 for( ; row <= sys->residuals.rng->high; row++ ) {
1088 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1089 sys->residuals.vec[row] = rel_residual(rel)*sys->weights.vec[row];
1090 }
1091 square_norm( &(sys->residuals) );
1092 sys->residuals.accurate = TRUE;
1093 #if DEBUG
1094 FPRINTF(LIF(sys),"Residuals: ");
1095 debug_out_vector(LIF(sys),sys,&(sys->residuals));
1096 #endif
1097 }
1098
1099 /*
1100 * Scales the previously calculated residuals.
1101 */
1102 static void scale_perturbed_residuals( slv5_system_t sys)
1103 {
1104 int32 row;
1105
1106 if( sys->perturbed_residuals.accurate ) return;
1107
1108 row = sys->perturbed_residuals.rng->low;
1109 for( ; row <= sys->perturbed_residuals.rng->high; row++ ) {
1110 sys->perturbed_residuals.vec[row] = sys->perturbed_residuals.vec[row]
1111 * sys->weights.vec[row];
1112 }
1113 square_norm( &(sys->perturbed_residuals) );
1114 sys->perturbed_residuals.accurate = TRUE;
1115 #if DEBUG
1116 FPRINTF(LIF(sys),"Perturbed Residuals: ");
1117 debug_out_vector(LIF(sys),sys,&(sys->perturbed_residuals));
1118 #endif
1119 }
1120
1121 /*
1122 * Calculates relnoms for all relations in sys
1123 * using variable nominals.
1124 */
1125 static void calc_relnoms(slv5_system_t sys)
1126 {
1127 int32 row, col;
1128 struct var_variable *var;
1129 struct rel_relation *rel;
1130 real64 *var_list;
1131 var_list = create_array(sys->cap,real64);
1132 col = 0;
1133 var = sys->vlist[col];
1134 /* store current variable values and
1135 set variable value to nominal value */
1136 while(var != NULL){
1137 var_list[col] = var_value(var);
1138 var_set_value(var, var_nominal(var));
1139 col++;
1140 var = sys->vlist[col];
1141 }
1142 row = 0;
1143 rel = sys->rlist[row];
1144 /* calculate relation nominal */
1145 while(rel != NULL){
1146 relman_scale(rel);
1147 row++;
1148 rel = sys->rlist[row];
1149 }
1150 col = 0;
1151 var = sys->vlist[col];
1152 /* restore variable values */
1153 while(var != NULL){
1154 var_set_value(var, var_list[col]);
1155 col++;
1156 var = sys->vlist[col];
1157 }
1158 destroy_array(var_list);
1159 }
1160
1161 /*
1162 * Returns the maximum ratio of magnitudes of any two nonzero
1163 * elements in the same column of mtx. Only considers elements
1164 * in region reg.
1165 */
1166 static real64 col_max_ratio(mtx_matrix_t *mtx, mtx_region_t *reg)
1167 {
1168 real64 ratio;
1169 real64 max_ratio;
1170 real64 num, denom, dummy;
1171 mtx_coord_t coord;
1172 max_ratio = 0;
1173 for(coord.col = reg->col.low;
1174 coord.col <= reg->col.high; coord.col++) {
1175 ratio = 0;
1176 num = mtx_col_max(*mtx,&(coord),&(reg->row),&(dummy));
1177 denom = mtx_col_min(*mtx,&(coord),&(reg->row),&(dummy),1e-7);
1178 if(denom >0){
1179 ratio = num/denom;
1180 }
1181 if(ratio > max_ratio){
1182 max_ratio = ratio;
1183 }
1184 }
1185 if(max_ratio == 0){
1186 max_ratio = 1;
1187 }
1188 return max_ratio;
1189 }
1190
1191
1192 /*
1193 * Returns the maximum ratio of magnitudes of any two nonzero
1194 * elements in the same row of mtx. Only considers elements
1195 * in region reg.
1196 */
1197 static real64 row_max_ratio(mtx_matrix_t *mtx,
1198 mtx_region_t *reg)
1199 {
1200 real64 ratio;
1201 real64 max_ratio;
1202 real64 num, denom, dummy;
1203 mtx_coord_t coord;
1204 max_ratio = 0;
1205 for(coord.row = reg->row.low;
1206 coord.row <= reg->row.high; coord.row++) {
1207 ratio = 0;
1208 num = mtx_row_max(*mtx,&(coord),&(reg->col),&(dummy));
1209 denom = mtx_row_min(*mtx,&(coord),&(reg->col),&(dummy),1e-7);
1210 if(denom >0){
1211 ratio = num/denom;
1212 }
1213 if(ratio > 10000000){
1214 /* FPRINTF(stderr,"HELPME\n");*/
1215 }
1216 if(ratio > max_ratio){
1217 max_ratio = ratio;
1218 }
1219 }
1220 if(max_ratio == 0){
1221 max_ratio = 1;
1222 }
1223 return max_ratio;
1224 }
1225
1226 /*
1227 * Calculates scaling factor suggested by Fourer.
1228 * For option = 0, returns scaling factor for
1229 * row number loc.
1230 * For option = 1, returns scaling factor for
1231 * col number loc.
1232 */
1233 static real64 calc_fourer_scale(mtx_matrix_t mtx,
1234 mtx_region_t reg,
1235 int32 loc,
1236 int32 option)
1237 {
1238 mtx_coord_t coord;
1239 real64 min, max, dummy;
1240 real64 scale;
1241 if(option == 0){
1242 if((loc < reg.row.low) || (loc > reg.row.high)){
1243 return 1;
1244 }
1245 coord.row = loc;
1246 min = mtx_row_min(mtx,&(coord),&(reg.col),&(dummy),1e-7);
1247 max = mtx_row_max(mtx,&(coord),&(reg.col),&(dummy));
1248 scale = min*max;
1249 if(scale > 0){
1250 scale = sqrt(scale);
1251 } else {
1252 scale = 1;
1253 }
1254 return scale;
1255 } else {
1256 if(loc < reg.col.low || loc > reg.col.high){
1257 return 1;
1258 }
1259 coord.col = loc;
1260 min = mtx_col_min(mtx,&(coord),&(reg.row),&(dummy),1e-7);
1261 max = mtx_col_max(mtx,&(coord),&(reg.row),&(dummy));
1262 scale = min*max;
1263 if(scale > 0){
1264 scale = sqrt(scale);
1265 } else {
1266 scale = 1;
1267 }
1268 return scale;
1269 }
1270 }
1271
1272 /*
1273 * This funcion is an implementation of the scaling
1274 * routine by Fourer on p304 of Mathematical Programing
1275 * vol 23, (1982).
1276 * This function will scale the Jacobian and store the scaling
1277 * factors in sys->nominals and sys->weights.
1278 * If the Jacobian has been previously scaled
1279 * by another method (during this iteration) then these vectors
1280 * should contain the scale factors used in that scaling.
1281 */
1282 static void scale_J_iterative(slv5_system_t sys)
1283 {
1284 real64 rho_col_old, rho_col_new;
1285 real64 rho_row_old, rho_row_new;
1286 int32 k;
1287 int32 done;
1288 int32 row, col;
1289 real64 *colvec = sys->nominals.vec;
1290 real64 *rowvec = sys->weights.vec;
1291 real64 rowscale, colscale;
1292 rho_col_old = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1293 rho_row_old = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1294 k = 0;
1295 done = 0;
1296 while(done == 0){
1297 k++;
1298 for(row = sys->J.reg.row.low;
1299 row <= sys->J.reg.row.high; row++){
1300 rowscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,row,0);
1301 mtx_mult_row(sys->J.mtx,row,rowscale,&(sys->J.reg.col));
1302 rowvec[row] *= rowscale;
1303 }
1304 for(col = sys->J.reg.col.low;
1305 col <= sys->J.reg.col.high; col++){
1306 colscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,col,1);
1307 mtx_mult_col(sys->J.mtx,col,colscale,&(sys->J.reg.row));
1308 colvec[col] *= colscale;
1309 }
1310 rho_col_new = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1311 rho_row_new = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1312 if((rho_col_new >= ITSCALETOL*rho_col_old &&
1313 rho_row_new >= ITSCALETOL*rho_row_old)
1314 || k >= ITSCALELIM){
1315 done = 1;
1316 } else {
1317 rho_row_old = rho_row_new;
1318 rho_col_old = rho_col_new;
1319 }
1320 }
1321 square_norm( &(sys->nominals) );
1322 sys->update.nominals = UPDATE_NOMINALS;
1323 sys->nominals.accurate = TRUE;
1324
1325 square_norm( &(sys->weights) );
1326 sys->update.weights = UPDATE_WEIGHTS;
1327 sys->residuals.accurate = FALSE;
1328 sys->weights.accurate = TRUE;
1329 }
1330
1331 /*
1332 * Scale system dependent on interface parameters
1333 */
1334 static void scale_perturbed_system( slv5_system_t sys )
1335 {
1336 if(strcmp(SCALEOPT,"NONE") == 0){
1337 scale_perturbed_residuals(sys);
1338 return;
1339 }
1340 if(strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
1341 strcmp(SCALEOPT,"RELNOM") == 0){
1342 scale_perturbed_residuals(sys);
1343 return;
1344 }
1345 }
1346
1347 /*
1348 * Scale system dependent on interface parameters
1349 */
1350 static void scale_system( slv5_system_t sys )
1351 {
1352 if(strcmp(SCALEOPT,"NONE") == 0){
1353 if(sys->J.accurate == FALSE){
1354 calc_nominals(sys);
1355 calc_weights(sys);
1356 jacobian_scaled(sys);
1357 }
1358 scale_variables(sys);
1359 scale_residuals(sys);
1360 return;
1361 }
1362 if(strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
1363 strcmp(SCALEOPT,"RELNOM") == 0){
1364 if(sys->J.accurate == FALSE){
1365 scale_J(sys);
1366 jacobian_scaled(sys);
1367 }
1368 scale_variables(sys);
1369 scale_residuals(sys);
1370 return;
1371 }
1372 if(strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0 ||
1373 strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0){
1374 if(sys->J.accurate == FALSE){
1375 --sys->update.iterative;
1376 if(sys->update.iterative <= 0) {
1377 scale_J(sys);
1378 scale_J_iterative(sys);
1379 sys->update.iterative =
1380 UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1381 } else {
1382 sys->weights.accurate = TRUE;
1383 sys->nominals.accurate = TRUE;
1384 scale_J(sys); /* will use current scaling vectors */
1385 }
1386 jacobian_scaled(sys);
1387 }
1388 scale_variables(sys);
1389 scale_residuals(sys);
1390 return;
1391 }
1392 if(strcmp(SCALEOPT,"ITERATIVE") == 0){
1393 if(sys->J.accurate == FALSE){
1394 --sys->update.iterative;
1395 if(sys->update.iterative <= 0) {
1396 calc_nominals(sys);
1397 calc_weights(sys);
1398 scale_J_iterative(sys);
1399 sys->update.iterative =
1400 UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1401 } else {
1402 sys->weights.accurate = TRUE;
1403 sys->nominals.accurate = TRUE;
1404 scale_J(sys); /* will use current scaling vectors */
1405 }
1406 jacobian_scaled(sys);
1407 }
1408 scale_variables(sys);
1409 scale_residuals(sys);
1410 }
1411 return;
1412 }
1413
1414
1415 /*
1416 * Obtain the equations and variables which
1417 * are able to be pivoted.
1418 * return value is the row rank deficiency, which we hope is 0.
1419 */
1420 static int calc_pivots(slv5_system_t sys)
1421 {
1422 int row_rank_defect=0, oldtiming;
1423 linsolqr_system_t lsys = sys->J.sys;
1424 FILE *fp = LIF(sys);
1425
1426 oldtiming = g_linsolqr_timing;
1427 g_linsolqr_timing =LINTIME;
1428 linsolqr_factor(lsys,sys->J.fm); /* factor */
1429 g_linsolqr_timing = oldtiming;
1430
1431 sys->J.rank = linsolqr_rank(lsys);
1432 sys->J.singular = FALSE;
1433 row_rank_defect = sys->J.reg.row.high -
1434 sys->J.reg.row.low+1 - sys->J.rank;
1435 if( row_rank_defect > 0 ) {
1436 int32 row,krow;
1437 mtx_sparse_t *uprows=NULL;
1438 sys->J.singular = TRUE;
1439 uprows = linsolqr_unpivoted_rows(lsys);
1440 if (uprows !=NULL) {
1441 for( krow=0; krow < uprows->len ; krow++ ) {
1442 int32 org_row;
1443 struct rel_relation *rel;
1444
1445 org_row = uprows->idata[krow];
1446 row = mtx_org_to_row(sys->J.mtx,org_row);
1447 rel = sys->rlist[org_row];
1448 FPRINTF(fp,"%-40s ---> ","Relation not pivoted");
1449 print_rel_name(fp,sys,rel);
1450 PUTC('\n',fp);
1451
1452 /*
1453 * assign zeros to the corresponding weights
1454 * so that subsequent calls to "scale_residuals"
1455 * will only measure the pivoted equations.
1456 */
1457 sys->weights.vec[row] = 0.0;
1458 sys->residuals.vec[row] = 0.0;
1459 sys->residuals.accurate = FALSE;
1460 sys->perturbed_residuals.vec[row] = 0.0;
1461 sys->perturbed_residuals.accurate = FALSE;
1462 sys->newton_residuals.vec[row] = 0.0;
1463 sys->newton_residuals.accurate = FALSE;
1464 mtx_mult_row(sys->J.mtx,row,0.0,&(sys->J.reg.col));
1465 }
1466 mtx_destroy_sparse(uprows);
1467 }
1468 if( !sys->residuals.accurate ) {
1469 square_norm( &(sys->residuals) );
1470 sys->residuals.accurate = TRUE;
1471 sys->update.weights = 0; /* re-compute weights next iteration. */
1472 }
1473 }
1474 if( sys->J.rank < sys->J.reg.col.high-sys->J.reg.col.low+1 ) {
1475 int32 col,kcol;
1476 mtx_sparse_t *upcols=NULL;
1477 if (NOTNULL(upcols)) {
1478 for( kcol=0; upcols != NULL && kcol < upcols->len ; kcol++ ) {
1479 int32 org_col;
1480 struct var_variable *var;
1481
1482 org_col = upcols->idata[kcol];
1483 col = mtx_org_to_col(sys->J.mtx,org_col);
1484 var = sys->vlist[org_col];
1485 FPRINTF(fp,"%-40s ---> ","Variable not pivoted");
1486 print_var_name(fp,sys,var);
1487 PUTC('\n',fp);
1488 /*
1489 * If we're not optimizing (everything should be
1490 * pivotable) or this was one of the dependent variables,
1491 * consider this variable as if it were fixed.
1492 */
1493 if( col <= sys->J.reg.col.high ) {
1494 mtx_mult_col(sys->J.mtx,col,0.0,&(sys->J.reg.row));
1495 }
1496 }
1497 mtx_destroy_sparse(upcols);
1498 }
1499 }
1500 if (SHOW_LESS_IMPT) {
1501 FPRINTF(LIF(sys),"%-40s ---> %d (%s)\n","Jacobian rank", sys->J.rank,
1502 sys->J.singular ? "deficient":"full");
1503 FPRINTF(LIF(sys),"%-40s ---> %g\n","Smallest pivot",
1504 linsolqr_smallest_pivot(sys->J.sys));
1505 }
1506 return row_rank_defect;
1507 }
1508
1509
1510 /*
1511 * Calculates just the jacobian RHS. This function should be used to
1512 * supplement calculation of the jacobian. The vector vec must
1513 * already be calculated and scaled so as to simply be added to the
1514 * rhs. Caller is responsible for initially zeroing the rhs vector.
1515 */
1516 static void calc_rhs(slv5_system_t sys, struct vector_data *vec,
1517 real64 scalar, boolean transpose)
1518 {
1519 if( transpose ) { /* vec is indexed by col */
1520 int32 col;
1521 for( col=vec->rng->low; col<=vec->rng->high; col++ ) {
1522 sys->J.rhs[mtx_col_to_org(sys->J.mtx,col)] += scalar*vec->vec[col];
1523 }
1524 } else { /* vec is indexed by row */
1525 int32 row;
1526 for( row=vec->rng->low; row<=vec->rng->high; row++ ) {
1527 sys->J.rhs[mtx_row_to_org(sys->J.mtx,row)] += scalar*vec->vec[row];
1528 }
1529 }
1530 linsolqr_rhs_was_changed(sys->J.sys,sys->J.rhs);
1531 }
1532
1533 /*
1534 * Calculate the gamma vector.
1535 */
1536 static void calc_gamma( slv5_system_t sys)
1537 {
1538 if( sys->gamma.accurate )
1539 return;
1540 if ( sys->mu > 0.0 ) {
1541 matrix_product(sys->J.mtx, &(sys->perturbed_residuals),
1542 &(sys->gamma), -1.0, TRUE);
1543 } else {
1544 matrix_product(sys->J.mtx, &(sys->residuals),
1545 &(sys->gamma), -1.0, TRUE);
1546 }
1547 square_norm( &(sys->gamma) );
1548 sys->gamma.accurate = TRUE;
1549 #if DEBUG
1550 FPRINTF(LIF(sys),"Gamma: ");
1551 debug_out_vector(LIF(sys),sys,&(sys->gamma));
1552 #endif
1553 }
1554
1555 /*
1556 * Calculate the Jgamma vector.
1557 */
1558 static void calc_Jgamma( slv5_system_t sys)
1559 {
1560 if( sys->Jgamma.accurate )
1561 return;
1562
1563 matrix_product(sys->J.mtx, &(sys->gamma),
1564 &(sys->Jgamma), 1.0, FALSE);
1565 square_norm( &(sys->Jgamma) );
1566 sys->Jgamma.accurate = TRUE;
1567 #if DEBUG
1568 FPRINTF(LIF(sys),"Jgamma: ");
1569 debug_out_vector(LIF(sys),sys,&(sys->Jgamma));
1570 #endif
1571 }
1572
1573 /*
1574 * Computes a step to solve the linearized equations.
1575 */
1576 static void calc_newton( slv5_system_t sys)
1577 {
1578 linsolqr_system_t lsys = sys->J.sys;
1579 int32 col;
1580
1581 if( sys->newton.accurate )
1582 return;
1583
1584 sys->J.rhs = linsolqr_get_rhs(lsys,0);
1585 mtx_zero_real64(sys->J.rhs,sys->cap);
1586 calc_rhs(sys, &(sys->residuals), -1.0, FALSE);
1587 linsolqr_solve(lsys,sys->J.rhs);
1588 col = sys->newton.rng->low;
1589 for( ; col <= sys->newton.rng->high; col++ ) {
1590 sys->newton.vec[col] =
1591 linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1592 }
1593 if (SAVLIN) {
1594 FILE *ldat;
1595 int32 ov;
1596 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1597 ldat=fopen(savlinfilename,"w");
1598 FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1599 sys->s.iteration);
1600 debug_write_array(ldat,sys->J.rhs,sys->cap);
1601 FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1602 for(ov=0 ; ov < sys->cap; ov++ )
1603 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1604 fclose(ldat);
1605 }
1606 square_norm( &(sys->newton) );
1607 sys->newton.accurate = TRUE;
1608 #if DEBUG
1609 FPRINTF(LIF(sys),"Newton: ");
1610 debug_out_vector(LIF(sys),sys,&(sys->newton));
1611 #endif
1612 }
1613
1614 /*
1615 * Computes a step to solve the linearized equations.
1616 */
1617 static void calc_perturbed_newton( slv5_system_t sys)
1618 {
1619 linsolqr_system_t lsys = sys->J.sys;
1620 int32 col;
1621
1622 if( sys->perturbed_newton.accurate )
1623 return;
1624
1625 sys->J.rhs = linsolqr_get_rhs(lsys,1);
1626 mtx_zero_real64(sys->J.rhs,sys->cap);
1627 calc_rhs(sys, &(sys->perturbed_residuals), -1.0, FALSE);
1628 linsolqr_solve(lsys,sys->J.rhs);
1629 col = sys->perturbed_newton.rng->low;
1630 for( ; col <= sys->perturbed_newton.rng->high; col++ ) {
1631 sys->perturbed_newton.vec[col] =
1632 linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1633 }
1634 if (SAVLIN) {
1635 FILE *ldat;
1636 int32 ov;
1637 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1638 ldat=fopen(savlinfilename,"w");
1639 FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1640 sys->s.iteration);
1641 debug_write_array(ldat,sys->J.rhs,sys->cap);
1642 FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1643 for(ov=0 ; ov < sys->cap; ov++ )
1644 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1645 fclose(ldat);
1646 }
1647 square_norm( &(sys->perturbed_newton) );
1648 sys->perturbed_newton.accurate = TRUE;
1649 #if DEBUG
1650 FPRINTF(LIF(sys),"Perturbed Newton: ");
1651 debug_out_vector(LIF(sys),sys,&(sys->perturbed_newton));
1652 #endif
1653 }
1654
1655
1656 /*
1657 * Calculate the 1st order descent direction for phi
1658 * in the variables.
1659 */
1660 static void calc_varstep1( slv5_system_t sys)
1661 {
1662 if( sys->varstep1.accurate ) {
1663 return;
1664 }
1665
1666 copy_vector(&(sys->gamma),&(sys->varstep1));
1667 sys->varstep1.norm2 = sys->gamma.norm2;
1668 sys->varstep1.accurate = TRUE;
1669 #if DEBUG
1670 FPRINTF(LIF(sys),"Varstep1: ");
1671 debug_out_vector(LIF(sys),sys,&(sys->varstep1));
1672 #endif
1673 }
1674
1675
1676
1677 /*
1678 * Calculate the 2nd order descent direction for phi
1679 * in the variables.
1680 */
1681 static void calc_varstep2( slv5_system_t sys)
1682 {
1683 if( sys->varstep2.accurate )
1684 return;
1685 if (sys->mu > 0.0) {
1686 copy_vector(&(sys->perturbed_newton),&(sys->varstep2));
1687 sys->varstep2.norm2 = sys->perturbed_newton.norm2;
1688 } else {
1689 copy_vector(&(sys->newton),&(sys->varstep2));
1690 sys->varstep2.norm2 = sys->newton.norm2;
1691 }
1692 sys->varstep2.accurate = TRUE;
1693 #if DEBUG
1694 FPRINTF(LIF(sys),"Varstep2: ");
1695 debug_out_vector(LIF(sys),sys,&(sys->varstep2));
1696 #endif
1697 }
1698
1699 /*
1700 * Calculate the hypothetical netown direction
1701 * in the variables.
1702 */
1703 static void calc_varnewstep( slv5_system_t sys)
1704 {
1705 if( sys->varnewstep.accurate )
1706 return;
1707 copy_vector(&(sys->newton),&(sys->varnewstep));
1708 sys->varnewstep.norm2 = sys->newton.norm2;
1709
1710 sys->varnewstep.accurate = TRUE;
1711 #if DEBUG
1712 FPRINTF(LIF(sys),"Varnewstep: ");
1713 debug_out_vector(LIF(sys),sys,&(sys->varnewstep));
1714 #endif
1715 }
1716
1717 /*
1718 * Computes the global minimizing function Phi.
1719 */
1720 static void calc_phi( slv5_system_t sys)
1721 {
1722 sys->phi = 0.5*sys->residuals.norm2;
1723
1724 if (sys->mu > 0.0) {
1725 sys->phi = 0.5 *sys->perturbed_residuals.norm2;
1726 }
1727
1728 }
1729
1730 /*
1731 * OK. Here's where we compute the actual step to be taken. It will
1732 * be some linear combination of the 1st order and 2nd order steps.
1733 */
1734
1735 typedef real64 sym_2x2_t[3]; /* Stores symmetric 2x2 matrices */
1736
1737 struct parms_t {
1738 real64 low,high,guess; /* Used to search for parameter */
1739 };
1740
1741 struct calc_step_vars {
1742 sym_2x2_t coef1, coef2;
1743 real64 rhs[2]; /* RHS for 2x2 system */
1744 struct parms_t parms;
1745 real64 alpha1, alpha2;
1746 real64 error; /* Error between step norm and sys->maxstep */
1747 };
1748
1749
1750 /*
1751 * Calculates 2x2 system (coef1,coef2,rhs).
1752 */
1753 static void calc_2x2_system(slv5_system_t sys, struct calc_step_vars *vars)
1754 {
1755 vars->coef1[0] = (2.0*sys->phi/sys->newton.norm2)*
1756 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1757 vars->coef1[1] = 1.0;
1758 vars->coef1[2] = (sys->Jgamma.norm2/sys->gamma.norm2)*
1759 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1760
1761 vars->coef2[0] = 1.0;
1762 vars->coef2[1] = 2.0*sys->phi/
1763 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1764 vars->coef2[2] = 1.0;
1765
1766 vars->rhs[0] = 2.0*sys->phi/
1767 sys->maxstep/calc_sqrt_D0(sys->gamma.norm2);
1768 vars->rhs[1] = calc_sqrt_D0(sys->newton.norm2)/sys->maxstep;
1769 }
1770
1771 /*
1772 * Determines alpha1 and alpha2 from the parameter (guess).
1773 */
1774 static void coefs_from_parm( slv5_system_t sys, struct calc_step_vars *vars)
1775 {
1776
1777 sym_2x2_t coef; /* Actual coefficient matrix */
1778 real64 det; /* Determinant of coefficient matrix */
1779 int i;
1780
1781 for( i=0 ; i<3 ; ++i ) coef[i] =
1782 vars->coef1[i] + vars->parms.guess * vars->coef2[i];
1783 det = coef[0]*coef[2] - coef[1]*coef[1];
1784 if( det < 0.0 )
1785 FPRINTF(MIF(sys),"%-40s ---> %g\n",
1786 " Unexpected negative determinant!",det);
1787 if( det <= DETZERO ) {
1788 /*
1789 * varstep2 and varstep1 are essentially parallel:
1790 * adjust length of either
1791 */
1792 vars->alpha2 = 0.0;
1793 vars->alpha1 = 1.0;
1794 } else {
1795 vars->alpha2 = (vars->rhs[0]*coef[2] - vars->rhs[1]*coef[1])/det;
1796 vars->alpha1 = (vars->rhs[1]*coef[0] - vars->rhs[0]*coef[1])/det;
1797 }
1798 }
1799
1800
1801 /*
1802 * Computes step vector length based on 1st order and 2nd order
1803 * vectors and their coefficients.
1804 */
1805 static real64 step_norm2( slv5_system_t sys, struct calc_step_vars *vars)
1806 {
1807 return sys->maxstep*sys->maxstep*
1808 (vars->alpha2 * vars->alpha2 +
1809 vars->alpha2 * vars->alpha1 * sys->phi/
1810 calc_sqrt_D0(sys->varstep2.norm2 )/
1811 calc_sqrt_D0(sys->varstep1.norm2 ) +
1812 vars->alpha1 * vars->alpha1);
1813 }
1814
1815 /*
1816 * Re-guesses the parameters based on
1817 * step size vs. target value.
1818 */
1819 static void adjust_parms( slv5_system_t sys, struct calc_step_vars *vars)
1820 {
1821 vars->error = (calc_sqrt_D0(step_norm2(sys,vars))/sys->maxstep) - 1.0;
1822 if( vars->error > 0.0 ) {
1823 /* Increase parameter (to decrease step length) */
1824 vars->parms.low = vars->parms.guess;
1825 vars->parms.guess = (vars->parms.high>3.0*vars->parms.guess)
1826 ? 2.0*vars->parms.guess
1827 : 0.5*(vars->parms.low + vars->parms.high);
1828 } else {
1829 /* Decrease parameter (to increase step norm) */
1830 vars->parms.high = vars->parms.guess;
1831 vars->parms.guess = 0.5*(vars->parms.low + vars->parms.high);
1832 }
1833 }
1834
1835 /*
1836 * Computes the step based on the coefficients in vars.
1837 */
1838 static void compute_step( slv5_system_t sys, struct calc_step_vars *vars)
1839 {
1840 int32 col;
1841 real64 tot1_norm2, tot2_norm2;
1842
1843 tot1_norm2 = sys->varstep1.norm2;
1844 tot2_norm2 = sys->varstep2.norm2;
1845 if( !sys->varstep.accurate ) {
1846 for( col=sys->varstep.rng->low ; col<=sys->varstep.rng->high ; ++col )
1847 if( (vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
1848 sys->varstep.vec[col] = sys->maxstep *
1849 sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2);
1850 } else if( (vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
1851 sys->varstep.vec[col] = sys->maxstep *
1852 sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2);
1853 } else if( (vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
1854 sys->varstep.vec[col] = sys->maxstep*
1855 (vars->alpha2*sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2) +
1856 vars->alpha1*sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2));
1857 }
1858 sys->varstep.accurate = TRUE;
1859 }
1860 #if DEBUG
1861 FPRINTF(LIF(sys),"Varstep: ");
1862 debug_out_vector(LIF(sys),sys,&(sys->varstep));
1863
1864 #endif
1865 }
1866
1867 /*
1868 * Calculates step vector, based on sys->maxstep, and the varstep2/
1869 * varstep1. Nothing is assumed to be
1870 * calculated, except the weights and the jacobian (scaled). Also,
1871 * the step is not checked for legitimacy.
1872 * NOTE: the step is scaled.
1873 */
1874 static void calc_step( slv5_system_t sys, int minor)
1875 {
1876
1877 struct calc_step_vars vars;
1878 real64 tot1_norm2, tot2_norm2;
1879
1880 if( sys->varstep.accurate )
1881 return;
1882 if (SHOW_LESS_IMPT) {
1883 FPRINTF(LIF(sys),"\n%-40s ---> %d\n", " Step trial",minor);
1884 }
1885
1886 tot1_norm2 = sys->varstep1.norm2;
1887 tot2_norm2 = sys->varstep2.norm2;
1888 if( (tot1_norm2 == 0.0) && (tot2_norm2 == 0.0) ) {
1889 /* Take no step at all */
1890 vars.alpha1 = 0.0;
1891 vars.alpha2 = 0.0;
1892 sys->maxstep = 0.0;
1893 sys->varstep.norm2 = 0.0;
1894
1895
1896 } else if( (tot2_norm2>0.0)&&(calc_sqrt_D0(tot2_norm2)<=sys->maxstep) ) {
1897 /* Attempt step in varstep2 direction */
1898 vars.alpha1 = 0.0;
1899 vars.alpha2 = 1.0;
1900 sys->maxstep = calc_sqrt_D0(tot2_norm2);
1901 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
1902 sys->varstep2.norm2/tot2_norm2;
1903
1904
1905 } else if( (tot2_norm2==0.0 || sys->s.block.current_size==1) &&
1906 (tot1_norm2 > 0.0) ) {
1907 /* Attempt step in varstep1 direction */
1908 vars.alpha1 = 1.0;
1909 vars.alpha2 = 0.0;
1910 if ( (sys->gamma.norm2/sys->Jgamma.norm2)*
1911 calc_sqrt_D0(sys->gamma.norm2) <= sys->maxstep )
1912 sys->maxstep = (sys->gamma.norm2/sys->Jgamma.norm2)*
1913 calc_sqrt_D0(sys->gamma.norm2);
1914 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
1915 sys->varstep1.norm2/tot1_norm2;
1916
1917
1918 } else {
1919 /* Attempt step in varstep1-varstep2 direction */
1920 vars.parms.low = 0.0;
1921 vars.parms.high = MAXDOUBLE;
1922 vars.parms.guess = 1.0;
1923 calc_2x2_system(sys,&vars);
1924 do {
1925 coefs_from_parm(sys, &vars);
1926 adjust_parms(sys, &vars);
1927 } while( fabs(vars.error) > STEPSIZEERR_MAX &&
1928 vars.parms.high - vars.parms.low > PARMRNG_MIN );
1929 if (SHOW_LESS_IMPT) {
1930 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1931 " parameter high", vars.parms.high);
1932 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1933 " parameter low", vars.parms.low);
1934 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1935 " Error in step length", vars.error);
1936 }
1937 sys->varstep.norm2 = step_norm2(sys, &vars);
1938
1939 }
1940
1941 if (SHOW_LESS_IMPT) {
1942 FPRINTF(LIF(sys),"%-40s ---> %g\n", " Alpha1 coefficient (normalized)",
1943 vars.alpha1);
1944 FPRINTF(LIF(sys),"%-40s ---> %g\n", " Alpha2 coefficient (normalized)",
1945 vars.alpha2);
1946 }
1947 compute_step(sys,&vars);
1948 return;
1949
1950 }
1951
1952
1953
1954 /*
1955 * Variable values maintenance
1956 * ---------------------------
1957 * restore_variables(sys)
1958 * coef = required_coef_to_stay_inbounds(sys)
1959 * apply_step(sys,coef)
1960 * step_accepted(sys)
1961 * change_maxstep(sys,maxstep)
1962 */
1963
1964 /*
1965 * Restores the values of the variables before applying
1966 * a step.
1967 */
1968 static void restore_variables( slv5_system_t sys)
1969 {
1970 int32 col;
1971 real64 *vec;
1972 vec = (sys->nominals.vec);
1973 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
1974 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1975 var_set_value(var,sys->variables.vec[col]*vec[col]);
1976 }
1977 }
1978
1979
1980 /*
1981 * Calculates the maximum fraction of the step which can be
1982 * taken without making negative the complementary vars.
1983 * It is assumed that the current complementary variable
1984 * is positive. The step must be calculated.
1985 */
1986 static real64 factor_for_complementary_vars( slv5_system_t sys, int32 v)
1987 {
1988 real64 factor, minfactor;
1989 struct var_variable *var;
1990 real64 dx,val,bnd;
1991 int32 col;
1992 struct vector_data step;
1993 real64 *vec;
1994
1995 vec = (sys->nominals.vec);
1996
1997 if (v == 1) {
1998 step = sys->varstep;
1999 } else {
2000 step = sys->varnewstep;
2001 }
2002
2003 minfactor = 1.0;
2004 factor = 1.0;
2005 for( col=step.rng->low; col <= step.rng->high; col++ ) {
2006 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2007 val = var_value(var);
2008 if (var_complementary(var) && (!var_fixed(var))) {
2009 dx = step.vec[col] * vec[col];
2010 bnd = 0.0;
2011 if( val + dx < bnd ) {
2012 factor = MIN((bnd-val)/dx, 1.0);
2013 if (factor < 1.0) {
2014 #if DEBUG_COMPLEMENTARY_VAR
2015 FPRINTF(ASCERR,"Negative Complementary Variable : \n");
2016 print_var_name(ASCERR,sys,var);
2017 FPRINTF(ASCERR,"\n");
2018 FPRINTF(ASCERR,"factor = %f \n",factor);
2019 #endif /* DEBUG_COMPLEMENTARY_VAR */
2020 }
2021 if( factor < minfactor ) {
2022 minfactor = factor;
2023 }
2024 }
2025 }
2026 }
2027 #if DEBUG_COMPLEMENTARY_VAR
2028 FPRINTF(ASCERR,"Minimum factor = %f \n",minfactor);
2029 #endif /* DEBUG_COMPLEMENTARY_VAR */
2030 return minfactor;
2031 }
2032
2033 /*
2034 * Calculates the maximum fraction of the step which can be
2035 * taken without going out of bounds. If the entire step can be
2036 * taken, 1.0 is returned. Otherwise a value less than 1 is
2037 * returned. It is assumed that the current variable values
2038 * are within their bounds. The step must be calculated.
2039 */
2040 static real64 required_coef_to_stay_inbounds( slv5_system_t sys, int32 v)
2041 {
2042 real64 mincoef;
2043 int32 col;
2044 real64 *vec;
2045 struct vector_data step;
2046 vec = (sys->nominals.vec);
2047
2048 if( sys->p.ignore_bounds )
2049 return(1.0);
2050
2051 if (v == 1) {
2052 step = sys->varstep;
2053 } else {
2054 step = sys->varnewstep;
2055 }
2056
2057 mincoef = 1.0;
2058 for( col=step.rng->low; col <= step.rng->high; col++ ) {
2059 struct var_variable *var;
2060 real64 coef,dx,val,bnd;
2061 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2062 coef = 1.0;
2063 dx = step.vec[col] * vec[col];
2064 bnd = var_upper_bound(var);
2065 if( (val=var_value(var)) + dx > bnd )
2066 coef = MIN((bnd-val)/dx, 1.0);
2067 bnd = var_lower_bound(var);
2068 if( val + dx < bnd )
2069 coef = MIN((bnd-val)/dx, 1.0);
2070 if( coef < mincoef )
2071 mincoef = coef;
2072 }
2073 #if DEBUG_COMPLEMENTARY_VAR
2074 FPRINTF(ASCERR,"Minimum coefficient to stay in bounds = %f \n",mincoef);
2075 #endif /* DEBUG_COMPLEMENTARY_VAR */
2076 return(mincoef);
2077 }
2078
2079 /*
2080 * Adds sys->varnewstep to the variable values in block: projecting
2081 * near bounds.
2082 */
2083 static void apply_newton_step( slv5_system_t sys)
2084 {
2085 FILE *lif = LIF(sys);
2086 int nproj = 0;
2087 real64 bounds_coef = 1.0;
2088 int32 col;
2089 real64 *vec;
2090 real64 factor;
2091
2092 vec = (sys->nominals.vec);
2093 factor = factor_for_complementary_vars(sys,0);
2094 if (factor >= 0.9995) {
2095 sys->sigma_const = sys->sigma_const/10.0;
2096 }
2097 factor = 0.995 * factor;
2098
2099 if (TRUNCATE && (!sys->p.ignore_bounds)) {
2100 bounds_coef = required_coef_to_stay_inbounds(sys,0);
2101 }
2102
2103 if ((bounds_coef < 1.0) && (factor < bounds_coef)) {
2104 bounds_coef = factor;
2105 }
2106
2107 for(col=sys->varnewstep.rng->low; col<=sys->varnewstep.rng->high;col++) {
2108 struct var_variable *var;
2109 real64 dx,val,bnd;
2110 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2111 dx = vec[col]*sys->varnewstep.vec[col];
2112 val = var_value(var);
2113 if (bounds_coef < 1.0) {
2114 dx = dx*TOWARD_BOUNDS*bounds_coef;
2115 sys->varnewstep.vec[col] = dx/vec[col];
2116 } else {
2117 if ((var_complementary(var)) && (!var_fixed(var))) {
2118 dx = dx*factor;
2119 sys->varnewstep.vec[col] = dx/vec[col];
2120 } else {
2121 if( !sys->p.ignore_bounds ) {
2122 if( val + dx > (bnd=var_upper_bound(var)) ) {
2123 dx = TOWARD_BOUNDS*(bnd-val);
2124 sys->varnewstep.vec[col] = dx/vec[col];
2125 if (SHOW_LESS_IMPT) {
2126 FPRINTF(lif,"%-40s ---> ",
2127 " Variable projected to upper bound");
2128 print_var_name(lif,sys,var); PUTC('\n',lif);
2129 }
2130 ++nproj;
2131 } else if( val + dx < (bnd=var_lower_bound(var)) ) {
2132 dx = TOWARD_BOUNDS*(bnd-val);
2133 sys->varnewstep.vec[col] = dx/vec[col];
2134 if (SHOW_LESS_IMPT) {
2135 FPRINTF(lif,"%-40s ---> ",
2136 " Variable projected to lower bound");
2137 print_var_name(lif,sys,var); PUTC('\n',lif);
2138 }
2139 ++nproj;
2140 }
2141 }
2142 }
2143 }
2144 var_set_value(var,val+dx);
2145 }
2146
2147 /* Allow weighted residuals to be recalculated at new point */
2148 sys->newton_residuals.accurate = FALSE;
2149 sys->perturbed_residuals.accurate = FALSE;
2150 sys->residuals.accurate = FALSE;
2151
2152 return;
2153 }
2154
2155
2156 /*
2157 * Adds sys->varstep to the variable values in block: projecting
2158 * near bounds.
2159 */
2160 static void apply_step( slv5_system_t sys)
2161 {
2162 FILE *lif = LIF(sys);
2163 int nproj = 0;
2164 real64 bounds_coef = 1.0;
2165 int32 col;
2166 real64 *vec;
2167 real64 factor;
2168 vec = (sys->nominals.vec);
2169
2170 if (TRUNCATE && (!sys->p.ignore_bounds))
2171 bounds_coef = required_coef_to_stay_inbounds(sys,1);
2172
2173 for( col=sys->varstep.rng->low; col <= sys->varstep.rng->high; col++ ) {
2174 struct var_variable *var;
2175 real64 dx,val,bnd;
2176 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2177 dx = vec[col]*sys->varstep.vec[col];
2178 val = var_value(var);
2179 if (bounds_coef < 1.0) {
2180 dx = dx*TOWARD_BOUNDS*bounds_coef;
2181 sys->varstep.vec[col] = dx/vec[col];
2182 } else {
2183 if( !sys->p.ignore_bounds ) {
2184 if( val + dx > (bnd=var_upper_bound(var)) ) {
2185 dx = TOWARD_BOUNDS*(bnd-val);
2186 sys->varstep.vec[col] = dx/vec[col];
2187 if (SHOW_LESS_IMPT) {
2188 FPRINTF(lif,"%-40s ---> ",
2189 " Variable projected to upper bound");
2190 print_var_name(lif,sys,var); PUTC('\n',lif);
2191 }
2192 ++nproj;
2193 } else if( val + dx < (bnd=var_lower_bound(var)) ) {
2194 dx = TOWARD_BOUNDS*(bnd-val);
2195 sys->varstep.vec[col] = dx/vec[col];
2196 if (SHOW_LESS_IMPT) {
2197 FPRINTF(lif,"%-40s ---> ",
2198 " Variable projected to lower bound");
2199 print_var_name(lif,sys,var); PUTC('\n',lif);
2200 }
2201 ++nproj;
2202 }
2203 }
2204 }
2205 var_set_value(var,val+dx);
2206 }
2207
2208 if( !sys->p.ignore_bounds ) {
2209 if (nproj > 0) {
2210 square_norm(&(sys->varstep));
2211 sys->progress = calc_sqrt_D0
2212 (calc_sqrt_D0((sys->varstep.norm2 )*
2213 (sys->varstep1.norm2 )));
2214 if (SHOW_LESS_IMPT) {
2215 FPRINTF(lif,"%-40s ---> %g\n", " Projected step length (scaled)",
2216 calc_sqrt_D0(sys->varstep.norm2));
2217 FPRINTF(lif,"%-40s ---> %g\n",
2218 " Projected progress", sys->progress);
2219 }
2220 }
2221 if (bounds_coef < 1.0) {
2222 square_norm(&(sys->varstep));
2223 if (SHOW_LESS_IMPT) {
2224 FPRINTF(lif,"%-40s ---> %g\n",
2225 " Truncated step length (scaled)",
2226 calc_sqrt_D0(sys->varstep.norm2 ));
2227 }
2228 sys->progress = calc_sqrt_D0
2229 (calc_sqrt_D0((sys->varstep.norm2 )*
2230 (sys->varstep1.norm2 )));
2231 (calc_sqrt_D0(sys->varstep.norm2*sys->varstep1.norm2));
2232 if (SHOW_LESS_IMPT) {
2233 FPRINTF(lif,"%-40s ---> %g\n",
2234 " Truncated progress", sys->progress);
2235 }
2236 }
2237 }
2238
2239 /* Allow weighted residuals to be recalculated at new point */
2240 sys->newton_residuals.accurate = FALSE;
2241 sys->perturbed_residuals.accurate = FALSE;
2242 sys->residuals.accurate = FALSE;
2243 return;
2244 }
2245
2246
2247 /*
2248 * This function should be called when the step is accepted.
2249 */
2250 static void step_accepted( slv5_system_t sys)
2251 {
2252 /* Maintain update status on jacobian and weights */
2253 if (--(sys->update.jacobian) <= 0) {
2254 sys->J.accurate = FALSE;
2255 }
2256
2257 sys->variables.accurate = FALSE;
2258 sys->newton_residuals.accurate = FALSE;
2259 sys->perturbed_residuals.accurate = FALSE;
2260 sys->newton.accurate = FALSE;
2261 sys->perturbed_newton.accurate = FALSE;
2262 sys->gamma.accurate = FALSE;
2263 sys->Jgamma.accurate = FALSE;
2264 sys->varstep1.accurate = FALSE;
2265 sys->varnewstep.accurate = FALSE;
2266 sys->varstep2.accurate = FALSE;
2267 sys->varstep.accurate = FALSE;
2268 }
2269
2270 /*
2271 * This function changes sys->maxstep to the given number and should be
2272 * called whenever sys->maxstep is to be changed.
2273 */
2274 static void change_maxstep( slv5_system_t sys, real64 maxstep)
2275 {
2276 sys->maxstep = maxstep;
2277 sys->varstep.accurate = FALSE;
2278 }
2279
2280
2281 /*
2282 * Block routines
2283 * --------------
2284 * feas = block_feasible(sys)
2285 * move_to_next_block(sys)
2286 * find_next_unconverged_block(sys)
2287 */
2288
2289 /*
2290 * Returns TRUE if the current block is feasible, FALSE otherwise.
2291 * It is assumed that the residuals have been computed.
2292 */
2293 static boolean block_feasible( slv5_system_t sys)
2294 {
2295 int32 row;
2296
2297 if( !sys->s.calc_ok )
2298 return(FALSE);
2299
2300 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2301 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2302 if( !rel_satisfied(rel) ) return FALSE;
2303 }
2304 return TRUE;
2305 }
2306
2307 /*
2308 * Moves on to the next block, updating all of the solver information.
2309 * To move to the first block, set sys->s.block.current_block to -1 before
2310 * calling. If already at the last block, then sys->s.block.current_block
2311 * will equal the number of blocks and the system will be declared
2312 * converged. Otherwise, the residuals for the new block will be computed
2313 * and sys->s.calc_ok set according.
2314 */
2315 static void move_to_next_block( slv5_system_t sys)
2316 {
2317 struct var_variable *var;
2318 struct rel_relation *rel;
2319 int32 row;
2320 int32 col;
2321 int32 ci;
2322
2323 if( sys->s.block.current_block >= 0 ) {
2324
2325 /* Record cost accounting info here. */
2326 ci=sys->s.block.current_block;
2327 sys->s.cost[ci].size = sys->s.block.current_size;
2328 sys->s.cost[ci].iterations = sys->s.block.iteration;
2329 sys->s.cost[ci].funcs = sys->s.block.funcs;
2330 sys->s.cost[ci].jacs = sys->s.block.jacs;
2331 sys->s.cost[ci].functime = sys->s.block.functime;
2332 sys->s.cost[ci].jactime = sys->s.block.jactime;
2333 sys->s.cost[ci].time = sys->s.block.cpu_elapsed;
2334 sys->s.cost[ci].resid = sys->s.block.residual;
2335
2336 /* De-initialize previous block */
2337 if (SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
2338 LIFDS)) {
2339 FPRINTF(LIF(sys),"Block %d converged.\n",
2340 sys->s.block.current_block);
2341 }
2342 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
2343 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2344 var_set_in_block(var,FALSE);
2345 }
2346 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2347 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2348 rel_set_in_block(rel,FALSE);
2349 }
2350 sys->s.block.previous_total_size += sys->s.block.current_size;
2351 }
2352
2353 sys->s.block.current_block++;
2354 if( sys->s.block.current_block < sys->s.block.number_of ) {
2355 boolean ok;
2356
2357 /* Initialize next block */
2358
2359 sys->J.reg =
2360 (slv_get_solvers_blocks(SERVER))->block[sys->s.block.current_block];
2361
2362
2363 row = sys->J.reg.row.high - sys->J.reg.row.low + 1;
2364 col = sys->J.reg.col.high - sys->J.reg.col.low + 1;
2365 sys->s.block.current_size = MAX(row,col);
2366
2367 sys->s.block.iteration = 0;
2368 sys->s.block.cpu_elapsed = 0.0;
2369 sys->s.block.functime = 0.0;
2370 sys->s.block.jactime = 0.0;
2371 sys->s.block.funcs = 0;
2372 sys->s.block.jacs = 0;
2373
2374 if(SHOW_LESS_IMPT && (LIFDS ||
2375 sys->s.block.current_size > 1)) {
2376 debug_delimiter(LIF(sys));
2377 debug_delimiter(LIF(sys));
2378 }
2379 if(SHOW_LESS_IMPT && LIFDS) {
2380 FPRINTF(LIF(sys),"\n%-40s ---> %d in [%d..%d]\n",
2381 "Current block number", sys->s.block.current_block,
2382 0, sys->s.block.number_of-1);
2383 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Current block size",
2384 sys->s.block.current_size);
2385 }
2386 sys->s.calc_ok = TRUE;
2387
2388 if( !(ok = calc_objective(sys)) ) {
2389 FPRINTF(MIF(sys),"Objective calculation errors detected.\n");
2390 }
2391 if(SHOW_LESS_IMPT && sys->obj) {
2392 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Objective", sys->objective);
2393 }
2394 sys->s.calc_ok = sys->s.calc_ok && ok;
2395
2396 if (!(sys->p.ignore_bounds) ) {
2397 slv_insure_bounds(SERVER, sys->J.reg.col.low,
2398 sys->J.reg.col.high,MIF(sys));
2399 }
2400
2401 sys->residuals.accurate = FALSE;
2402 if( !(ok = calc_residuals(sys)) ) {
2403 FPRINTF(MIF(sys),
2404 "Residual calculation errors detected in move_to_next_block.\n");
2405 }
2406
2407 /*
2408 * Update number of complementary equations
2409 */
2410 calc_comp(sys);
2411 sys->sigma_const = SIGMA;
2412
2413 if( SHOW_LESS_IMPT &&
2414 (sys->s.block.current_size >1 ||
2415 LIFDS) ) {
2416 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
2417 sys->s.block.residual);
2418 }
2419 sys->s.calc_ok = sys->s.calc_ok && ok;
2420
2421 /* Must be updated as soon as required */
2422 sys->J.accurate = FALSE;
2423 sys->update.weights = 0;
2424 sys->update.nominals = 0;
2425 sys->update.relnoms = 0;
2426 sys->update.iterative = 0;
2427 sys->variables.accurate = FALSE;
2428 sys->newton_residuals.accurate = FALSE;
2429 sys->perturbed_residuals.accurate = FALSE;
2430 sys->newton.accurate = FALSE;
2431 sys->perturbed_newton.accurate = FALSE;
2432 sys->gamma.accurate = FALSE;
2433 sys->Jgamma.accurate = FALSE;
2434 sys->varstep1.accurate = FALSE;
2435 sys->varstep2.accurate = FALSE;
2436 sys->varnewstep.accurate = FALSE;
2437 sys->varstep.accurate = FALSE;
2438
2439 } else {
2440 boolean ok;
2441 /*
2442 * Before we claim convergence, we must check if we left behind
2443 * some unassigned relations. If and only if they happen to be
2444 * satisfied at the current point, convergence has been obtained.
2445 *
2446 * Also insures that all included relations have valid residuals.
2447 * Included inequalities will have correct residuals.
2448 * Unsatisfied included inequalities cause inconsistency.
2449 *
2450 * This of course ignores that fact an objective function might
2451 * be present. Then, feasibility isn't enough, is it now.
2452 */
2453 if( sys->s.struct_singular ) {
2454 /* black box w/singletons provoking bug here, maybe */
2455 sys->s.block.current_size = sys->rused - sys->rank;
2456 if(SHOW_LESS_IMPT) {
2457 debug_delimiter(LIF(sys));
2458 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Unassigned Relations",
2459 sys->s.block.current_size);
2460 }
2461 sys->J.reg.row.low = sys->J.reg.col.low = sys->rank;
2462 sys->J.reg.row.high = sys->J.reg.col.high = sys->rused - 1;
2463 sys->residuals.accurate = FALSE;
2464 if( !(ok=calc_residuals(sys)) ) {
2465 FPRINTF(MIF(sys),
2466 "Residual calculation errors detected in leftover equations.\n");
2467 }
2468 if(SHOW_LESS_IMPT) {
2469 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
2470 sys->s.block.residual);
2471 }
2472 if( block_feasible(sys) ) {
2473 if(SHOW_LESS_IMPT) {
2474 FPRINTF(LIF(sys),"\nUnassigned relations ok. You lucked out.\n");
2475 }
2476 sys->s.converged = TRUE;
2477 } else {
2478 if(SHOW_LESS_IMPT) {
2479 FPRINTF(LIF(sys),"\nProblem inconsistent: %s.\n",
2480 "Unassigned relations not satisfied");
2481 }
2482 sys->s.inconsistent = TRUE;
2483 }
2484 if(SHOW_LESS_IMPT) {
2485 debug_delimiter(LIF(sys));
2486 }
2487 } else {
2488 sys->s.converged = TRUE;
2489 }
2490 /* nearly done checking. Must verify included inequalities if
2491 we think equalities are ok. */
2492 if (sys->s.converged) {
2493 sys->s.inconsistent=(!calc_inequalities(sys));
2494 }
2495 }
2496 }
2497
2498
2499 /*
2500 * Calls the appropriate reorder function on a block
2501 */
2502 static void reorder_new_block(slv5_system_t sys)
2503 {
2504 int32 method;
2505 if( sys->s.block.current_block < sys->s.block.number_of ) {
2506 if (strcmp(REORDER_OPTION,"SPK1") == 0) {
2507 method = 2;
2508 } else {
2509 method = 1;
2510 }
2511
2512 if( sys->s.block.current_block <= sys->s.block.current_reordered_block &&
2513 sys->s.cost[sys->s.block.current_block].reorder_method == method &&
2514 sys->s.block.current_block >= 0 ) {
2515 #if DEBUG
2516 FPRINTF(stderr,"YOU JUST AVOIDED A REORDERING\n");
2517 #endif
2518 slv_set_up_block(SERVER,sys->s.block.current_block);
2519 /* tell linsol to bless it and get on with things */
2520 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
2521 return; /*must have been reordered since last system build*/
2522 }
2523
2524 /* Let the slv client function take care of reordering things
2525 * and setting in block flags.
2526 */
2527 if (strcmp(REORDER_OPTION,"SPK1") == 0) {
2528 sys->s.cost[sys->s.block.current_block].reorder_method = 2;
2529 slv_spk1_reorder_block(SERVER,sys->s.block.current_block,1);
2530 } else if (strcmp(REORDER_OPTION,"TEAR_DROP") == 0) {
2531 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2532 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2533 CUTOFF,
2534 0,mtx_SPK1);
2535 /* khack: try tspk1 for transpose case */
2536 } else if (strcmp(REORDER_OPTION,"OVER_TEAR") == 0) {
2537 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2538 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2539 CUTOFF,
2540 1,mtx_SPK1);
2541 } else {
2542 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2543 FPRINTF(MIF(sys),"IPSlv called with unknown reorder option\n");
2544 FPRINTF(MIF(sys),"IPSlv using single edge tear drop (TEAR_DROP).\n");
2545 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2546 CUTOFF,0,mtx_SPK1);
2547 }
2548 /* tell linsol to bless it and get on with things */
2549 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
2550 if (sys->s.block.current_block > sys->s.block.current_reordered_block) {
2551 sys->s.block.current_reordered_block = sys->s.block.current_block;
2552 }
2553 }
2554 }
2555
2556 /*
2557 * Moves to next unconverged block, assuming that the current block has
2558 * converged (or is -1, to start).
2559 */
2560 static void find_next_unconverged_block( slv5_system_t sys)
2561 {
2562 do {
2563 move_to_next_block(sys);
2564 #if DEBUG
2565 debug_out_var_values(stderr,sys);
2566 debug_out_rel_residuals(stderr,sys);
2567 #endif
2568 } while( !sys->s.converged && block_feasible(sys));
2569 reorder_new_block(sys);
2570 }
2571
2572
2573 /*
2574 * Iteration begin/end routines
2575 * ----------------------------
2576 * iteration_begins(sys)
2577 * iteration_ends(sys)
2578 */
2579
2580 /*
2581 * Prepares sys for entering an iteration, increasing the iteration counts
2582 * and starting the clock.
2583 */
2584 static void iteration_begins( slv5_system_t sys)
2585 {
2586 sys->clock = tm_cpu_time();
2587 ++(sys->s.block.iteration);
2588 ++(sys->s.iteration);
2589 if(SHOW_LESS_IMPT&& (sys->s.block.current_size >1 ||
2590 LIFDS)) {
2591 FPRINTF(LIF(sys),"\n%-40s ---> %d\n",
2592 "Iteration", sys->s.block.iteration);
2593 FPRINTF(LIF(sys),"%-40s ---> %d\n",
2594 "Total iteration", sys->s.iteration);
2595 }
2596 }
2597
2598 /*
2599 * Prepares sys for exiting an iteration, stopping the clock and recording
2600 * the cpu time.
2601 */
2602 static void iteration_ends( slv5_system_t sys)
2603 {
2604 double cpu_elapsed; /* elapsed this iteration */
2605
2606 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
2607 sys->s.block.cpu_elapsed += cpu_elapsed;
2608 sys->s.cpu_elapsed += cpu_elapsed;
2609 if(SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
2610 LIFDS)) {
2611 FPRINTF(LIF(sys),"%-40s ---> %g\n",
2612 "Elapsed time", sys->s.block.cpu_elapsed);
2613 FPRINTF(LIF(sys),"%-40s ---> %g\n",
2614 "Total elapsed time", sys->s.cpu_elapsed);
2615 }
2616 }
2617
2618 /*
2619 * Updates the solver status.
2620 */
2621 static void update_status( slv5_system_t sys)
2622 {
2623 boolean unsuccessful;
2624
2625 if( !sys->s.converged ) {
2626 sys->s.time_limit_exceeded =
2627 (sys->s.block.cpu_elapsed >= TIME_LIMIT);
2628 sys->s.iteration_limit_exceeded =
2629 (sys->s.block.iteration >= ITER_LIMIT);
2630 }
2631
2632 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
2633 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
2634
2635 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
2636 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
2637 }
2638
2639 static
2640 int32 slv5_get_default_parameters(slv_system_t server, SlvClientToken asys,
2641 slv_parameters_t *parameters)
2642 {
2643 slv5_system_t sys;
2644 union parm_arg lo,hi,val;
2645 struct slv_parameter *new_parms = NULL;
2646 int32 make_macros = 0;
2647
2648 static char *factor_names[] = {
2649 "SPK1/RANKI","SPK1/RANKI+ROW",
2650 "Fast-SPK1/RANKI","Fast-SPK1/RANKI+ROW",
2651 "Fastest-SPK1/MR-RANKI","CondQR","CPQR"
2652 /* ,"GAUSS","GAUSS_EASY" currently only works for ken */
2653 };
2654 static char *reorder_names[] = {
2655 "SPK1","TEAR_DROP","OVER_TEAR"
2656 };
2657 static char *converge_names[] = {
2658 "ABSOLUTE","RELNOM_SCALE"
2659 };
2660 static char *scaling_names[] = {
2661 "NONE","ROW_2NORM","RELNOM"
2662 };
2663
2664 if (server != NULL && asys != NULL) {
2665 sys = SLV5(asys);
2666 make_macros = 1;
2667 }
2668
2669 #ifndef NDEBUG /* keep purify from whining on UMR */
2670 lo.argr = hi.argr = val.argr = 0.0;
2671 #endif
2672
2673 if (parameters->parms == NULL) {
2674 /* an external client wants our parameter list.
2675 * an instance of slv5_system_structure has this pointer
2676 * already set in slv5_create
2677 */
2678 new_parms = (struct slv_parameter *)
2679 ascmalloc((slv5_PA_SIZE)*sizeof(struct slv_parameter));
2680 if (new_parms == NULL) {
2681 return -1;
2682 }
2683 parameters->parms = new_parms;
2684 parameters->dynamic_parms = 1;
2685 }
2686 parameters->num_parms = 0;
2687
2688 /* begin defining parameters */
2689
2690 slv_define_parm(parameters, bool_parm,
2691 "ignorebounds","ignore bounds?","ignore bounds?",
2692 U_p_bool(val,0),U_p_bool(lo,0),U_p_bool(hi,1),-1);
2693 SLV_BPARM_MACRO(IGNORE_BOUNDS_PTR,parameters);
2694
2695 slv_define_parm(parameters, bool_parm,
2696 "showmoreimportant", "showmoreimportant", "showmoreimportant",
2697 U_p_bool(val,1),U_p_bool(lo,0),U_p_bool(hi,1),-1);
2698 SLV_BPARM_MACRO(SHOW_MORE_IMPT_PTR,parameters);
2699
2700 slv_define_parm(parameters, real_parm,
2701 "rho", "penalty parameter", "penalty parameter",
2702 U_p_real(val,100),U_p_real(lo, 0),U_p_real(hi,10e100), 1);
2703 SLV_RPARM_MACRO(RHO_PTR,parameters);
2704
2705 slv_define_parm(parameters, bool_parm,
2706 "partition", "partitioning enabled", "partitioning enabled",
2707 U_p_bool(val, 1),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2708 SLV_BPARM_MACRO(PARTITION_PTR,parameters);
2709
2710 slv_define_parm(parameters, bool_parm,
2711 "sigmaconstant", "sigma constant", "sigma constant",
2712 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2713 SLV_BPARM_MACRO(SIGMA_CONSTANT_PTR,parameters);
2714
2715 slv_define_parm(parameters, bool_parm,
2716 "mumaximum", "mu maximum", "mu maximum",
2717 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2718 SLV_BPARM_MACRO(MU_MAXIMUM_PTR,parameters);
2719
2720 slv_define_parm(parameters, bool_parm,
2721 "showlessimportant", "detailed solving info",
2722 "detailed solving info",
2723 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2724 SLV_BPARM_MACRO(SHOW_LESS_IMPT_PTR,parameters);
2725
2726 slv_define_parm(parameters, bool_parm,
2727 "autoresolve", "auto-resolve", "auto-resolve",
2728 U_p_bool(val,1),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2729 SLV_BPARM_MACRO(AUTO_RESOLVE_PTR,parameters);
2730
2731 slv_define_parm(parameters, int_parm,
2732 "timelimit", "time limit (CPU sec/block)", "time limit (CPU sec/block)",
2733 U_p_int(val,1500),U_p_int(lo, 1),U_p_int(hi,20000),1);
2734 SLV_IPARM_MACRO(TIME_LIMIT_PTR,parameters);
2735
2736 slv_define_parm(parameters, int_parm,
2737 "iterationlimit", "max iterations/block", "max iterations/block",
2738 U_p_int(val, 30),U_p_int(lo, 1),U_p_int(hi,20000),1);
2739 SLV_IPARM_MACRO(ITER_LIMIT_PTR,parameters);
2740
2741 slv_define_parm(parameters,real_parm,
2742 "stattol","stattol" ,"stattol" ,
2743 U_p_real(val,1e-6),U_p_real(lo,0),U_p_real(hi,1.0),-1);
2744 SLV_RPARM_MACRO(STAT_TOL_PTR,parameters);
2745
2746 slv_define_parm(parameters,real_parm,
2747 "termtol","termtol" ,"termtol" ,
2748 U_p_real(val,1e-12),U_p_real(lo,0),U_p_real(hi,1.0),-1);
2749 SLV_RPARM_MACRO(TERM_TOL_PTR,parameters);
2750
2751 slv_define_parm(parameters, real_parm,
2752 "singtol", "epsilon (min pivot)", "epsilon (min pivot)",
2753 U_p_real(val, 1e-12),U_p_real(lo, 1e-12),U_p_real(hi,1.0),1);
2754 SLV_RPARM_MACRO(SING_TOL_PTR,parameters);
2755
2756 slv_define_parm(parameters, real_parm,
2757 "sigmaconst", "sigma (centering parameter)",
2758 "sigma (centering parameter)",
2759 U_p_real(val, 0.5),U_p_real(lo, 0.0),U_p_real(hi,1.0),1);
2760 SLV_RPARM_MACRO(SIGMA_PTR,parameters);
2761
2762 slv_define_parm(parameters, real_parm,
2763 "mumax", "mu maximum (centering parameter)",
2764 "mu maximum (centering parameter)",
2765 U_p_real(val, 0.9),U_p_real(lo, 0.0),U_p_real(hi,1.0),1);
2766 SLV_RPARM_MACRO(MUK_PTR,parameters);
2767
2768 slv_define_parm(parameters, real_parm,
2769 "pivottol", "condition tolerance", "condition tolerance",
2770 U_p_real(val, 0.5),U_p_real(lo, 0),U_p_real(hi, 1),1);
2771 SLV_RPARM_MACRO(PIVOT_TOL_PTR,parameters);
2772
2773 slv_define_parm(parameters, real_parm,
2774 "feastol", "max residual (absolute)", "max residual (absolute)",
2775 U_p_real(val, 1e-8),U_p_real(lo, 1e-13),U_p_real(hi,100000.5),1);
2776 SLV_RPARM_MACRO(FEAS_TOL_PTR,parameters);
2777
2778 slv_define_parm(parameters, bool_parm,
2779 "lifds", "show singletons details", IEX(0),
2780 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2781 SLV_BPARM_MACRO(LIFDS_PTR,parameters);
2782
2783 slv_define_parm(parameters, bool_parm,
2784 "savlin", "write to file SlvLinsol.dat", IEX(1),
2785 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2786 SLV_BPARM_MACRO(SAVLIN_PTR,parameters);
2787
2788 slv_define_parm(parameters, bool_parm,
2789 "safe_calc", "safe calculations", IEX(12),
2790 U_p_bool(val, 1),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2791
2792 SLV_BPARM_MACRO(SAFE_CALC_PTR,parameters);
2793
2794 slv_define_parm(parameters, bool_parm,
2795 "relnomscale", "calc rel nominals", IEX(2),
2796 U_p_bool(val, 1),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2797 SLV_BPARM_MACRO(RELNOMSCALE_PTR,parameters);
2798
2799 slv_define_parm(parameters, int_parm,
2800 "cutoff", "block size cutoff (MODEL-based)", IEX(3),
2801 U_p_int(val, 500),U_p_int(lo,0),U_p_int(hi,20000), 2);
2802 SLV_IPARM_MACRO(CUTOFF_PTR,parameters);
2803
2804
2805 slv_define_parm(parameters, int_parm,
2806 "upjac", "Jacobian update frequency", IEX(4),
2807 U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
2808 SLV_IPARM_MACRO(UPDATE_JACOBIAN_PTR,parameters);
2809
2810 slv_define_parm(parameters, int_parm,
2811 "upwts", "Row scaling update frequency", IEX(5),
2812 U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
2813 SLV_IPARM_MACRO(UPDATE_WEIGHTS_PTR,parameters);
2814
2815 slv_define_parm(parameters, int_parm,
2816 "upnom", "Column scaling update frequency", IEX(6),
2817 U_p_int(val, 1000),U_p_int(lo,0),U_p_int(hi,20000), 3);
2818 SLV_IPARM_MACRO(UPDATE_NOMINALS_PTR,parameters);
2819
2820 slv_define_parm(parameters, int_parm,
2821 "uprelnom", "Relation nominal update frequency", IEX(13),
2822 U_p_int(val, 5),U_p_int(lo,0),U_p_int(hi,20000), 3);
2823 SLV_IPARM_MACRO(UPDATE_RELNOMS_PTR,parameters);
2824
2825 slv_define_parm(parameters, int_parm,
2826 "itscalelim", "Iteration lim for iterative scale", IEX(14),
2827 U_p_int(val, 10),U_p_int(lo,0),U_p_int(hi,20000), 3);
2828 SLV_IPARM_MACRO(ITSCALELIM_PTR,parameters);
2829
2830 slv_define_parm(parameters, char_parm,
2831 "convopt", "convergence test", "convergence test",
2832 U_p_string(val,converge_names[0]),
2833 U_p_strings(lo,converge_names),
2834 U_p_int(hi,sizeof(converge_names)/sizeof(char *)),1);
2835 SLV_CPARM_MACRO(CONVOPT_PTR,parameters);
2836
2837 slv_define_parm(parameters, char_parm,
2838 "scaleopt", "jacobian scaling option", IEX(15),
2839 U_p_string(val,scaling_names[1]),
2840 U_p_strings(lo,scaling_names),
2841 U_p_int(hi,sizeof(scaling_names)/sizeof(char *)),1);
2842 SLV_CPARM_MACRO(SCALEOPT_PTR,parameters);
2843
2844 slv_define_parm(parameters, bool_parm,
2845 "reduce", "step reduction on?", IEX(7),
2846 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2847 SLV_BPARM_MACRO(REDUCE_PTR,parameters);
2848
2849 slv_define_parm(parameters, bool_parm,
2850 "exact", "exact line search", IEX(8),
2851 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2852 SLV_BPARM_MACRO(EXACT_LINE_SEARCH_PTR,parameters);
2853
2854 slv_define_parm(parameters, bool_parm,
2855 "cncols", "Check poorly scaled columns", IEX(9),
2856 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2857 SLV_BPARM_MACRO(DUMPCNORM_PTR,parameters);
2858
2859 slv_define_parm(parameters, bool_parm,
2860 "lintime", "Enable linsolqr timing",
2861 "Enable linsolqr factor timing",
2862 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2863 SLV_BPARM_MACRO(LINTIME_PTR,parameters);
2864
2865 slv_define_parm(parameters, bool_parm,
2866 "btrunc", "truncate whole step vector", IEX(10),
2867 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2868 SLV_BPARM_MACRO(TRUNCATE_PTR,parameters);
2869
2870 slv_define_parm(parameters, char_parm,
2871 "reorder", "reorder method", IEX(11),
2872 U_p_string(val,reorder_names[0]),
2873 U_p_strings(lo,reorder_names),
2874 U_p_int(hi,sizeof(reorder_names)/sizeof(char *)),1);
2875 SLV_CPARM_MACRO(REORDER_OPTION_PTR,parameters);
2876
2877
2878 slv_define_parm(parameters, real_parm,
2879 "toosmall", "default for zero nominal", REX(0),
2880 U_p_real(val, 1e-8),U_p_real(lo, 1e-12),U_p_real(hi,1.5), 3);
2881 SLV_RPARM_MACRO(TOO_SMALL_PTR,parameters);
2882
2883 slv_define_parm(parameters, real_parm,
2884 "cnlow", "smallest allowable column norm", REX(1),
2885 U_p_real(val, 0.01),U_p_real(lo, 0),U_p_real(hi,100000000.5), 3);
2886 SLV_RPARM_MACRO(CNLOW_PTR,parameters);
2887
2888 slv_define_parm(parameters, real_parm,
2889 "cnhigh", "largest allowable column norm", REX(2),
2890 U_p_real(val, 100.0),U_p_real(lo,0),U_p_real(hi,100000000.5), 3);
2891 SLV_RPARM_MACRO(CNHIGH_PTR,parameters);
2892
2893 slv_define_parm(parameters, real_parm,
2894 "tobnds", "fraction move to bounds", REX(3),
2895 U_p_real(val, 0.95),U_p_real(lo, 0),U_p_real(hi,1.0), 3);
2896 SLV_RPARM_MACRO(TOWARD_BOUNDS_PTR,parameters);
2897
2898 slv_define_parm(parameters, real_parm,
2899 "posdef", "Positive Definite Hessian Check", REX(4),
2900 U_p_real(val, 0.01),U_p_real(lo,0),U_p_real(hi,1.0), 3);
2901 SLV_RPARM_MACRO(POSITIVE_DEFINITE_PTR,parameters);
2902
2903 slv_define_parm(parameters, real_parm,
2904 "detzero", "Min newt/grad determinant ||", REX(5),
2905 U_p_real(val, 1e-8),U_p_real(lo,0),U_p_real(hi,1.0), 3);
2906 SLV_RPARM_MACRO(DETZERO_PTR,parameters);
2907
2908 slv_define_parm(parameters, real_parm,
2909 "steperrmax", "Step size precision", REX(6),
2910 U_p_real(val, 1e-4),U_p_real(lo, 1e-10),U_p_real(hi,1.0), 3);
2911 SLV_RPARM_MACRO(STEPSIZEERR_MAX_PTR,parameters);
2912
2913 slv_define_parm(parameters, real_parm,
2914 "prngmin", "Parameter range tolerance (in-loop)", REX(7),
2915 U_p_real(val, 1e-12),U_p_real(lo, 1e-12),U_p_real(hi,1.0), 3);
2916 SLV_RPARM_MACRO(PARMRNG_MIN_PTR,parameters);
2917
2918 slv_define_parm(parameters, real_parm,
2919 "mincoef", "'Largest' drop in maxstep allowed", REX(8),
2920 U_p_real(val, 0.05),U_p_real(lo, 1e-5),U_p_real(hi,1.0), 3);
2921 SLV_RPARM_MACRO(MIN_COEF_PTR,parameters);
2922
2923 slv_define_parm(parameters, real_parm,
2924 "maxcoef", "'Smallest' drop in maxstep allowed", REX(9),
2925 U_p_real(val, 0.99999),U_p_real(lo,0),U_p_real(hi,1.0), 3);
2926 SLV_RPARM_MACRO(MAX_COEF_PTR,parameters);
2927
2928 slv_define_parm(parameters, real_parm,
2929 "itscaletol", "Iterative scaling tolerance", REX(10),
2930 U_p_real(val, 0.99999),U_p_real(lo,0),U_p_real(hi,1.0), 3);
2931 SLV_RPARM_MACRO(ITSCALETOL_PTR,parameters);
2932
2933 slv_define_parm(parameters, char_parm,
2934 "bppivoting","linear method","linear method choice",
2935 U_p_string(val,factor_names[4]),
2936 U_p_strings(lo,factor_names),
2937 U_p_int(hi,sizeof(factor_names)/sizeof(char *)),1);
2938 SLV_CPARM_MACRO(FACTOR_OPTION_PTR,parameters);
2939
2940 return 1;
2941 }
2942
2943 /*
2944 * External routines
2945 * -----------------
2946 * See slv_client.h
2947 */
2948
2949 static SlvClientToken slv5_create(slv_system_t server, int *statusindex)
2950 {
2951 slv5_system_t sys;
2952
2953 sys = (slv5_system_t)asccalloc(1, sizeof(struct slv5_system_structure) );
2954 if (sys==NULL) {
2955 *statusindex = 1;
2956 return sys;
2957 }
2958 SERVER = server;
2959 sys->p.parms = sys->pa;
2960 sys->p.dynamic_parms = 0;
2961 slv5_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
2962 sys->integrity = OK;
2963 sys->presolved = 0;
2964 sys->p.output.more_important = stdout;
2965 sys->p.output.less_important = stdout;
2966 sys->J.old_partition = TRUE;
2967 sys->p.whose = (*statusindex);
2968 sys->s.ok = TRUE;
2969 sys->s.calc_ok = TRUE;
2970 sys->s.costsize = 0;
2971 sys->s.cost = NULL; /*redundant, but sanity preserving */
2972 sys->vlist = slv_get_solvers_var_list(server);
2973 sys->rlist = slv_get_solvers_rel_list(server);
2974 sys->obj = slv_get_obj_relation(server);
2975 if (sys->vlist == NULL) {
2976 ascfree(sys);
2977 FPRINTF(stderr,"IPSlv called with no variables.\n");
2978 *statusindex = -2;
2979 return NULL; /*_prolly leak here */
2980 }
2981 if (sys->rlist == NULL && sys->obj == NULL) {
2982 ascfree(sys);
2983 FPRINTF(stderr,"IPSlv called with no relations or objective.\n");
2984 *statusindex = -1;
2985 return NULL; /*_prolly leak here */
2986 }
2987 /* we don't give a damn about the objective list or the pars or
2988 * bounds or extrels or any of the other crap.
2989 */
2990 slv_check_var_initialization(server);
2991 *statusindex = 0;
2992 return((SlvClientToken)sys);
2993
2994 }
2995
2996 static void destroy_matrices( slv5_system_t sys)
2997 {
2998 if( sys->J.sys ) {
2999 int count = linsolqr_number_of_rhs(sys->J.sys)-1;
3000 for( ; count >= 0; count-- ) {
3001 destroy_array(linsolqr_get_rhs(sys->J.sys,count));
3002 }
3003 mtx_destroy(linsolqr_get_matrix(sys->J.sys));
3004 linsolqr_set_matrix(sys->J.sys,NULL);
3005 linsolqr_destroy(sys->J.sys);
3006 sys->J.