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