/[ascend]/trunk/ascend4/solver/slv3.c
ViewVC logotype

Contents of /trunk/ascend4/solver/slv3.c

Parent Directory Parent Directory | Revision Log Revision Log


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