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