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

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