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