/[ascend]/trunk/ipslv-temp/slv4.comp.c
ViewVC logotype

Contents of /trunk/ipslv-temp/slv4.comp.c

Parent Directory Parent Directory | Revision Log Revision Log


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