/[ascend]/branches/ksenija2/solvers/qrslv/qrslv.c
ViewVC logotype

Contents of /branches/ksenija2/solvers/qrslv/qrslv.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2869 - (show annotations) (download) (as text)
Tue Mar 24 13:25:23 2015 UTC (7 years, 3 months ago) by jpye
File MIME type: text/x-csrc
File size: 135022 byte(s)
add debug output of QRSlv vars are rels list (active,included etc) to idaboundary.c.
change solardynamics model for minimum-size boundary NLA problem. now it doesn't work.
some pointless tinkering with formatting.

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