/[ascend]/trunk/base/generic/solver/slv3.c
ViewVC logotype

Contents of /trunk/base/generic/solver/slv3.c

Parent Directory Parent Directory | Revision Log Revision Log


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