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

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

Parent Directory Parent Directory | Revision Log Revision Log


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