/[ascend]/trunk/solvers/ipslv/slv5.scaled.c
ViewVC logotype

Contents of /trunk/solvers/ipslv/slv5.scaled.c

Parent Directory Parent Directory | Revision Log Revision Log


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