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