/[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 1034 - (show annotations) (download) (as text)
Thu Jan 4 05:37:55 2007 UTC (13 years, 4 months ago) by johnpye
File MIME type: text/x-csrc
File size: 137519 byte(s)
Switch 'void slv_*' function to output a 0-on-success error flag.
Added some debug output to ensure that such output gets through.
Fixed TestBlackBox.testfail1 to work, using above changes.
Simulation::solve now throws an exception on failure (will need to modify PyGTK GUI accordingly)
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 boolean 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;
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 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
553 for (i = 0; i < len; i++) {
554 if (rel_apply_filter(rlist[i],&rfilter)) {
555 relman_eval(rlist[i],&calc_ok,SAFE_CALC);
556 #if DEBUG
557 if (calc_ok == FALSE) {
558 FPRINTF(stderr,"error in calc_objectives\n");
559 calc_ok = TRUE;
560 }
561 #endif
562 }
563 }
564 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
565 return calc_ok;
566 }
567
568
569 /**
570 Calculates all of the residuals of included inequalities.
571 Returns true iff (calculations preceded without error and
572 all inequalities are satisfied.)
573 */
574 static boolean calc_inequalities( slv3_system_t sys){
575 struct rel_relation **rp;
576 boolean satisfied=TRUE;
577 static rel_filter_t rfilter;
578 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
579 rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
580
581 boolean calc_ok = TRUE;
582 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
583 for (rp=sys->rlist;*rp != NULL; rp++) {
584 if (rel_apply_filter(*rp,&rfilter)) {
585 relman_eval(*rp,&calc_ok,SAFE_CALC);
586 satisfied= satisfied &&
587 relman_calc_satisfied(*rp,FEAS_TOL);
588 }
589 }
590 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
591 return (calc_ok && satisfied);
592 }
593
594 /**
595 Calculates all of the residuals in the current block and computes
596 the residual norm for block status.
597
598 @return 0 on failure, non-zero on success
599 */
600 static boolean calc_residuals( slv3_system_t sys){
601 int32 row;
602 struct rel_relation *rel;
603 double time0;
604
605 if( sys->residuals.accurate ) return TRUE;
606
607 boolean calc_ok = TRUE;
608 boolean calc_ok_1;
609 row = sys->residuals.rng->low;
610 time0=tm_cpu_time();
611 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
612 for( ; row <= sys->residuals.rng->high; row++ ) {
613 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
614 #if DEBUG
615 if (!rel) {
616 int r;
617 r=mtx_row_to_org(sys->J.mtx,row);
618 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR
619 ,"QRSlv::calc_residuals: NULL relation found at ropw %d rel %d !"
620 ,(int)row,r
621 );
622 }
623 #endif
624 sys->residuals.vec[row] = relman_eval(rel,&calc_ok_1,SAFE_CALC);
625 if(!calc_ok_1){
626 calc_ok = FALSE;
627 CONSOLE_DEBUG("error calculating residual for row %d",row);
628 }
629
630 if (strcmp(CONVOPT,"ABSOLUTE") == 0) {
631 relman_calc_satisfied(rel,FEAS_TOL);
632 } else if (strcmp(CONVOPT,"RELNOM_SCALE") == 0) {
633 relman_calc_satisfied_scaled(rel,FEAS_TOL);
634 }
635 }
636 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
637 sys->s.block.functime += (tm_cpu_time() -time0);
638 sys->s.block.funcs++;
639 square_norm( &(sys->residuals) );
640 sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
641 if(!calc_ok)CONSOLE_DEBUG("error calculating residuals");
642 return(calc_ok);
643 }
644
645
646 /**
647 Calculates the current block of the jacobian.
648 It is initially unscaled.
649 */
650 static boolean calc_J( slv3_system_t sys){
651 int32 row;
652 var_filter_t vfilter;
653 double time0;
654 real64 resid;
655
656 if( sys->J.accurate )
657 return TRUE;
658
659 calc_ok = TRUE;
660 vfilter.matchbits = (VAR_INBLOCK | VAR_ACTIVE);
661 vfilter.matchvalue = (VAR_INBLOCK | VAR_ACTIVE);
662 time0=tm_cpu_time();
663 mtx_clear_region(sys->J.mtx,&(sys->J.reg));
664 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
665 struct rel_relation *rel;
666 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
667 relman_diffs(rel,&vfilter,sys->J.mtx,&resid,SAFE_CALC);
668 }
669 sys->s.block.jactime += (tm_cpu_time() - time0);
670 sys->s.block.jacs++;
671
672 if( --(sys->update.nominals) <= 0 ) sys->nominals.accurate = FALSE;
673 if( --(sys->update.weights) <= 0 ) sys->weights.accurate = FALSE;
674
675 linsolqr_matrix_was_changed(sys->J.sys);
676 return(calc_ok);
677 }
678
679
680 /**
681 Retrieves the nominal values of all of the block variables,
682 insuring that they are all strictly positive.
683 */
684 static void calc_nominals( slv3_system_t sys){
685 int32 col;
686 FILE *fp = MIF(sys);
687
688 if( sys->nominals.accurate ) return;
689 fp = MIF(sys);
690 col = sys->nominals.rng->low;
691 if(strcmp(SCALEOPT,"NONE") == 0 ||
692 strcmp(SCALEOPT,"ITERATIVE") == 0){
693 for( ; col <= sys->nominals.rng->high; col++ ) {
694 sys->nominals.vec[col] = 1;
695 }
696 } else {
697 for( ; col <= sys->nominals.rng->high; col++ ) {
698 struct var_variable *var;
699 real64 n;
700 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
701 n = var_nominal(var);
702 if( n <= 0.0 ) {
703 if( n == 0.0 ) {
704 n = TOO_SMALL;
705
706 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);
707 FPRINTF(fp,"QRSlv::calc_nominals: Variable '");
708 print_var_name(fp,sys,var);
709 FPRINTF(fp,"' has nominal value of zero. Resetting to %g.",n);
710 error_reporter_end_flush();
711
712 var_set_nominal(var,n);
713 } else {
714 n = -n;
715
716 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);
717 FPRINTF(fp,"QRSlv::calc_nominals Variable ");
718 print_var_name(fp,sys,var);
719 FPRINTF(fp,"has negative nominal value. Resetting to %g.",n);
720 error_reporter_end_flush();
721
722 var_set_nominal(var,n);
723 }
724 }
725 #if DEBUG
726 FPRINTF(fp,"Column %d is");
727 print_var_name(fp,sys,var);
728 FPRINTF(fp,"\nScaling of column %d is %g\n",col,n);
729 #endif
730 sys->nominals.vec[col] = n;
731 }
732 }
733 square_norm( &(sys->nominals) );
734 sys->update.nominals = UPDATE_NOMINALS;
735 sys->nominals.accurate = TRUE;
736 }
737
738 /**
739 Calculates the weights of all of the block relations
740 to scale the rows of the Jacobian.
741 */
742 static void calc_weights( slv3_system_t sys){
743 mtx_coord_t nz;
744 real64 sum;
745
746 if( sys->weights.accurate )
747 return;
748
749 nz.row = sys->weights.rng->low;
750 if(strcmp(SCALEOPT,"NONE") == 0 ||
751 strcmp(SCALEOPT,"ITERATIVE") == 0) {
752 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
753 sys->weights.vec[nz.row] = 1;
754 }
755 } else if (strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
756 strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0) {
757 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
758 sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
759 sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
760 }
761 } else if (strcmp(SCALEOPT,"RELNOM") == 0 ||
762 strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0) {
763 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
764 sys->weights.vec[nz.row] =
765 1.0/rel_nominal(sys->rlist[mtx_row_to_org(sys->J.mtx,nz.row)]);
766 }
767 }
768 square_norm( &(sys->weights) );
769 sys->update.weights = UPDATE_WEIGHTS;
770 sys->residuals.accurate = FALSE;
771 sys->weights.accurate = TRUE;
772 }
773
774 /**
775 Scales the jacobian.
776 */
777 static void scale_J( slv3_system_t sys){
778 int32 row;
779 int32 col;
780
781 if( sys->J.accurate ) return;
782
783 calc_nominals(sys);
784 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ )
785 mtx_mult_col(sys->J.mtx,col,sys->nominals.vec[col],&(sys->J.reg.row));
786
787 calc_weights(sys);
788 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ )
789 mtx_mult_row(sys->J.mtx,row,sys->weights.vec[row],&(sys->J.reg.col));
790 }
791
792 /**
793 ...?
794 */
795 static void jacobian_scaled(slv3_system_t sys){
796 int32 col;
797 if (DUMPCNORM) {
798 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
799 real64 cnorm;
800 cnorm =
801 calc_sqrt_D0(mtx_sum_sqrs_in_col(sys->J.mtx,col,&(sys->J.reg.row)));
802 if (cnorm >CNHIGH || cnorm <CNLOW) {
803 FPRINTF(stderr,"[col %d org %d] %g\n", col,
804 mtx_col_to_org(sys->J.mtx,col), cnorm);
805 }
806 }
807 }
808
809 sys->update.jacobian = UPDATE_JACOBIAN;
810 sys->J.accurate = TRUE;
811 sys->J.singular = FALSE; /* yet to be determined */
812 #if DEBUG
813 ERROR_REPORTER_START_HERE(ASC_PROG_NOTE);
814 FPRINTF(ASCERR,"Jacobian:\n");
815 debug_out_jacobian(stderr,sys);
816 error_reporter_end_flush();
817 #endif
818 }
819
820 /**
821 ...?
822 */
823 static void scale_variables( slv3_system_t sys){
824 int32 col;
825
826 if( sys->variables.accurate ) return;
827
828 col = sys->variables.rng->low;
829 for( ; col <= sys->variables.rng->high; col++ ) {
830 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
831 sys->variables.vec[col] = var_value(var)/sys->nominals.vec[col];
832 }
833 square_norm( &(sys->variables) );
834 sys->variables.accurate = TRUE;
835 #if DEBUG
836 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Variables: ");
837 debug_out_vector(LIF(sys),sys,&(sys->variables));
838 #endif
839 }
840
841 /**
842 Scales the previously calculated residuals.
843 */
844 static void scale_residuals( slv3_system_t sys){
845 int32 row;
846
847 if( sys->residuals.accurate ) return;
848
849 row = sys->residuals.rng->low;
850 for( ; row <= sys->residuals.rng->high; row++ ) {
851 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
852 sys->residuals.vec[row] = rel_residual(rel)*sys->weights.vec[row];
853 }
854 square_norm( &(sys->residuals) );
855 sys->residuals.accurate = TRUE;
856 #if DEBUG
857 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Residuals: ");
858 debug_out_vector(LIF(sys),sys,&(sys->residuals));
859 #endif
860 }
861
862 /**
863 Calculates relnoms for all relations in sys
864 using variable nominals.
865 */
866 static void calc_relnoms(slv3_system_t sys){
867 int32 row, col;
868 struct var_variable *var;
869 struct rel_relation *rel;
870 real64 *var_list;
871
872 /* CONSOLE_DEBUG("Begin 'calc_relnoms'"); */
873
874 var_list = ASC_NEW_ARRAY_OR_NULL(real64,sys->cap);
875 col = 0;
876 var = sys->vlist[col];
877 /* store current variable values and
878 set variable value to nominal value */
879 while(var != NULL){
880 var_list[col] = var_value(var);
881 var_set_value(var, var_nominal(var));
882 col++;
883 var = sys->vlist[col];
884 }
885 row = 0;
886 rel = sys->rlist[row];
887 /* calculate relation nominal */
888 while(rel != NULL){
889 relman_scale(rel);
890 row++;
891 rel = sys->rlist[row];
892 }
893 col = 0;
894 var = sys->vlist[col];
895 /* restore variable values */
896 while(var != NULL){
897 var_set_value(var, var_list[col]);
898 col++;
899 var = sys->vlist[col];
900 }
901 destroy_array(var_list);
902
903 /* CONSOLE_DEBUG("End 'calc_relnoms'"); */
904 }
905
906
907 /**
908 Returns the maximum ratio of magnitudes of any two nonzero
909 elements in the same column of mtx. Only considers elements
910 in region reg.
911 */
912 static real64 col_max_ratio(mtx_matrix_t *mtx,
913 mtx_region_t *reg
914 ){
915 real64 ratio;
916 real64 max_ratio;
917 real64 num, denom, dummy;
918 mtx_coord_t coord;
919
920 max_ratio = 0;
921 for(coord.col = reg->col.low;
922 coord.col <= reg->col.high; coord.col++) {
923 ratio = 0;
924 num = mtx_col_max(*mtx,&(coord),&(reg->row),&(dummy));
925 denom = mtx_col_min(*mtx,&(coord),&(reg->row),&(dummy),1e-7);
926 if(denom >0){
927 ratio = num/denom;
928 }
929 if(ratio > 10000000){
930 /* FPRINTF(stderr,"HELPME\n");*/
931 }
932 if(ratio > max_ratio){
933 max_ratio = ratio;
934 }
935 }
936 if(max_ratio == 0){
937 max_ratio = 1;
938 }
939 return max_ratio;
940 }
941
942 /**
943 Returns the maximum ratio of magnitudes of any two nonzero
944 elements in the same row of mtx. Only considers elements
945 in region reg.
946 */
947 static real64 row_max_ratio(mtx_matrix_t *mtx,
948 mtx_region_t *reg
949 ){
950 real64 ratio;
951 real64 max_ratio;
952 real64 num, denom, dummy;
953 mtx_coord_t coord;
954 max_ratio = 0;
955
956 for(coord.row = reg->row.low;
957 coord.row <= reg->row.high; coord.row++) {
958 ratio = 0;
959 num = mtx_row_max(*mtx,&(coord),&(reg->col),&(dummy));
960 denom = mtx_row_min(*mtx,&(coord),&(reg->col),&(dummy),1e-7);
961 if(denom >0){
962 ratio = num/denom;
963 }
964 if(ratio > 10000000){
965 /* FPRINTF(stderr,"HELPME\n");*/
966 }
967 if(ratio > max_ratio){
968 max_ratio = ratio;
969 }
970 }
971 if(max_ratio == 0){
972 max_ratio = 1;
973 }
974 return max_ratio;
975 }
976
977 /**
978 Calculates scaling factor suggested by Fourer.
979 For option = 0, returns scaling factor for
980 row number loc.
981 For option = 1, returns scaling factor for
982 col number loc.
983 */
984 static real64 calc_fourer_scale(mtx_matrix_t mtx,
985 mtx_region_t reg,
986 int32 loc,
987 int32 option
988 ){
989 mtx_coord_t coord;
990 real64 min, max, dummy;
991 real64 scale;
992
993 if(option == 0){
994 if((loc < reg.row.low) || (loc > reg.row.high)){
995 return 1;
996 }
997 coord.row = loc;
998 min = mtx_row_min(mtx,&(coord),&(reg.col),&(dummy),1e-7);
999 max = mtx_row_max(mtx,&(coord),&(reg.col),&(dummy));
1000 scale = min*max;
1001 if(scale > 0){
1002 scale = sqrt(scale);
1003 } else {
1004 scale = 1;
1005 }
1006 return scale;
1007 } else {
1008 if(loc < reg.col.low || loc > reg.col.high){
1009 return 1;
1010 }
1011 coord.col = loc;
1012 min = mtx_col_min(mtx,&(coord),&(reg.row),&(dummy),1e-7);
1013 max = mtx_col_max(mtx,&(coord),&(reg.row),&(dummy));
1014 scale = min*max;
1015 if(scale > 0){
1016 scale = sqrt(scale);
1017 } else {
1018 scale = 1;
1019 }
1020 return scale;
1021 }
1022 }
1023
1024 /**
1025 This funcion is an implementation of the scaling
1026 routine by Fourer on p304 of Mathematical Programing
1027 vol 23, (1982).
1028 This function will scale the Jacobian and store the scaling
1029 factors in sys->nominals and sys->weights.
1030 If the Jacobian has been previously scaled
1031 by another method (during this iteration) then these vectors
1032 should contain the scale factors used in that scaling.
1033 */
1034 static void scale_J_iterative(slv3_system_t sys){
1035 real64 rho_col_old, rho_col_new;
1036 real64 rho_row_old, rho_row_new;
1037 int32 k;
1038 int32 done;
1039 int32 row, col;
1040
1041 real64 *colvec = sys->nominals.vec;
1042 real64 *rowvec = sys->weights.vec;
1043 real64 rowscale, colscale;
1044
1045 rho_col_old = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1046 rho_row_old = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1047 k = 0;
1048 done = 0;
1049 while(done == 0){
1050 k++;
1051 for(row = sys->J.reg.row.low;
1052 row <= sys->J.reg.row.high; row++){
1053 rowscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,row,0);
1054 mtx_mult_row(sys->J.mtx,row,rowscale,&(sys->J.reg.col));
1055 rowvec[row] *= rowscale;
1056 }
1057 for(col = sys->J.reg.col.low;
1058 col <= sys->J.reg.col.high; col++){
1059 colscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,col,1);
1060 mtx_mult_col(sys->J.mtx,col,colscale,&(sys->J.reg.row));
1061 colvec[col] *= colscale;
1062 }
1063 rho_col_new = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1064 rho_row_new = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1065 if((rho_col_new >= ITSCALETOL*rho_col_old &&
1066 rho_row_new >= ITSCALETOL*rho_row_old)
1067 || k >= ITSCALELIM){
1068 done = 1;
1069 /* FPRINTF(stderr,"%d ITERATIVE SCALING ITERATIONS\n",k);*/
1070 } else {
1071 rho_row_old = rho_row_new;
1072 rho_col_old = rho_col_new;
1073 }
1074 }
1075 square_norm( &(sys->nominals) );
1076 sys->update.nominals = UPDATE_NOMINALS;
1077 sys->nominals.accurate = TRUE;
1078
1079 square_norm( &(sys->weights) );
1080 sys->update.weights = UPDATE_WEIGHTS;
1081 sys->residuals.accurate = FALSE;
1082 sys->weights.accurate = TRUE;
1083 }
1084
1085 /**
1086 Scale system dependent on interface parameters
1087 */
1088 static void scale_system( slv3_system_t sys ){
1089 if(strcmp(SCALEOPT,"NONE") == 0){
1090 if(sys->J.accurate == FALSE){
1091 calc_nominals(sys);
1092 calc_weights(sys);
1093 jacobian_scaled(sys);
1094 }
1095 scale_variables(sys);
1096 scale_residuals(sys);
1097 return;
1098 }
1099 if(strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
1100 strcmp(SCALEOPT,"RELNOM") == 0){
1101 if(sys->J.accurate == FALSE){
1102 scale_J(sys);
1103 jacobian_scaled(sys);
1104 }
1105 scale_variables(sys);
1106 scale_residuals(sys);
1107 return;
1108 }
1109 if(strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0 ||
1110 strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0){
1111 if(sys->J.accurate == FALSE){
1112 --sys->update.iterative;
1113 if(sys->update.iterative <= 0) {
1114 scale_J(sys);
1115 scale_J_iterative(sys);
1116 sys->update.iterative =
1117 UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1118 } else {
1119 sys->weights.accurate = TRUE;
1120 sys->nominals.accurate = TRUE;
1121 scale_J(sys); /* will use current scaling vectors */
1122 }
1123 jacobian_scaled(sys);
1124 }
1125 scale_variables(sys);
1126 scale_residuals(sys);
1127 return;
1128 }
1129 if(strcmp(SCALEOPT,"ITERATIVE") == 0){
1130 if(sys->J.accurate == FALSE){
1131 --sys->update.iterative;
1132 if(sys->update.iterative <= 0) {
1133 calc_nominals(sys);
1134 calc_weights(sys);
1135 scale_J_iterative(sys);
1136 sys->update.iterative =
1137 UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1138 } else {
1139 sys->weights.accurate = TRUE;
1140 sys->nominals.accurate = TRUE;
1141 scale_J(sys); /* will use current scaling vectors */
1142 }
1143 jacobian_scaled(sys);
1144 }
1145 scale_variables(sys);
1146 scale_residuals(sys);
1147 }
1148 return;
1149 }
1150
1151 /**
1152 Calculate scaled gradient of the objective function.
1153
1154 @TODO This entire function needs to be reimplemented with relman_diffs.
1155 */
1156 static boolean calc_gradient(slv3_system_t sys){
1157
1158 if( sys->gradient.accurate ) return TRUE;
1159
1160 calc_ok = TRUE;
1161 if ( !OPTIMIZING(sys) ) {
1162 zero_vector(&(sys->gradient));
1163 sys->gradient.norm2 = 0.0;
1164 } else {
1165 Asc_Panic(2, "calc_gradient", "Not implemented");
1166 #if CANOPTIMIZE
1167 real64 pd;
1168 const struct var_variable **vp;
1169 var_filter_t vfilter;
1170 vfilter.matchbits = (VAR_INBLOCK | VAR_SVAR | VAR_ACTIVE);
1171 vfilter.matchvalue = (VAR_INBLOCK | VAR_SVAR | VAR_ACTIVE);
1172 zero_vector(&(sys->gradient));
1173 /* the next line will core dump anyway since vp not null-terminated*/
1174 for( vp = rel_incidence_list(sys->obj) ; *vp != NULL ; ++vp ) {
1175 int32 col;
1176 col = mtx_org_to_col(sys->J.mtx,var_sindex(*vp));
1177 if( var_apply_filter(*vp,&vfilter) ) {
1178 /* the next line will core dump anyway since _diff not implemented */
1179 relman_diff(sys->obj,*vp,&pd,SAFE_CALC); /* barf */
1180 sys->gradient.vec[col] = sys->nominals.vec[col]*pd;
1181 }
1182 }
1183 #endif
1184 square_norm( &(sys->gradient) );
1185 }
1186 sys->gradient.accurate = TRUE;
1187 #if DEBUG
1188 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Gradient: ");
1189 debug_out_vector(LIF(sys),sys,&(sys->gradient));
1190 #endif
1191 return calc_ok;
1192 }
1193
1194 /**
1195 Create a new hessian_data structure for storing
1196 latest update information.
1197 */
1198 static void create_update(slv3_system_t sys){
1199 struct hessian_data *update;
1200
1201 if( !OPTIMIZING(sys) )
1202 return;
1203
1204 update = ASC_NEW(struct hessian_data);
1205 update->y.vec = ASC_NEW_ARRAY_OR_NULL(real64,sys->cap);
1206 update->y.rng = &(sys->J.reg.col);
1207 update->y.accurate = FALSE;
1208 update->Bs.vec = ASC_NEW_ARRAY_OR_NULL(real64,sys->cap);
1209 update->Bs.rng = &(sys->J.reg.col);
1210 update->Bs.accurate = FALSE;
1211 update->next = sys->B;
1212 sys->B = update;
1213 }
1214
1215
1216 /**
1217 Computes a rank 2 BFGS update to the hessian matrix
1218 B which accumulates curvature.
1219 */
1220 static void calc_B( slv3_system_t sys){
1221 if( sys->s.block.iteration > 1 ) {
1222 create_update(sys);
1223 } else {
1224 if( sys->B ) {
1225 struct hessian_data *update;
1226 for( update=sys->B; update != NULL; ) {
1227 struct hessian_data *handle;
1228 handle = update;
1229 update = update->next;
1230 destroy_array(handle->y.vec);
1231 destroy_array(handle->Bs.vec);
1232 ascfree(handle);
1233 }
1234 sys->B = NULL;
1235 }
1236 }
1237 if( sys->B ) {
1238 real64 theta;
1239 /*
1240 * The y vector
1241 */
1242 if( !sys->B->y.accurate ) {
1243 int32 col;
1244 matrix_product(sys->J.mtx, &(sys->multipliers),
1245 &(sys->B->y), 1.0, TRUE);
1246 col = sys->B->y.rng->low;
1247 for( ; col <= sys->B->y.rng->high; col++ ) {
1248 sys->B->y.vec[col] += sys->gradient.vec[col] -
1249 sys->stationary.vec[col];
1250 }
1251 square_norm( &(sys->B->y) );
1252 sys->B->y.accurate = TRUE;
1253 }
1254
1255 /*
1256 * The Bs vector
1257 */
1258 if( !sys->B->Bs.accurate ) {
1259 struct hessian_data *update;
1260 copy_vector(&(sys->varstep),&(sys->B->Bs));
1261 for( update=sys->B->next; update != NULL; update = update->next ) {
1262 int32 col;
1263 real64 ys = inner_product( &(update->y),&(sys->varstep) );
1264 real64 sBs = inner_product( &(update->Bs),&(sys->varstep) );
1265 col = sys->B->Bs.rng->low;
1266 for( ; col<=sys->B->Bs.rng->high; col++) {
1267 sys->B->Bs.vec[col] += update->ys > 0.0 ?
1268 (update->y.vec[col])*ys/update->ys : 0.0;
1269 sys->B->Bs.vec[col] -= update->sBs > 0.0 ?
1270 (update->Bs.vec[col])*sBs/update->sBs : 0.0;
1271 }
1272 }
1273 square_norm( &(sys->B->Bs) );
1274 sys->B->Bs.accurate = TRUE;
1275 }
1276
1277 sys->B->ys = inner_product( &(sys->B->y),&(sys->varstep) );
1278 sys->B->sBs = inner_product( &(sys->B->Bs),&(sys->varstep) );
1279
1280 if( sys->B->ys == 0.0 && sys->B->sBs == 0.0 ) {
1281 theta = 0.0;
1282 } else {
1283 theta = sys->B->ys < POSITIVE_DEFINITE*sys->B->sBs ?
1284 (1.0-POSITIVE_DEFINITE)*sys->B->sBs/(sys->B->sBs - sys->B->ys):1.0;
1285 }
1286 #if DEBUG
1287 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"ys, sBs, PD, theta = %g, %g, %g, %g\n",
1288 sys->B->ys,
1289 sys->B->sBs,
1290 POSITIVE_DEFINITE,
1291 theta);
1292 #endif
1293 if( theta < 1.0 ) {
1294 int32 col;
1295 col = sys->B->y.rng->low;
1296 for( ; col <= sys->B->y.rng->high; col++ )
1297 sys->B->y.vec[col] = theta*sys->B->y.vec[col] +
1298 (1.0-theta)*sys->B->Bs.vec[col];
1299 square_norm( &(sys->B->y) );
1300 sys->B->ys = theta*sys->B->ys + (1.0-theta)*sys->B->sBs;
1301 }
1302 }
1303 }
1304
1305
1306 /**
1307 Obtain the equations and variables which
1308 are able to be pivoted.
1309 return value is the row rank deficiency, which we hope is 0.
1310 */
1311 static int calc_pivots(slv3_system_t sys){
1312 int row_rank_defect=0, oldtiming;
1313 linsolqr_system_t lsys = sys->J.sys;
1314 FILE *fp = LIF(sys);
1315
1316 oldtiming = g_linsolqr_timing;
1317 g_linsolqr_timing =LINTIME;
1318 linsolqr_factor(lsys,sys->J.fm); /* factor */
1319 g_linsolqr_timing = oldtiming;
1320
1321 if (OPTIMIZING(sys)) {
1322 /* need things for nullspace move. don't care about
1323 * dependency coefficiency in any circumstances at present.
1324 */
1325 linsolqr_calc_col_dependencies(lsys);
1326 set_null(sys->J.relpivots,sys->cap);
1327 set_null(sys->J.varpivots,sys->cap);
1328 linsolqr_get_pivot_sets(lsys,sys->J.relpivots,sys->J.varpivots);
1329 }
1330
1331 sys->J.rank = linsolqr_rank(lsys);
1332 sys->J.singular = FALSE;
1333 row_rank_defect = sys->J.reg.row.high - sys->J.reg.row.low+1 - sys->J.rank;
1334 if( row_rank_defect > 0 ) {
1335 int32 row,krow;
1336 mtx_sparse_t *uprows=NULL;
1337 sys->J.singular = TRUE;
1338 uprows = linsolqr_unpivoted_rows(lsys);
1339 if (uprows !=NULL) {
1340 for( krow=0; krow < uprows->len ; krow++ ) {
1341 int32 org_row;
1342 struct rel_relation *rel;
1343
1344 org_row = uprows->idata[krow];
1345 row = mtx_org_to_row(sys->J.mtx,org_row);
1346 rel = sys->rlist[org_row];
1347
1348 ERROR_REPORTER_START_HERE(ASC_PROG_ERROR);
1349 FPRINTF(stderr,"Relation '");
1350 print_rel_name(stderr,sys,rel);
1351 FPRINTF(stderr,"' not pivoted.");
1352 error_reporter_end_flush();
1353
1354 /*
1355 * assign zeros to the corresponding weights
1356 * so that subsequent calls to "scale_residuals"
1357 * will only measure the pivoted equations.
1358 */
1359 sys->weights.vec[row] = 0.0;
1360 sys->residuals.vec[row] = 0.0;
1361 sys->residuals.accurate = FALSE;
1362 mtx_mult_row(sys->J.mtx,row,0.0,&(sys->J.reg.col));
1363 }
1364 mtx_destroy_sparse(uprows);
1365 }
1366 if( !sys->residuals.accurate ) {
1367 square_norm( &(sys->residuals) );
1368 sys->residuals.accurate = TRUE;
1369 sys->update.weights = 0; /* re-compute weights next iteration. */
1370 }
1371 }
1372 if( sys->J.rank < sys->J.reg.col.high-sys->J.reg.col.low+1 ) {
1373 int32 col,kcol;
1374 mtx_sparse_t *upcols=NULL;
1375 if (NOTNULL(upcols)) {
1376 for( kcol=0; upcols != NULL && kcol < upcols->len ; kcol++ ) {
1377 int32 org_col;
1378 struct var_variable *var;
1379
1380 org_col = upcols->idata[kcol];
1381 col = mtx_org_to_col(sys->J.mtx,org_col);
1382 var = sys->vlist[org_col];
1383 FPRINTF(fp,"%-40s ---> ","Variable not pivoted");
1384 print_var_name(fp,sys,var);
1385 PUTC('\n',fp);
1386 /*
1387 * If we're not optimizing (everything should be
1388 * pivotable) or this was one of the dependent variables,
1389 * consider this variable as if it were fixed.
1390 */
1391 if( col <= sys->J.reg.col.high - sys->ZBZ.order ) {
1392 mtx_mult_col(sys->J.mtx,col,0.0,&(sys->J.reg.row));
1393 }
1394 }
1395 mtx_destroy_sparse(upcols);
1396 }
1397 }
1398 if (SHOW_LESS_IMPT) {
1399 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d (%s)\n","Jacobian rank", sys->J.rank,
1400 sys->J.singular ? "deficient":"full");
1401 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n","Smallest pivot",
1402 linsolqr_smallest_pivot(sys->J.sys));
1403 }
1404 return row_rank_defect;
1405 }
1406
1407 /**
1408 Updates the reduced hessian matrix.
1409 if !OPTIMIZING just sets zbz.accurate true and returns.
1410 */
1411 static void calc_ZBZ(slv3_system_t sys){
1412 mtx_coord_t nz;
1413
1414 if( sys->ZBZ.accurate ) return;
1415
1416 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1417 for( nz.col = 0; nz.col <= nz.row; nz.col++ ) {
1418 int32 col, depr, depc;
1419 col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
1420 depr = mtx_col_to_org(sys->J.mtx,col);
1421 col = nz.col+sys->J.reg.col.high+1-sys->ZBZ.order;
1422 depc = mtx_col_to_org(sys->J.mtx,col);
1423 sys->ZBZ.mtx[nz.row][nz.col] = (nz.row==nz.col ? 1.0 : 0.0);
1424 col = sys->J.reg.col.low;
1425 for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
1426 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1427 if( set_is_member(sys->J.varpivots,ind) )
1428 sys->ZBZ.mtx[nz.row][nz.col] +=
1429 (-linsolqr_org_col_dependency(sys->J.sys,depr,ind))*
1430 (-linsolqr_org_col_dependency(sys->J.sys,depc,ind));
1431 }
1432 if( nz.row != nz.col ) {
1433 sys->ZBZ.mtx[nz.col][nz.row] =
1434 sys->ZBZ.mtx[nz.row][nz.col];
1435 }
1436 }
1437 }
1438 if( OPTIMIZING(sys) ) {
1439 struct hessian_data *update;
1440 for( update=sys->B; update != NULL; update = update->next ) {
1441 for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1442 int32 col, dep;
1443 col = nz.row + sys->J.reg.col.high + 1 - sys->ZBZ.order;
1444 dep = mtx_col_to_org(sys->J.mtx,col);
1445 sys->ZBZ.Zy[nz.row] = update->y.vec[col];
1446 sys->ZBZ.ZBs[nz.row] = update->Bs.vec[col];
1447 col = sys->J.reg.col.low;
1448 for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
1449 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1450 if( set_is_member(sys->J.varpivots,ind) ) {
1451 sys->ZBZ.Zy[nz.row] += update->y.vec[col]*
1452 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1453 sys->ZBZ.ZBs[nz.row] += update->Bs.vec[col]*
1454 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1455 }
1456 }
1457 for( nz.col=0; nz.col <= nz.row; nz.col++ ) {
1458 sys->ZBZ.mtx[nz.row][nz.col] += update->ys > 0.0 ?
1459 sys->ZBZ.Zy[nz.row]*sys->ZBZ.Zy[nz.col]/update->ys : 0.0;
1460 sys->ZBZ.mtx[nz.row][nz.col] -= update->sBs > 0.0 ?
1461 sys->ZBZ.ZBs[nz.row]*sys->ZBZ.ZBs[nz.col]/update->sBs : 0.0;
1462 if( nz.row != nz.col ) {
1463 sys->ZBZ.mtx[nz.col][nz.row] =
1464 sys->ZBZ.mtx[nz.row][nz.col];
1465 }
1466 }
1467 }
1468 }
1469 }
1470 sys->ZBZ.accurate = TRUE;
1471 #if DEBUG
1472 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nReduced Hessian: \n");
1473 debug_out_hessian(LIF(sys),sys);
1474 #endif
1475 }
1476
1477
1478 /**
1479 Calculates just the jacobian RHS. This function should be used to
1480 supplement calculation of the jacobian. The vector vec must
1481 already be calculated and scaled so as to simply be added to the
1482 rhs. Caller is responsible for initially zeroing the rhs vector.
1483 */
1484 static void calc_rhs(slv3_system_t sys, struct vector_data *vec,
1485 real64 scalar, boolean transpose
1486 ){
1487 if( transpose ) { /* vec is indexed by col */
1488 int32 col;
1489 for( col=vec->rng->low; col<=vec->rng->high; col++ ) {
1490 sys->J.rhs[mtx_col_to_org(sys->J.mtx,col)] += scalar*vec->vec[col];
1491 }
1492 } else { /* vec is indexed by row */
1493 int32 row;
1494 for( row=vec->rng->low; row<=vec->rng->high; row++ ) {
1495 sys->J.rhs[mtx_row_to_org(sys->J.mtx,row)] += scalar*vec->vec[row];
1496 }
1497 }
1498 linsolqr_rhs_was_changed(sys->J.sys,sys->J.rhs);
1499 }
1500
1501
1502 /**
1503 Computes the lagrange multipliers for the equality constraints.
1504 */
1505 static void calc_multipliers(slv3_system_t sys){
1506 if( sys->multipliers.accurate )
1507 return;
1508
1509 if ( !OPTIMIZING(sys) ) {
1510 zero_vector(&(sys->multipliers));
1511 sys->multipliers.norm2 = 0.0;
1512 } else {
1513 linsolqr_system_t lsys = sys->J.sys;
1514 int32 row;
1515 sys->J.rhs = linsolqr_get_rhs(lsys,0);
1516 mtx_zero_real64(sys->J.rhs,sys->cap);
1517 calc_rhs(sys, &(sys->gradient), -1.0, TRUE );
1518 linsolqr_solve(lsys,sys->J.rhs);
1519 row = sys->multipliers.rng->low;
1520 for( ; row <= sys->multipliers.rng->high; row++ ) {
1521 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1522 sys->multipliers.vec[row] = linsolqr_var_value
1523 (lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row));
1524 rel_set_multiplier(rel,sys->multipliers.vec[row]*
1525 sys->weights.vec[row]);
1526
1527 }
1528 if (SAVLIN) {
1529 FILE *ldat;
1530 int32 ov;
1531 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1532 ldat=fopen(savlinfilename,"w");
1533 FPRINTF(ldat,
1534 "================= multipliersrhs (orgcoled) itn %d =====\n",
1535 sys->s.iteration);
1536 debug_write_array(ldat,sys->J.rhs,sys->cap);
1537 FPRINTF(ldat,
1538 "================= multipliers (orgrowed) ============\n");
1539 for(ov=0 ; ov < sys->cap; ov++ )
1540 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1541 fclose(ldat);
1542 }
1543 square_norm( &(sys->multipliers) );
1544 }
1545 sys->multipliers.accurate = TRUE;
1546 #if DEBUG
1547 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Multipliers: ");
1548 debug_out_vector(LIF(sys),sys,&(sys->multipliers));
1549 #endif
1550 }
1551
1552
1553 /**
1554 Computes the gradient of the lagrangian which
1555 should be zero at the optimum solution.
1556 */
1557 static void calc_stationary( slv3_system_t sys){
1558 if( sys->stationary.accurate )
1559 return;
1560
1561 if ( !OPTIMIZING(sys) ) {
1562 zero_vector(&(sys->stationary));
1563 sys->stationary.norm2 = 0.0;
1564 } else {
1565 int32 col;
1566 matrix_product(sys->J.mtx, &(sys->multipliers),
1567 &(sys->stationary), 1.0, TRUE);
1568 col = sys->stationary.rng->low;
1569 for( ; col <= sys->stationary.rng->high; col++ )
1570 sys->stationary.vec[col] += sys->gradient.vec[col];
1571 square_norm( &(sys->stationary) );
1572 }
1573 sys->stationary.accurate = TRUE;
1574 #if DEBUG
1575 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Stationary: ");
1576 debug_out_vector(LIF(sys),sys,&(sys->stationary));
1577 #endif
1578 }
1579
1580
1581 /**
1582 Calculate the gamma vector.
1583 */
1584 static void calc_gamma( slv3_system_t sys){
1585 if( sys->gamma.accurate )
1586 return;
1587
1588 matrix_product(sys->J.mtx, &(sys->residuals),
1589 &(sys->gamma), -1.0, TRUE);
1590 square_norm( &(sys->gamma) );
1591 sys->gamma.accurate = TRUE;
1592 #if DEBUG
1593 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Gamma: ");
1594 debug_out_vector(LIF(sys),sys,&(sys->gamma));
1595 #endif
1596 }
1597
1598 /**
1599 Calculate the Jgamma vector.
1600 */
1601 static void calc_Jgamma( slv3_system_t sys){
1602 if( sys->Jgamma.accurate )
1603 return;
1604
1605 matrix_product(sys->J.mtx, &(sys->gamma),
1606 &(sys->Jgamma), 1.0, FALSE);
1607 square_norm( &(sys->Jgamma) );
1608 sys->Jgamma.accurate = TRUE;
1609 #if DEBUG
1610 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Jgamma: ");
1611 debug_out_vector(LIF(sys),sys,&(sys->Jgamma));
1612 #endif
1613 }
1614
1615
1616 /**
1617 Computes a step to solve the linearized equations.
1618 */
1619 static void calc_newton( slv3_system_t sys){
1620 linsolqr_system_t lsys = sys->J.sys;
1621 int32 col;
1622
1623 if( sys->newton.accurate )
1624 return;
1625
1626 sys->J.rhs = linsolqr_get_rhs(lsys,1);
1627 mtx_zero_real64(sys->J.rhs,sys->cap);
1628 calc_rhs(sys, &(sys->residuals), -1.0, FALSE);
1629 linsolqr_solve(lsys,sys->J.rhs);
1630 col = sys->newton.rng->low;
1631 for( ; col <= sys->newton.rng->high; col++ ) {
1632 sys->newton.vec[col] =
1633 linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1634 }
1635 if (SAVLIN) {
1636 FILE *ldat;
1637 int32 ov;
1638 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1639 ldat=fopen(savlinfilename,"w");
1640 FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1641 sys->s.iteration);
1642 debug_write_array(ldat,sys->J.rhs,sys->cap);
1643 FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1644 for(ov=0 ; ov < sys->cap; ov++ )
1645 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1646 fclose(ldat);
1647 }
1648 square_norm( &(sys->newton) );
1649 sys->newton.accurate = TRUE;
1650 #if DEBUG
1651 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Newton: ");
1652 debug_out_vector(LIF(sys),sys,&(sys->newton));
1653 #endif
1654 }
1655
1656
1657 /**
1658 Computes an update to the product B and newton.
1659 */
1660 static void calc_Bnewton( slv3_system_t sys){
1661 if( sys->Bnewton.accurate )
1662 return;
1663
1664 if ( !OPTIMIZING(sys) ) {
1665 zero_vector(&(sys->Bnewton));
1666 sys->Bnewton.norm2 = 0.0;
1667 } else {
1668 struct hessian_data *update;
1669 copy_vector(&(sys->newton),&(sys->Bnewton));
1670 for( update=sys->B; update != NULL; update = update->next ) {
1671 int32 col;
1672 real64 Yn = inner_product( &(update->y),&(sys->newton) );
1673 real64 sBn = inner_product( &(update->Bs),&(sys->newton) );
1674 col = sys->Bnewton.rng->low;
1675 for( ; col <= sys->Bnewton.rng->high; col++ ) {
1676 sys->Bnewton.vec[col] += update->ys > 0.0 ?
1677 (update->y.vec[col])*Yn/update->ys : 0.0;
1678 sys->Bnewton.vec[col] -= update->sBs > 0.0 ?
1679 (update->Bs.vec[col])*sBn/update->sBs : 0.0;
1680 }
1681 }
1682 square_norm( &(sys->Bnewton) );
1683 }
1684 sys->Bnewton.accurate = TRUE;
1685 #if DEBUG
1686 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Bnewton: ");
1687 debug_out_vector(LIF(sys),sys,&(sys->Bnewton));
1688 #endif
1689 }
1690
1691
1692 /**
1693 Calculate the nullspace move if OPTIMIZING.
1694 */
1695 static void calc_nullspace( slv3_system_t sys){
1696 if( sys->nullspace.accurate )
1697 return;
1698
1699 if( !OPTIMIZING(sys) ) {
1700 zero_vector(&(sys->nullspace));
1701 sys->nullspace.norm2 = 0.0;
1702 } else {
1703 mtx_coord_t nz;
1704 zero_vector(&(sys->nullspace));
1705 for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1706 int32 col, dep, ndx;
1707 col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
1708 dep = mtx_col_to_org(sys->J.mtx,col);
1709 sys->nullspace.vec[col] = -sys->stationary.vec[col] -
1710 sys->Bnewton.vec[col];
1711 ndx = sys->J.reg.col.low;
1712 for( ; ndx <= sys->J.reg.col.high - sys->ZBZ.order; ndx++ ) {
1713 int32 ind = mtx_col_to_org(sys->J.mtx,ndx);
1714 if( set_is_member(sys->J.varpivots,ind) )
1715 sys->nullspace.vec[col] -=
1716 (sys->stationary.vec[ndx] + sys->Bnewton.vec[ndx])*
1717 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1718 }
1719 }
1720 /*
1721 * Must invert ZBZ first. It's symmetric so
1722 * can find Cholesky factors. Essentially, find
1723 * the "square root" of the matrix such that
1724 *
1725 * T
1726 * L U = U U = ZBZ, where U is an upper triangular
1727 * matrix.
1728 */
1729 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1730 for( nz.col = nz.row; nz.col < sys->ZBZ.order; nz.col++ ) {
1731 int32 col;
1732 for( col = 0; col < nz.row; col++ )
1733 sys->ZBZ.mtx[nz.row][nz.col] -=
1734 sys->ZBZ.mtx[nz.row][col]*
1735 sys->ZBZ.mtx[col][nz.col];
1736 if( nz.row == nz.col )
1737 sys->ZBZ.mtx[nz.row][nz.col] =
1738 calc_sqrt_D0(sys->ZBZ.mtx[nz.row][nz.col]);
1739 else {
1740 sys->ZBZ.mtx[nz.row][nz.col] /=
1741 sys->ZBZ.mtx[nz.row][nz.row];
1742 sys->ZBZ.mtx[nz.col][nz.row] =
1743 sys->ZBZ.mtx[nz.row][nz.col];
1744 }
1745 }
1746 }
1747 #if DEBUG
1748 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nInverse Reduced Hessian: \n");
1749 debug_out_hessian(LIF(sys),sys);
1750 #endif
1751 /*
1752 * forward substitute
1753 */
1754 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1755 int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1756 for( nz.col = 0; nz.col < nz.row; nz.col++ ) {
1757 sys->nullspace.vec[nz.row+offset] -=
1758 sys->nullspace.vec[nz.col+offset]*
1759 sys->ZBZ.mtx[nz.row][nz.col];
1760 }
1761 sys->nullspace.vec[nz.row+offset] /=
1762 sys->ZBZ.mtx[nz.row][nz.row];
1763 }
1764
1765 /*
1766 * backward substitute
1767 */
1768 for( nz.row = sys->ZBZ.order-1; nz.row >= 0; nz.row-- ) {
1769 int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1770 for( nz.col = nz.row+1; nz.col < sys->ZBZ.order; nz.col++ ) {
1771 sys->nullspace.vec[nz.row+offset] -=
1772 sys->nullspace.vec[nz.col+offset]*
1773 sys->ZBZ.mtx[nz.row][nz.col];
1774 }
1775 sys->nullspace.vec[nz.row+offset] /=
1776 sys->ZBZ.mtx[nz.row][nz.row];
1777 }
1778 square_norm( &(sys->nullspace) );
1779 }
1780 sys->nullspace.accurate = TRUE;
1781 #if DEBUG
1782 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Nullspace: ");
1783 debug_out_vector(LIF(sys),sys,&(sys->nullspace));
1784 #endif
1785 }
1786
1787 /**
1788 Calculate the 1st order descent direction for phi
1789 in the variables.
1790 */
1791 static void calc_varstep1( slv3_system_t sys){
1792 if( sys->varstep1.accurate )
1793 return;
1794
1795 if( !OPTIMIZING(sys) ) {
1796 copy_vector(&(sys->gamma),&(sys->varstep1));
1797 sys->varstep1.norm2 = sys->gamma.norm2;
1798 } else {
1799 int32 col;
1800 col = sys->varstep1.rng->low;
1801 for( ; col <= sys->varstep1.rng->high; col++ )
1802 sys->varstep1.vec[col] = RHO*sys->gamma.vec[col] -
1803 sys->stationary.vec[col];
1804 square_norm( &(sys->varstep1) );
1805 }
1806 sys->varstep1.accurate = TRUE;
1807 #if DEBUG
1808 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Varstep1: ");
1809 debug_out_vector(LIF(sys),sys,&(sys->varstep1));
1810 #endif
1811 }
1812
1813
1814 /**
1815 Computes an update to the product B and varstep1.
1816 */
1817 static void calc_Bvarstep1( slv3_system_t sys){
1818 if( sys->Bvarstep1.accurate )
1819 return;
1820
1821 if ( !OPTIMIZING(sys) ) {
1822 zero_vector(&(sys->Bvarstep1));
1823 sys->Bvarstep1.norm2 = 0.0;
1824 } else {
1825 struct hessian_data *update;
1826 copy_vector(&(sys->varstep1),&(sys->Bvarstep1));
1827 for( update=sys->B; update != NULL; update = update->next ) {
1828 int32 col;
1829 real64 yv = inner_product( &(update->y),&(sys->varstep1) );
1830 real64 sBv = inner_product( &(update->Bs),&(sys->varstep1) );
1831 col = sys->Bvarstep1.rng->low;
1832 for( ; col <= sys->Bvarstep1.rng->high; col++ ) {
1833 sys->Bvarstep1.vec[col] += update->ys > 0.0 ?
1834 (update->y.vec[col])*yv/update->ys : 0.0;
1835 sys->Bvarstep1.vec[col] -= update->sBs > 0.0 ?
1836 (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1837 }
1838 }
1839 square_norm( &(sys->Bvarstep1) );
1840 }
1841 sys->Bvarstep1.accurate = TRUE;
1842 #if DEBUG
1843 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Bvarstep1: ");
1844 debug_out_vector(LIF(sys),sys,&(sys->Bvarstep1));
1845 #endif
1846 }
1847
1848
1849 /**
1850 Calculate the 2nd order descent direction for phi
1851 in the variables.
1852 */
1853 static void calc_varstep2( slv3_system_t sys){
1854 if( sys->varstep2.accurate )
1855 return;
1856
1857 if( !OPTIMIZING(sys) ) {
1858 copy_vector(&(sys->newton),&(sys->varstep2));
1859 sys->varstep2.norm2 = sys->newton.norm2;
1860 } else {
1861 int32 col;
1862 col = sys->varstep2.rng->low;
1863 for( ; col <= sys->varstep2.rng->high - sys->ZBZ.order ; ++col ) {
1864 int32 dep;
1865 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1866 sys->varstep2.vec[col] = sys->newton.vec[col];
1867 if( set_is_member(sys->J.varpivots,ind) ) {
1868 dep = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1869 for( ; dep <= sys->varstep2.rng->high; dep++ )
1870 sys->varstep2.vec[col] += sys->nullspace.vec[dep]*
1871 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1872 }
1873 }
1874 col = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1875 for( ; col <= sys->varstep2.rng->high; ++col )
1876 sys->varstep2.vec[col] = sys->nullspace.vec[col] +
1877 sys->newton.vec[col];
1878 square_norm( &(sys->varstep2) );
1879 }
1880 sys->varstep2.accurate = TRUE;
1881 #if DEBUG
1882 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Varstep2: ");
1883 debug_out_vector(LIF(sys),sys,&(sys->varstep2));
1884 #endif
1885 }
1886
1887
1888 /**
1889 Computes an update to the product B and varstep2.
1890 */
1891 static void calc_Bvarstep2( slv3_system_t sys){
1892 if( sys->Bvarstep2.accurate )
1893 return;
1894
1895 if ( !OPTIMIZING(sys) ) {
1896 zero_vector(&(sys->Bvarstep2));
1897 sys->Bvarstep2.norm2 = 0.0;
1898 } else {
1899 struct hessian_data *update;
1900 copy_vector(&(sys->varstep2),&(sys->Bvarstep2));
1901 for( update=sys->B; update != NULL; update = update->next ) {
1902 int32 col;
1903 real64 yv = inner_product( &(update->y),&(sys->varstep2) );
1904 real64 sBv = inner_product( &(update->Bs),&(sys->varstep2) );
1905 col = sys->Bvarstep2.rng->low;
1906 for( ; col <= sys->Bvarstep2.rng->high; col++ ) {
1907 sys->Bvarstep2.vec[col] += update->ys > 0.0 ?
1908 (update->y.vec[col])*yv/update->ys : 0.0;
1909 sys->Bvarstep2.vec[col] -= update->sBs > 0.0 ?
1910 (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1911 }
1912 }
1913 square_norm( &(sys->Bvarstep2) );
1914 }
1915 sys->Bvarstep2.accurate = TRUE;
1916 #if DEBUG
1917 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Bvarstep2: ");
1918 debug_out_vector(LIF(sys),sys,&(sys->Bvarstep2));
1919 #endif
1920 }
1921
1922
1923 /**
1924 Calculate the negative gradient direction of phi in the
1925 multipliers.
1926 */
1927 static void calc_mulstep1( slv3_system_t sys){
1928 if( sys->mulstep1.accurate )
1929 return;
1930
1931 if( !OPTIMIZING(sys) ) {
1932 zero_vector(&(sys->mulstep1));
1933 sys->mulstep1.norm2 = 0.0;
1934 } else {
1935 int32 row;
1936 row = sys->mulstep1.rng->low;
1937 for( ; row <= sys->mulstep1.rng->high; row++ )
1938 sys->mulstep1.vec[row] = -sys->residuals.vec[row];
1939 square_norm( &(sys->mulstep1) );
1940 }
1941 sys->mulstep1.accurate = TRUE;
1942 #if DEBUG
1943 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Mulstep1: ");
1944 debug_out_vector(LIF(sys),sys,&(sys->mulstep1));
1945 #endif
1946 }
1947
1948
1949 /**
1950 Calculate the mulstep2 direction of phi in the
1951 multipliers.
1952 */
1953 static void calc_mulstep2( slv3_system_t sys){
1954 if( sys->mulstep2.accurate )
1955 return;
1956
1957 if( !OPTIMIZING(sys) ) {
1958 zero_vector(&(sys->mulstep2));
1959 sys->mulstep2.norm2 = 0.0;
1960 } else {
1961 linsolqr_system_t lsys = sys->J.sys;
1962 int32 row;
1963 sys->J.rhs = linsolqr_get_rhs(lsys,2);
1964 mtx_zero_real64(sys->J.rhs,sys->cap);
1965 calc_rhs(sys, &(sys->Bvarstep2), -1.0, TRUE);
1966 calc_rhs(sys, &(sys->stationary), -1.0, TRUE);
1967 linsolqr_solve(lsys,sys->J.rhs);
1968 row = sys->mulstep2.rng->low;
1969 for( ; row <= sys->mulstep2.rng->high; row++ )
1970 sys->mulstep2.vec[row] = linsolqr_var_value
1971 (lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row));
1972 if (SAVLIN) {
1973 FILE *ldat;
1974 int32 ov;
1975 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1976 ldat=fopen(savlinfilename,"w");
1977 FPRINTF(ldat,
1978 "================= mulstep2rhs (orgcoled) itn %d =======\n",
1979 sys->s.iteration);
1980 debug_write_array(ldat,sys->J.rhs,sys->cap);
1981 FPRINTF(ldat,
1982 "================= mulstep2vars (orgrowed) ============\n");
1983 for(ov=0 ; ov < sys->cap; ov++ )
1984 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1985 fclose(ldat);
1986 }
1987 square_norm( &(sys->mulstep2) );
1988 }
1989 sys->mulstep2.accurate = TRUE;
1990 #if DEBUG
1991 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Mulstep2: ");
1992 debug_out_vector(LIF(sys),sys,&(sys->mulstep2));
1993 #endif
1994 }
1995
1996
1997 /**
1998 Computes the global minimizing function Phi.
1999 */
2000 static void calc_phi( slv3_system_t sys){
2001 if( !OPTIMIZING(sys) )
2002 sys->phi = 0.5*sys->residuals.norm2;
2003 else {
2004 sys->phi = sys->objective;
2005 sys->phi += inner_product( &(sys->multipliers),&(sys->residuals) );
2006 sys->phi += 0.5*RHO*sys->residuals.norm2;
2007 }
2008 }
2009
2010 /*------------------------------------------------------------------------------
2011 STEP CALCULATION STUFF
2012
2013 * OK. Here's where we compute the actual step to be taken. It will
2014 * be some linear combination of the 1st order and 2nd order steps.
2015 */
2016
2017 typedef real64 sym_2x2_t[3]; /* Stores symmetric 2x2 matrices */
2018
2019 struct parms_t {
2020 real64 low,high,guess; /* Used to search for parameter */
2021 };
2022
2023 struct calc_step_vars {
2024 sym_2x2_t coef1, coef2;
2025 real64 rhs[2]; /* RHS for 2x2 system */
2026 struct parms_t parms;
2027 real64 alpha1, alpha2;
2028 real64 error; /* Error between step norm and sys->maxstep */
2029 };
2030
2031 /**
2032 Calculates 2x2 system (coef1,coef2,rhs).
2033 */
2034 static void calc_2x2_system(slv3_system_t sys, struct calc_step_vars *vars){
2035 vars->coef1[0] = (2.0*sys->phi/sys->newton.norm2)*
2036 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
2037 vars->coef1[1] = 1.0;
2038 vars->coef1[2] = (sys->Jgamma.norm2/sys->gamma.norm2)*
2039 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
2040
2041 vars->coef2[0] = 1.0;
2042 vars->coef2[1] = 2.0*sys->phi/
2043 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
2044 vars->coef2[2] = 1.0;
2045
2046 vars->rhs[0] = 2.0*sys->phi/
2047 sys->maxstep/calc_sqrt_D0(sys->gamma.norm2);
2048 vars->rhs[1] = calc_sqrt_D0(sys->newton.norm2)/sys->maxstep;
2049 }
2050
2051 /**
2052 Determines alpha1 and alpha2 from the parameter (guess).
2053 */
2054 static void coefs_from_parm( slv3_system_t sys, struct calc_step_vars *vars){
2055
2056 sym_2x2_t coef; /* Actual coefficient matrix */
2057 real64 det; /* Determinant of coefficient matrix */
2058 int i;
2059
2060 for( i=0 ; i<3 ; ++i ) coef[i] =
2061 vars->coef1[i] + vars->parms.guess * vars->coef2[i];
2062 det = coef[0]*coef[2] - coef[1]*coef[1];
2063 if( det < 0.0 )
2064
2065 /* ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Unexpected negative determinant %f.", det); */
2066 fprintf(stderr,"Unexpected negative determinant %f.\n",det);
2067
2068 /* FPRINTF(MIF(sys),"%-40s ---> %g\n",
2069 " Unexpected negative determinant!",det); */
2070
2071 if( det <= DETZERO ) {
2072 /*
2073 * varstep2 and varstep1 are essentially parallel:
2074 * adjust length of either
2075 */
2076 vars->alpha2 = 0.0;
2077 vars->alpha1 = 1.0;
2078 } else {
2079 vars->alpha2 = (vars->rhs[0]*coef[2] - vars->rhs[1]*coef[1])/det;
2080 vars->alpha1 = (vars->rhs[1]*coef[0] - vars->rhs[0]*coef[1])/det;
2081 }
2082 }
2083
2084 /**
2085 Computes step vector length based on 1st order and 2nd order
2086 vectors and their coefficients.
2087 */
2088 static real64 step_norm2( slv3_system_t sys, struct calc_step_vars *vars){
2089 return sys->maxstep*sys->maxstep*
2090 (vars->alpha2 * vars->alpha2 +
2091 vars->alpha2 * vars->alpha1 * sys->phi/
2092 calc_sqrt_D0(sys->varstep2.norm2 + sys->mulstep2.norm2)/
2093 calc_sqrt_D0(sys->varstep1.norm2 + sys->mulstep1.norm2) +
2094 vars->alpha1 * vars->alpha1);
2095 }
2096
2097 /**
2098 Re-guesses the parameters based on step size vs. target value.
2099 */
2100 static void adjust_parms( slv3_system_t sys, struct calc_step_vars *vars){
2101 vars->error = (calc_sqrt_D0(step_norm2(sys,vars))/sys->maxstep) - 1.0;
2102 if( vars->error > 0.0 ) {
2103 /* Increase parameter (to decrease step length) */
2104 vars->parms.low = vars->parms.guess;
2105 vars->parms.guess = (vars->parms.high>3.0*vars->parms.guess)
2106 ? 2.0*vars->parms.guess
2107 : 0.5*(vars->parms.low + vars->parms.high);
2108 } else {
2109 /* Decrease parameter (to increase step norm) */
2110 vars->parms.high = vars->parms.guess;
2111 vars->parms.guess = 0.5*(vars->parms.low + vars->parms.high);
2112 }
2113 }
2114
2115 /**
2116 Computes the step based on the coefficients in vars.
2117 */
2118 static void compute_step( slv3_system_t sys, struct calc_step_vars *vars){
2119 int32 row,col;
2120 real64 tot1_norm2, tot2_norm2;
2121
2122 tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
2123 tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
2124 if( !sys->varstep.accurate ) {
2125 for( col=sys->varstep.rng->low ; col<=sys->varstep.rng->high ; ++col )
2126 if( (vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
2127 sys->varstep.vec[col] = sys->maxstep *
2128 sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2);
2129 } else if( (vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
2130 sys->varstep.vec[col] = sys->maxstep *
2131 sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2);
2132 } else if( (vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
2133 sys->varstep.vec[col] = sys->maxstep*
2134 (
2135 vars->alpha2*sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2) +
2136 vars->alpha1*sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2)
2137 );
2138 }
2139 sys->varstep.accurate = TRUE;
2140 }
2141 if( !sys->mulstep.accurate ) {
2142 for( row=sys->mulstep.rng->low ; row<=sys->mulstep.rng->high ; ++row )
2143 if( (vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
2144 sys->mulstep.vec[row] = sys->maxstep *
2145 sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2);
2146 } else if( (vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
2147 sys->mulstep.vec[row] = sys->maxstep *
2148 sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2);
2149 } else if( (vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
2150 sys->mulstep.vec[row] = sys->maxstep*
2151 (
2152 vars->alpha2*sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2) +
2153 vars->alpha1*sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2)
2154 );
2155 }
2156 sys->mulstep.accurate = TRUE;
2157 }
2158 #if DEBUG
2159 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Varstep: ");
2160 debug_out_vector(LIF(sys),sys,&(sys->varstep));
2161 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Mulstep: ");
2162 debug_out_vector(LIF(sys),sys,&(sys->mulstep));
2163 #endif
2164 }
2165
2166
2167 /**
2168 Calculates step vector, based on sys->maxstep, and the varstep2/
2169 varstep1 and mulstep2/mulstep1 vectors. Nothing is assumed to be
2170 calculated, except the weights and the jacobian (scaled). Also,
2171 the step is not checked for legitimacy.
2172 NOTE: the step is scaled.
2173 */
2174 static void calc_step( slv3_system_t sys, int minor){
2175
2176 struct calc_step_vars vars;
2177 real64 tot1_norm2, tot2_norm2;
2178
2179 if( sys->varstep.accurate && sys->mulstep.accurate )
2180 return;
2181 if (SHOW_LESS_IMPT) {
2182 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\n%-40s ---> %d\n", " Step trial",minor);
2183 }
2184
2185 tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
2186 tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
2187 if( (tot1_norm2 == 0.0) && (tot2_norm2 == 0.0) ) {
2188 /* Take no step at all */
2189 vars.alpha1 = 0.0;
2190 vars.alpha2 = 0.0;
2191 sys->maxstep = 0.0;
2192 sys->varstep.norm2 = 0.0;
2193 sys->mulstep.norm2 = 0.0;
2194
2195 } else if( (tot2_norm2 > 0.0) && OPTIMIZING(sys) ) {
2196 /* Stay in varstep2 direction */
2197 vars.alpha1 = 0.0;
2198 vars.alpha2 = 1.0;
2199 sys->maxstep = MIN(sys->maxstep,calc_sqrt_D0(tot2_norm2));
2200 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2201 sys->varstep2.norm2/tot2_norm2;
2202 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2203 sys->mulstep2.norm2/tot2_norm2;
2204
2205 } else if( (tot2_norm2>0.0)&&(calc_sqrt_D0(tot2_norm2)<=sys->maxstep) ) {
2206 /* Attempt step in varstep2 direction */
2207 vars.alpha1 = 0.0;
2208 vars.alpha2 = 1.0;
2209 sys->maxstep = calc_sqrt_D0(tot2_norm2);
2210 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2211 sys->varstep2.norm2/tot2_norm2;
2212 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2213 sys->mulstep2.norm2/tot2_norm2;
2214
2215 } else if( (tot2_norm2==0.0 || sys->s.block.current_size==1) &&
2216 (tot1_norm2 > 0.0) ) {
2217 /* Attempt step in varstep1 direction */
2218 vars.alpha1 = 1.0;
2219 vars.alpha2 = 0.0;
2220 if ( (sys->gamma.norm2/sys->Jgamma.norm2)*
2221 calc_sqrt_D0(sys->gamma.norm2) <= sys->maxstep )
2222 sys->maxstep = (sys->gamma.norm2/sys->Jgamma.norm2)*
2223 calc_sqrt_D0(sys->gamma.norm2);
2224 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2225 sys->varstep1.norm2/tot1_norm2;
2226 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2227 sys->mulstep1.norm2/tot1_norm2;
2228
2229 } else {
2230 /* Attempt step in varstep1-varstep2 direction */
2231 vars.parms.low = 0.0;
2232 vars.parms.high = MAXDOUBLE;
2233 vars.parms.guess = 1.0;
2234 calc_2x2_system(sys,&vars);
2235 do {
2236 coefs_from_parm(sys, &vars);
2237 adjust_parms(sys, &vars);
2238 } while( fabs(vars.error) > STEPSIZEERR_MAX &&
2239 vars.parms.high - vars.parms.low > PARMRNG_MIN );
2240 if (SHOW_LESS_IMPT) {
2241 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2242 " parameter high", vars.parms.high);
2243 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2244 " parameter low", vars.parms.low);
2245 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2246 " Error in step length", vars.error);
2247 }
2248 sys->varstep.norm2 = step_norm2(sys, &vars);
2249 sys->mulstep.norm2 = 0.0;
2250 }
2251
2252 if (SHOW_LESS_IMPT) {
2253 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", " Alpha1 coefficient (normalized)",
2254 vars.alpha1);
2255 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", " Alpha2 coefficient (normalized)",
2256 vars.alpha2);
2257 }
2258 compute_step(sys,&vars);
2259 return;
2260
2261 }
2262
2263 /*------------------------------------------------------------------------------
2264 VARIABLE VALUES MAINTENANCE
2265 */
2266
2267 /**
2268 Restores the values of the variables before applying
2269 a step.
2270 */
2271 static void restore_variables( slv3_system_t sys){
2272 int32 col;
2273 real64 *vec;
2274 vec = (sys->nominals.vec);
2275 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
2276 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2277 var_set_value(var,sys->variables.vec[col]*vec[col]);
2278 }
2279 }
2280
2281
2282 /**
2283 Calculates the maximum fraction of the step which can be
2284 taken without going out of bounds. If the entire step can be
2285 taken, 1.0 is returned. Otherwise a value less than 1 is
2286 returned. It is assumed that the current variable values
2287 are within their bounds. The step must be calculated.
2288 */
2289 static real64 required_coef_to_stay_inbounds( slv3_system_t sys){
2290 real64 mincoef;
2291 int32 col;
2292 real64 *vec;
2293 vec = (sys->nominals.vec);
2294
2295 if( sys->p.ignore_bounds )
2296 return(1.0);
2297
2298 mincoef = 1.0;
2299 for( col=sys->varstep.rng->low; col <= sys->varstep.rng->high; col++ ) {
2300 struct var_variable *var;
2301 real64 coef,dx,val,bnd;
2302 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2303 coef = 1.0;
2304 dx = sys->varstep.vec[col] * vec[col];
2305 bnd = var_upper_bound(var);
2306 if( (val=var_value(var)) + dx > bnd )
2307 coef = MIN((bnd-val)/dx, 1.0);
2308 bnd = var_lower_bound(var);
2309 if( val + dx < bnd )
2310 coef = MIN((bnd-val)/dx, 1.0);
2311 if( coef < mincoef )
2312 mincoef = coef;
2313 }
2314 return(mincoef);
2315 }
2316
2317
2318 /**
2319 Adds sys->varstep to the variable values in block: projecting
2320 near bounds.
2321 */
2322 static void apply_step( slv3_system_t sys){
2323 FILE *lif = LIF(sys);
2324 int nproj = 0;
2325 real64 bounds_coef = 1.0;
2326 int32 col;
2327 real64 *vec;
2328 vec = (sys->nominals.vec);
2329
2330 if (TRUNCATE && (!sys->p.ignore_bounds))
2331 bounds_coef = required_coef_to_stay_inbounds(sys);
2332
2333 for( col=sys->varstep.rng->low; col <= sys->varstep.rng->high; col++ ) {
2334 struct var_variable *var;
2335 real64 dx,val,bnd;
2336 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2337 dx = vec[col]*sys->varstep.vec[col];
2338 val = var_value(var);
2339 if (bounds_coef < 1.0) {
2340 dx = dx*TOWARD_BOUNDS*bounds_coef;
2341 sys->varstep.vec[col] = dx/vec[col];
2342 } else {
2343 if( !sys->p.ignore_bounds ) {
2344 if( val + dx > (bnd=var_upper_bound(var)) ) {
2345 dx = TOWARD_BOUNDS*(bnd-val);
2346 sys->varstep.vec[col] = dx/vec[col];
2347 if (SHOW_LESS_IMPT) {
2348 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> ",
2349 " Variable projected to upper bound");
2350 print_var_name(lif,sys,var); PUTC('\n',lif);
2351 }
2352 ++nproj;
2353 } else if( val + dx < (bnd=var_lower_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 lower bound");
2359 print_var_name(lif,sys,var); PUTC('\n',lif);
2360 }
2361 ++nproj;
2362 }
2363 }
2364 }
2365 var_set_value(var,val+dx);
2366 }
2367
2368 if( !sys->p.ignore_bounds ) {
2369 if (nproj > 0) {
2370 square_norm(&(sys->varstep));
2371 sys->progress = calc_sqrt_D0
2372 (calc_sqrt_D0((sys->varstep.norm2 + sys->mulstep.norm2)*
2373 (sys->varstep1.norm2 + sys->mulstep1.norm2)));
2374 if (SHOW_LESS_IMPT) {
2375 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", " Projected step length (scaled)",
2376 calc_sqrt_D0(sys->varstep.norm2 + sys->mulstep.norm2));
2377 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2378 " Projected progress", sys->progress);
2379 }
2380 }
2381 if (bounds_coef < 1.0) {
2382 square_norm(&(sys->varstep));
2383 if (SHOW_LESS_IMPT) {
2384 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2385 " Truncated step length (scaled)",
2386 calc_sqrt_D0(sys->varstep.norm2 + sys->mulstep.norm2));
2387 }
2388 sys->progress = calc_sqrt_D0
2389 (calc_sqrt_D0((sys->varstep.norm2 + sys->mulstep.norm2)*
2390 (sys->varstep1.norm2 + sys->mulstep1.norm2)));
2391 if (SHOW_LESS_IMPT) {
2392 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2393 " Truncated progress", sys->progress);
2394 }
2395 }
2396 }
2397
2398 /* Allow weighted residuals to be recalculated at new point */
2399 sys->residuals.accurate = FALSE;
2400
2401 return;
2402 }
2403
2404 /**
2405 This function should be called when the step is accepted.
2406 */
2407 static void step_accepted( slv3_system_t sys){
2408 /* Maintain update status on jacobian and weights */
2409 if (--(sys->update.jacobian) <= 0)
2410 sys->J.accurate = FALSE;
2411
2412 sys->ZBZ.accurate = FALSE;
2413 sys->variables.accurate = FALSE;
2414 sys->gradient.accurate = FALSE;
2415 sys->multipliers.accurate = FALSE;
2416 sys->stationary.accurate = FALSE;
2417 sys->newton.accurate = FALSE;
2418 sys->Bnewton.accurate = FALSE;
2419 sys->nullspace.accurate = FALSE;
2420 sys->gamma.accurate = FALSE;
2421 sys->Jgamma.accurate = FALSE;
2422 sys->varstep1.accurate = FALSE;
2423 sys->Bvarstep1.accurate = FALSE;
2424 sys->varstep2.accurate = FALSE;
2425 sys->Bvarstep2.accurate = FALSE;
2426 sys->mulstep1.accurate = FALSE;
2427 sys->mulstep2.accurate = FALSE;
2428 sys->varstep.accurate = FALSE;
2429 sys->mulstep.accurate = FALSE;
2430
2431 if( !OPTIMIZING(sys) ) {
2432 sys->ZBZ.accurate = TRUE;
2433 sys->gradient.accurate = TRUE;
2434 sys->multipliers.accurate = TRUE;
2435 sys->stationary.accurate = TRUE;
2436 sys->Bnewton.accurate = TRUE;
2437 sys->nullspace.accurate = TRUE;
2438 sys->Bvarstep1.accurate = TRUE;
2439 sys->Bvarstep2.accurate = TRUE;
2440 }
2441 }
2442
2443 /**
2444 This function changes sys->maxstep to the given number and should be
2445 called whenever sys->maxstep is to be changed.
2446 */
2447 static void change_maxstep( slv3_system_t sys, real64 maxstep){
2448 sys->maxstep = maxstep;
2449 sys->varstep.accurate = FALSE;
2450 if( OPTIMIZING(sys) ) sys->mulstep.accurate = FALSE;
2451 }
2452
2453
2454 /*------------------------------------------------------------------------------
2455 BLOCK ROUTINES
2456 */
2457
2458 /**
2459 Returns TRUE if the current block is feasible, FALSE otherwise.
2460 It is assumed that the residuals have been computed.
2461 */
2462 static boolean block_feasible( slv3_system_t sys){
2463 int32 row;
2464
2465 if( !sys->s.calc_ok )
2466 return(FALSE);
2467
2468 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2469 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2470 if( !rel_satisfied(rel) ) return FALSE;
2471 }
2472 return TRUE;
2473 }
2474
2475 /**
2476 Moves on to the next block, updating all of the solver information.
2477 To move to the first block, set sys->s.block.current_block to -1 before
2478 calling. If already at the last block, then sys->s.block.current_block
2479 will equal the number of blocks and the system will be declared
2480 converged. Otherwise, the residuals for the new block will be computed
2481 and sys->s.calc_ok set according.
2482 */
2483 static void move_to_next_block( slv3_system_t sys){
2484 struct var_variable *var;
2485 struct rel_relation *rel;
2486 int32 row;
2487 int32 col;
2488 int32 ci;
2489
2490 if( sys->s.block.current_block >= 0 ) {
2491
2492
2493 /* Record cost accounting info here. */
2494 ci=sys->s.block.current_block;
2495 sys->s.cost[ci].size = sys->s.block.current_size;
2496 sys->s.cost[ci].iterations = sys->s.block.iteration;
2497 sys->s.cost[ci].funcs = sys->s.block.funcs;
2498 sys->s.cost[ci].jacs = sys->s.block.jacs;
2499 sys->s.cost[ci].functime = sys->s.block.functime;
2500 sys->s.cost[ci].jactime = sys->s.block.jactime;
2501 sys->s.cost[ci].time = sys->s.block.cpu_elapsed;
2502 sys->s.cost[ci].resid = sys->s.block.residual;
2503
2504 /* De-initialize previous block */
2505 if (SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
2506 LIFDS)) {
2507 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Block %d converged.\n",
2508 sys->s.block.current_block);
2509 }
2510 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
2511 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2512 var_set_in_block(var,FALSE);
2513 }
2514 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2515 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2516 rel_set_in_block(rel,FALSE);
2517 }
2518 sys->s.block.previous_total_size += sys->s.block.current_size;
2519 }
2520
2521 sys->s.block.current_block++;
2522 if( sys->s.block.current_block < sys->s.block.number_of ) {
2523 boolean ok;
2524
2525 /* Initialize next block */
2526 if( OPTIMIZING(sys) ) {
2527 mtx_region(&(sys->J.reg), 0, sys->rank-1, 0, sys->vused-1 );
2528 } else {
2529 sys->J.reg =
2530 (slv_get_solvers_blocks(SERVER))->block[sys->s.block.current_block];
2531 }
2532
2533 row = sys->J.reg.row.high - sys->J.reg.row.low + 1;
2534 col = sys->J.reg.col.high - sys->J.reg.col.low + 1;
2535 sys->s.block.current_size = MAX(row,col);
2536
2537 sys->s.block.iteration = 0;
2538 sys->s.block.cpu_elapsed = 0.0;
2539 sys->s.block.functime = 0.0;
2540 sys->s.block.jactime = 0.0;
2541 sys->s.block.funcs = 0;
2542 sys->s.block.jacs = 0;
2543
2544 if(SHOW_LESS_IMPT && (LIFDS ||
2545 sys->s.block.current_size > 1)) {
2546 debug_delimiter(LIF(sys));
2547 debug_delimiter(LIF(sys));
2548 }
2549 if(SHOW_LESS_IMPT && LIFDS) {
2550 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\n%-40s ---> %d in [%d..%d]\n",
2551 "Current block number", sys->s.block.current_block,
2552 0, sys->s.block.number_of-1);
2553 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d\n", "Current block size",
2554 sys->s.block.current_size);
2555 }
2556 sys->s.calc_ok = TRUE;
2557
2558 if( !(ok = calc_objective(sys)) ) {
2559 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);
2560 FPRINTF(MIF(sys),"Objective calculation errors detected.\n");
2561 error_reporter_end_flush();
2562 }
2563 if(SHOW_LESS_IMPT && sys->obj) {
2564 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", "Objective", sys->objective);
2565 }
2566 sys->s.calc_ok = sys->s.calc_ok && ok;
2567
2568 if (!(sys->p.ignore_bounds) ) {
2569 slv_insure_bounds(SERVER, sys->J.reg.col.low,
2570 sys->J.reg.col.high,MIF(sys));
2571 }
2572
2573 sys->residuals.accurate = FALSE;
2574 if( !(ok = calc_residuals(sys)) ) {
2575 /* error_reporter will have been called somewhere else already */
2576 CONSOLE_DEBUG("Residual calculation errors detected in move_to_next_block.");
2577 }
2578 if( SHOW_LESS_IMPT &&
2579 (sys->s.block.current_size >1 ||
2580 LIFDS) ) {
2581 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", "Residual norm (unscaled)",
2582 sys->s.block.residual);
2583 }
2584 sys->s.calc_ok = sys->s.calc_ok && ok;
2585
2586 /* Must be updated as soon as required */
2587 sys->J.accurate = FALSE;
2588 sys->update.weights = 0;
2589 sys->update.nominals = 0;
2590 sys->update.relnoms = 0;
2591 sys->update.iterative = 0;
2592 sys->ZBZ.accurate = FALSE;
2593 sys->variables.accurate = FALSE;
2594 sys->gradient.accurate = FALSE;
2595 sys->multipliers.accurate = FALSE;
2596 sys->stationary.accurate = FALSE;
2597 sys->newton.accurate = FALSE;
2598 sys->Bnewton.accurate = FALSE;
2599 sys->nullspace.accurate = FALSE;
2600 sys->gamma.accurate = FALSE;
2601 sys->Jgamma.accurate = FALSE;
2602 sys->varstep1.accurate = FALSE;
2603 sys->Bvarstep1.accurate = FALSE;
2604 sys->varstep2.accurate = FALSE;
2605 sys->Bvarstep2.accurate = FALSE;
2606 sys->mulstep1.accurate = FALSE;
2607 sys->mulstep2.accurate = FALSE;
2608 sys->varstep.accurate = FALSE;
2609 sys->mulstep.accurate = FALSE;
2610
2611 if( !OPTIMIZING(sys) ) {
2612 sys->ZBZ.accurate = TRUE;
2613 sys->gradient.accurate = TRUE;
2614 sys->multipliers.accurate = TRUE;
2615 sys->stationary.accurate = TRUE;
2616 sys->Bnewton.accurate = TRUE;
2617 sys->nullspace.accurate = TRUE;
2618 sys->Bvarstep1.accurate = TRUE;
2619 sys->Bvarstep2.accurate = TRUE;
2620 }
2621
2622 } else {
2623 boolean ok;
2624 /*
2625 * Before we claim convergence, we must check if we left behind
2626 * some unassigned relations. If and only if they happen to be
2627 * satisfied at the current point, convergence has been obtained.
2628 *
2629 * Also insures that all included relations have valid residuals.
2630 * Included inequalities will have correct residuals.
2631 * Unsatisfied included inequalities cause inconsistency.
2632 *
2633 * This of course ignores that fact an objective function might
2634 * be present. Then, feasibility isn't enough, is it now.
2635 */
2636 if( sys->s.struct_singular ) {
2637 /* black box w/singletons provoking bug here, maybe */
2638 sys->s.block.current_size = sys->rused - sys->rank;
2639 if(SHOW_LESS_IMPT) {
2640 debug_delimiter(LIF(sys));
2641 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d\n", "Unassigned Relations",
2642 sys->s.block.current_size);
2643 }
2644 sys->J.reg.row.low = sys->J.reg.col.low = sys->rank;
2645 sys->J.reg.row.high = sys->J.reg.col.high = sys->rused - 1;
2646 sys->residuals.accurate = FALSE;
2647 if( !(ok=calc_residuals(sys)) ) {
2648 FPRINTF(MIF(sys),
2649 "Residual calculation errors detected in leftover equations.\n");
2650 }
2651 if(SHOW_LESS_IMPT) {
2652 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", "Residual norm (unscaled)",
2653 sys->s.block.residual);
2654 }
2655 if( block_feasible(sys) ) {
2656 if(SHOW_LESS_IMPT) {
2657 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nUnassigned relations ok. You lucked out.\n");
2658 }
2659 sys->s.converged = TRUE;
2660 } else {
2661 if(SHOW_LESS_IMPT) {
2662 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nProblem inconsistent: %s.\n",
2663 "Unassigned relations not satisfied");
2664 }
2665 sys->s.inconsistent = TRUE;
2666 }
2667 if(SHOW_LESS_IMPT) {
2668 debug_delimiter(LIF(sys));
2669 }
2670 } else {
2671 sys->s.converged = TRUE;
2672 }
2673 /* nearly done checking. Must verify included inequalities if
2674 we think equalities are ok. */
2675 if (sys->s.converged) {
2676 sys->s.inconsistent=(!calc_inequalities(sys));
2677 }
2678 }
2679 }
2680
2681 /**
2682 Calls the appropriate reorder function on a block
2683 */
2684 static void reorder_new_block(slv3_system_t sys){
2685 int32 method;
2686 if( sys->s.block.current_block < sys->s.block.number_of ) {
2687 if (strcmp(REORDER_OPTION,"SPK1") == 0) {
2688 method = 2;
2689 } else {
2690 method = 1;
2691 }
2692
2693 if( sys->s.block.current_block <= sys->s.block.current_reordered_block &&
2694 sys->s.cost[sys->s.block.current_block].reorder_method == method &&
2695 sys->s.block.current_block >= 0 ) {
2696 #if DEBUG
2697 FPRINTF(stderr,"YOU JUST AVOIDED A REORDERING\n");
2698 #endif
2699 slv_set_up_block(SERVER,sys->s.block.current_block);
2700 /* tell linsol to bless it and get on with things */
2701 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
2702 return; /*must have been reordered since last system build*/
2703 }
2704
2705 /* Let the slv client function take care of reordering things
2706 * and setting in block flags.
2707 */
2708 if (strcmp(REORDER_OPTION,"SPK1") == 0) {
2709 sys->s.cost[sys->s.block.current_block].reorder_method = 2;
2710 slv_spk1_reorder_block(SERVER,sys->s.block.current_block,1);
2711 } else if (strcmp(REORDER_OPTION,"TEAR_DROP") == 0) {
2712 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2713 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2714 CUTOFF,
2715 0,mtx_SPK1);
2716 /* khack: try tspk1 for transpose case */
2717 } else if (strcmp(REORDER_OPTION,"OVER_TEAR") == 0) {
2718 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2719 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2720 CUTOFF,
2721 1,mtx_SPK1);
2722 } else {
2723 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2724 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);
2725 FPRINTF(MIF(sys),"QRSlv called with unknown reorder option\n");
2726 FPRINTF(MIF(sys),"QRSlv using single edge tear drop (TEAR_DROP).\n");
2727 error_reporter_end_flush();
2728
2729 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2730 CUTOFF,0,mtx_SPK1);
2731 }
2732 /* tell linsol to bless it and get on with things */
2733 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
2734 if (sys->s.block.current_block > sys->s.block.current_reordered_block) {
2735 sys->s.block.current_reordered_block = sys->s.block.current_block;
2736 }
2737 }
2738 }
2739
2740 /**
2741 Moves to next unconverged block, assuming that the current block has
2742 converged (or is -1, to start).
2743 */
2744 static void find_next_unconverged_block( slv3_system_t sys){
2745
2746 do{
2747 move_to_next_block(sys);
2748 #if DEBUG
2749 debug_out_var_values(stderr,sys);
2750 debug_out_rel_residuals(stderr,sys);
2751 #endif
2752 }while( !sys->s.converged && block_feasible(sys) && !OPTIMIZING(sys));
2753
2754 reorder_new_block(sys);
2755 }
2756
2757 /*------------------------------------------------------------------------------
2758 ITERATION BEGIN/END ROUTINES
2759 */
2760
2761 /**
2762 Prepares sys for entering an iteration, increasing the iteration counts
2763 and starting the clock.
2764 */
2765 static void iteration_begins( slv3_system_t sys){
2766 sys->clock = tm_cpu_time();
2767 ++(sys->s.block.iteration);
2768 ++(sys->s.iteration);
2769 if(SHOW_LESS_IMPT&& (sys->s.block.current_size >1 ||
2770 LIFDS)) {
2771 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\n%-40s ---> %d\n",
2772 "Iteration", sys->s.block.iteration);
2773 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d\n",
2774 "Total iteration", sys->s.iteration);
2775 }
2776 }
2777
2778 /**
2779 Prepares sys for exiting an iteration, stopping the clock and recording
2780 the cpu time.
2781 */
2782 static void iteration_ends( slv3_system_t sys){
2783 double cpu_elapsed; /* elapsed this iteration */
2784
2785 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
2786 sys->s.block.cpu_elapsed += cpu_elapsed;
2787 sys->s.cpu_elapsed += cpu_elapsed;
2788 if(SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
2789 LIFDS)) {
2790 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2791 "Elapsed time", sys->s.block.cpu_elapsed);
2792 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2793 "Total elapsed time", sys->s.cpu_elapsed);
2794 }
2795 }
2796
2797 /**
2798 Updates the solver status.
2799 */
2800 static void update_status( slv3_system_t sys){
2801 boolean unsuccessful;
2802
2803 if( !sys->s.converged ) {
2804 sys->s.time_limit_exceeded =
2805 (sys->s.block.cpu_elapsed >= TIME_LIMIT);
2806 sys->s.iteration_limit_exceeded =
2807 (sys->s.block.iteration >= ITER_LIMIT);
2808 }
2809
2810 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
2811 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
2812
2813 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
2814 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
2815 }
2816
2817 static
2818 int32 slv3_get_default_parameters(slv_system_t server, SlvClientToken asys,
2819 slv_parameters_t *parameters)
2820 {
2821 slv3_system_t sys = NULL;
2822 union parm_arg lo,hi,val;
2823 struct slv_parameter *new_parms = NULL;
2824 int32 make_macros = 0;
2825
2826 static char *factor_names[] = {
2827 "SPK1/RANKI","SPK1/RANKI+ROW",
2828 "Fast-SPK1/RANKI","Fast-SPK1/RANKI+ROW",
2829 "Fastest-SPK1/MR-RANKI","CondQR","CPQR"
2830 /* ,"GAUSS","GAUSS_EASY" currently only works for ken */
2831 };
2832 static char *reorder_names[] = {
2833 "SPK1","TEAR_DROP","OVER_TEAR"
2834 };
2835 static char *converge_names[] = {
2836 "ABSOLUTE","RELNOM_SCALE"
2837 };
2838 static char *scaling_names[] = {
2839 "NONE","ROW_2NORM","RELNOM"
2840 };
2841
2842 if (server != NULL && asys != NULL) {
2843 sys = SLV3(asys);
2844 make_macros = 1;
2845 }
2846
2847 #ifndef NDEBUG /* keep purify from whining on UMR */
2848 lo.argr = hi.argr = val.argr = 0.0;
2849 #endif
2850
2851 if (parameters->parms == NULL) {
2852 /* an external client wants our parameter list.
2853 * an instance of slv3_system_structure has this pointer
2854 * already set in slv3_create
2855 */
2856 new_parms = ASC_NEW_ARRAY_OR_NULL(struct slv_parameter,slv3_PA_SIZE);
2857 if (new_parms == NULL) {
2858 return