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