/[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 1039 - (show annotations) (download) (as text)
Thu Jan 4 23:21:20 2007 UTC (12 years, 11 months ago) by johnpye
File MIME type: text/x-csrc
File size: 138290 byte(s)
Fixed up some #includes in compiler .[ch] files.
Switched instantiate.c to using 'asc_assert' instead of 'assert'.
Added some missing GPL headers in C++ code.
Silenced some slv3.c debug output.
Switch void-return to int-return in slv9_presolve etc (slv9.c)
Attemping to fix solvernotes.py for the commandline environment (browser==None)
Removed redundant solve(SELF) in thermalequilibrium2.a4c.
Some error reporting from addone_calc (extfntest.c).
Expanded test size in extrelfor.a4c.
Big rearrangement of bboxtest.c for top-down style.
Fixed TestFreesteam.testintegrator, added end-value checks.

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