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

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

Parent Directory Parent Directory | Revision Log Revision Log


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