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

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

Parent Directory Parent Directory | Revision Log Revision Log


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