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

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

Parent Directory Parent Directory | Revision Log Revision Log


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