/[ascend]/trunk/solvers/qrslv/qrslv.c
ViewVC logotype

Contents of /trunk/solvers/qrslv/qrslv.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1555 - (show annotations) (download) (as text)
Tue Jul 24 14:16:45 2007 UTC (13 years, 3 months ago) by jpye
File MIME type: text/x-csrc
File size: 134811 byte(s)
Bit more work on IPOPT presolve and hessian calculation.
Added export of relman_jacobian_count.
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3 Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
4 Copyright (C) 1993 Joseph Zaher
5 Copyright (C) 1990 Karl Michael Westerberg
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2, or (at your option)
10 any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 59 Temple Place - Suite 330,
20 Boston, MA 02111-1307, USA.
21 *//**
22 @file
23 QRSLV solver module for ASCEND.
24 *//*
25 by Karl Michael Westerberg
26 Created: 2/6/90
27 Last *CVS* version ballan 2000/01/25 02:27:32
28 */
29
30 #include <math.h>
31 #include <stdarg.h>
32
33 #define ASC_BUILDING_INTERFACE
34
35 #include <utilities/config.h>
36 #include <utilities/ascConfig.h>
37 #ifdef ASC_SIGNAL_TRAPS
38 # include <utilities/ascSignal.h>
39 #endif
40
41 #include <utilities/ascMalloc.h>
42 #include <utilities/set.h>
43 #include <general/mathmacros.h>
44 #include <general/tm_time.h>
45 #include <utilities/mem.h>
46 #include <utilities/ascPanic.h>
47 #include <general/list.h>
48
49 #include <linear/mtx_vector.h>
50
51 #include <system/calc.h>
52 #include <system/slv_stdcalls.h>
53 #include <system/relman.h>
54 #include <system/block.h>
55 #include <solver/solver.h>
56
57 #define CANOPTIMIZE FALSE
58 /**< TRUE iff optimization code completed, meaning relman_diff fixed. */
59
60 #define DEBUG FALSE
61 /**< makes lots of extra spew */
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 FILE *fmtx = NULL;
1247
1248 linsolqr_system_t lsys = sys->J.sys;
1249 FILE *fp = LIF(sys);
1250
1251 oldtiming = g_linsolqr_timing;
1252 g_linsolqr_timing =SLV_PARAM_BOOL(&(sys->p),LINTIME);
1253 linsolqr_factor(lsys,sys->J.fm); /* factor */
1254 g_linsolqr_timing = oldtiming;
1255
1256 if(OPTIMIZING(sys)){
1257 CONSOLE_DEBUG("OPTIMISING");
1258 /* need things for nullspace move. don't care about
1259 * dependency coefficiency in any circumstances at present.
1260 */
1261 linsolqr_calc_col_dependencies(lsys);
1262 set_null(sys->J.relpivots,sys->cap);
1263 set_null(sys->J.varpivots,sys->cap);
1264 linsolqr_get_pivot_sets(lsys,sys->J.relpivots,sys->J.varpivots);
1265 }
1266
1267 sys->J.rank = linsolqr_rank(lsys);
1268 sys->J.singular = FALSE;
1269 row_rank_defect = sys->J.reg.row.high - sys->J.reg.row.low+1 - sys->J.rank;
1270 if(row_rank_defect > 0) {
1271 int32 row,krow;
1272 mtx_sparse_t *uprows=NULL;
1273 sys->J.singular = TRUE;
1274 uprows = linsolqr_unpivoted_rows(lsys);
1275 if(uprows !=NULL){
1276 for( krow=0; krow < uprows->len ; krow++ ) {
1277 int32 org_row;
1278 struct rel_relation *rel;
1279
1280 org_row = uprows->idata[krow];
1281 row = mtx_org_to_row(sys->J.mtx,org_row);
1282 rel = sys->rlist[org_row];
1283
1284 ERROR_REPORTER_START_HERE(ASC_PROG_WARNING);
1285 FPRINTF(ASCERR,"Relation '");
1286 print_rel_name(stderr,sys,rel);
1287 FPRINTF(ASCERR,"' is not pivoted.");
1288 error_reporter_end_flush();
1289
1290 /*
1291 * assign zeros to the corresponding weights
1292 * so that subsequent calls to "scale_residuals"
1293 * will only measure the pivoted equations.
1294 */
1295 sys->weights.vec[row] = 0.0;
1296 sys->residuals.vec[row] = 0.0;
1297 sys->residuals.accurate = FALSE;
1298 mtx_mult_row(sys->J.mtx,row,0.0,&(sys->J.reg.col));
1299 }
1300 mtx_destroy_sparse(uprows);
1301 }
1302 if(!sys->residuals.accurate ) {
1303 square_norm( &(sys->residuals) );
1304 sys->residuals.accurate = TRUE;
1305 sys->update.weights = 0; /* re-compute weights next iteration. */
1306 }
1307
1308 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Row rank defect = %d (block = %d rows, rank = %d)"
1309 ,row_rank_defect
1310 ,sys->J.reg.row.high - sys->J.reg.row.low+1
1311 ,sys->J.rank
1312 );
1313
1314 #ifdef ASC_WITH_MMIO
1315 #define QRSLV_MMIO_FILE "qrslvmmio.mtx"
1316 /* #define QRSLV_MMIO_WHOLE */
1317 if((fmtx = fopen(QRSLV_MMIO_FILE,"w"))){
1318 #ifdef QRSLV_MMIO_WHOLE
1319 mtx_write_region_mmio(fmtx, sys->J.mtx, mtx_ENTIRE_MATRIX);
1320 #else
1321 mtx_write_region_mmio(fmtx, sys->J.mtx, &(sys->J.reg));
1322 #endif
1323 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Wrote matrix to '%s' (EXPERIMENTAL!)",QRSLV_MMIO_FILE);
1324 fclose(fmtx);
1325 }else{
1326 ERROR_REPORTER_HERE(ASC_PROG_ERR,
1327 "Unable to write matrix to '%s' (couldn't open for writing)",QRSLV_MMIO_FILE
1328 );
1329 }
1330 #endif
1331 }
1332
1333 if(sys->J.rank < sys->J.reg.col.high-sys->J.reg.col.low+1 ) {
1334 int32 col,kcol;
1335 mtx_sparse_t *upcols=NULL;
1336 if(NOTNULL(upcols)) {
1337 for( kcol=0; upcols != NULL && kcol < upcols->len ; kcol++ ) {
1338 int32 org_col;
1339 struct var_variable *var;
1340
1341 org_col = upcols->idata[kcol];
1342 col = mtx_org_to_col(sys->J.mtx,org_col);
1343 var = sys->vlist[org_col];
1344 FPRINTF(fp,"%-40s ---> ","Variable not pivoted");
1345 print_var_name(fp,sys,var);
1346 PUTC('\n',fp);
1347 /*
1348 * If we're not optimizing (everything should be
1349 * pivotable) or this was one of the dependent variables,
1350 * consider this variable as if it were fixed.
1351 */
1352 if(col <= sys->J.reg.col.high - sys->ZBZ.order ) {
1353 mtx_mult_col(sys->J.mtx,col,0.0,&(sys->J.reg.row));
1354 }
1355 }
1356 mtx_destroy_sparse(upcols);
1357 }
1358 }
1359 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
1360 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d (%s)\n","Jacobian rank", sys->J.rank,
1361 sys->J.singular ? "deficient":"full");
1362 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n","Smallest pivot",
1363 linsolqr_smallest_pivot(sys->J.sys));
1364 }
1365 return row_rank_defect;
1366 }
1367
1368 /**
1369 Updates the reduced hessian matrix.
1370 if !OPTIMIZING just sets zbz.accurate true and returns.
1371 */
1372 static void calc_ZBZ(qrslv_system_t sys){
1373 mtx_coord_t nz;
1374
1375 if(sys->ZBZ.accurate ) return;
1376
1377 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1378 for( nz.col = 0; nz.col <= nz.row; nz.col++ ) {
1379 int32 col, depr, depc;
1380 col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
1381 depr = mtx_col_to_org(sys->J.mtx,col);
1382 col = nz.col+sys->J.reg.col.high+1-sys->ZBZ.order;
1383 depc = mtx_col_to_org(sys->J.mtx,col);
1384 sys->ZBZ.mtx[nz.row][nz.col] = (nz.row==nz.col ? 1.0 : 0.0);
1385 col = sys->J.reg.col.low;
1386 for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
1387 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1388 if(set_is_member(sys->J.varpivots,ind) )
1389 sys->ZBZ.mtx[nz.row][nz.col] +=
1390 (-linsolqr_org_col_dependency(sys->J.sys,depr,ind))*
1391 (-linsolqr_org_col_dependency(sys->J.sys,depc,ind));
1392 }
1393 if(nz.row != nz.col ) {
1394 sys->ZBZ.mtx[nz.col][nz.row] =
1395 sys->ZBZ.mtx[nz.row][nz.col];
1396 }
1397 }
1398 }
1399 if(OPTIMIZING(sys)){
1400 struct hessian_data *update;
1401 for( update=sys->B; update != NULL; update = update->next ) {
1402 for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1403 int32 col, dep;
1404 col = nz.row + sys->J.reg.col.high + 1 - sys->ZBZ.order;
1405 dep = mtx_col_to_org(sys->J.mtx,col);
1406 sys->ZBZ.Zy[nz.row] = update->y.vec[col];
1407 sys->ZBZ.ZBs[nz.row] = update->Bs.vec[col];
1408 col = sys->J.reg.col.low;
1409 for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
1410 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1411 if(set_is_member(sys->J.varpivots,ind) ) {
1412 sys->ZBZ.Zy[nz.row] += update->y.vec[col]*
1413 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1414 sys->ZBZ.ZBs[nz.row] += update->Bs.vec[col]*
1415 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1416 }
1417 }
1418 for( nz.col=0; nz.col <= nz.row; nz.col++ ) {
1419 sys->ZBZ.mtx[nz.row][nz.col] += update->ys > 0.0 ?
1420 sys->ZBZ.Zy[nz.row]*sys->ZBZ.Zy[nz.col]/update->ys : 0.0;
1421 sys->ZBZ.mtx[nz.row][nz.col] -= update->sBs > 0.0 ?
1422 sys->ZBZ.ZBs[nz.row]*sys->ZBZ.ZBs[nz.col]/update->sBs : 0.0;
1423 if(nz.row != nz.col ) {
1424 sys->ZBZ.mtx[nz.col][nz.row] =
1425 sys->ZBZ.mtx[nz.row][nz.col];
1426 }
1427 }
1428 }
1429 }
1430 }
1431 sys->ZBZ.accurate = TRUE;
1432 #if DEBUG
1433 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nReduced Hessian: \n");
1434 debug_out_hessian(LIF(sys),sys);
1435 #endif
1436 }
1437
1438
1439 /**
1440 Calculates just the jacobian RHS. This function should be used to
1441 supplement calculation of the jacobian. The vector vec must
1442 already be calculated and scaled so as to simply be added to the
1443 rhs. Caller is responsible for initially zeroing the rhs vector.
1444 */
1445 static void calc_rhs(qrslv_system_t sys, struct vec_vector *vec,
1446 real64 scalar, boolean transpose
1447 ){
1448 if(transpose ) { /* vec is indexed by col */
1449 int32 col;
1450 for( col=vec->rng->low; col<=vec->rng->high; col++ ) {
1451 sys->J.rhs[mtx_col_to_org(sys->J.mtx,col)] += scalar*vec->vec[col];
1452 }
1453 }else{ /* vec is indexed by row */
1454 int32 row;
1455 for( row=vec->rng->low; row<=vec->rng->high; row++ ) {
1456 sys->J.rhs[mtx_row_to_org(sys->J.mtx,row)] += scalar*vec->vec[row];
1457 }
1458 }
1459 linsolqr_rhs_was_changed(sys->J.sys,sys->J.rhs);
1460 }
1461
1462
1463 /**
1464 Computes the lagrange multipliers for the equality constraints.
1465 */
1466 static void calc_multipliers(qrslv_system_t sys){
1467
1468 if(sys->multipliers.accurate)return;
1469
1470 if(!OPTIMIZING(sys)){
1471 zero_vector(&(sys->multipliers));
1472 sys->multipliers.norm2 = 0.0;
1473 }else{
1474 linsolqr_system_t lsys = sys->J.sys;
1475 int32 row;
1476 sys->J.rhs = linsolqr_get_rhs(lsys,0);
1477 mtx_zero_real64(sys->J.rhs,sys->cap);
1478 calc_rhs(sys, &(sys->gradient), -1.0, TRUE );
1479 linsolqr_solve(lsys,sys->J.rhs);
1480 row = sys->multipliers.rng->low;
1481 for( ; row <= sys->multipliers.rng->high; row++ ) {
1482 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1483 sys->multipliers.vec[row] = linsolqr_var_value(
1484 lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row)
1485 );
1486 rel_set_multiplier(rel,sys->multipliers.vec[row]*
1487 sys->weights.vec[row]);
1488
1489 }
1490 if(SLV_PARAM_BOOL(&(sys->p),SAVLIN)) {
1491 FILE *ldat;
1492 int32 ov;
1493 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1494 ldat=fopen(savlinfilename,"w");
1495 FPRINTF(ldat,
1496 "================= multipliersrhs (orgcoled) itn %d =====\n",
1497 sys->s.iteration);
1498 debug_write_array(ldat,sys->J.rhs,sys->cap);
1499 FPRINTF(ldat,
1500 "================= multipliers (orgrowed) ============\n");
1501 for(ov=0 ; ov < sys->cap; ov++ )
1502 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1503 fclose(ldat);
1504 }
1505 square_norm( &(sys->multipliers) );
1506 }
1507 sys->multipliers.accurate = TRUE;
1508 #if DEBUG
1509 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Multipliers: ");
1510 debug_out_vector(LIF(sys),sys,&(sys->multipliers));
1511 #endif
1512 }
1513
1514
1515 /**
1516 Computes the gradient of the lagrangian which
1517 should be zero at the optimum solution.
1518 */
1519 static void calc_stationary( qrslv_system_t sys){
1520 if(sys->stationary.accurate )
1521 return;
1522
1523 if(!OPTIMIZING(sys)){
1524 zero_vector(&(sys->stationary));
1525 sys->stationary.norm2 = 0.0;
1526 }else{
1527 int32 col;
1528 matrix_product(sys->J.mtx, &(sys->multipliers),
1529 &(sys->stationary), 1.0, TRUE);
1530 col = sys->stationary.rng->low;
1531 for( ; col <= sys->stationary.rng->high; col++ )
1532 sys->stationary.vec[col] += sys->gradient.vec[col];
1533 square_norm( &(sys->stationary) );
1534 }
1535 sys->stationary.accurate = TRUE;
1536 #if DEBUG
1537 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Stationary: ");
1538 debug_out_vector(LIF(sys),sys,&(sys->stationary));
1539 #endif
1540 }
1541
1542
1543 /**
1544 Calculate the gamma vector.
1545 */
1546 static void calc_gamma( qrslv_system_t sys){
1547 if(sys->gamma.accurate)return;
1548
1549 matrix_product(sys->J.mtx, &(sys->residuals),
1550 &(sys->gamma), -1.0, TRUE);
1551 square_norm( &(sys->gamma) );
1552 sys->gamma.accurate = TRUE;
1553 #if DEBUG
1554 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Gamma: ");
1555 debug_out_vector(LIF(sys),sys,&(sys->gamma));
1556 #endif
1557 }
1558
1559 /**
1560 Calculate the Jgamma vector.
1561 */
1562 static void calc_Jgamma( qrslv_system_t sys){
1563 if(sys->Jgamma.accurate)return;
1564
1565 matrix_product(sys->J.mtx, &(sys->gamma),
1566 &(sys->Jgamma), 1.0, FALSE);
1567 square_norm( &(sys->Jgamma) );
1568 sys->Jgamma.accurate = TRUE;
1569 #if DEBUG
1570 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Jgamma: ");
1571 debug_out_vector(LIF(sys),sys,&(sys->Jgamma));
1572 #endif
1573 }
1574
1575
1576 /**
1577 Computes a step to solve the linearized equations.
1578 */
1579 static void calc_newton( qrslv_system_t sys){
1580 linsolqr_system_t lsys = sys->J.sys;
1581 int32 col;
1582
1583 if(sys->newton.accurate)return;
1584
1585 sys->J.rhs = linsolqr_get_rhs(lsys,1);
1586 mtx_zero_real64(sys->J.rhs,sys->cap);
1587 calc_rhs(sys, &(sys->residuals), -1.0, FALSE);
1588 linsolqr_solve(lsys,sys->J.rhs);
1589 col = sys->newton.rng->low;
1590 for( ; col <= sys->newton.rng->high; col++ ) {
1591 sys->newton.vec[col] =
1592 linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1593 }
1594 if(SLV_PARAM_BOOL(&(sys->p),SAVLIN)) {
1595 FILE *ldat;
1596 int32 ov;
1597 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1598 ldat=fopen(savlinfilename,"w");
1599 FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1600 sys->s.iteration);
1601 debug_write_array(ldat,sys->J.rhs,sys->cap);
1602 FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1603 for(ov=0 ; ov < sys->cap; ov++ )
1604 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1605 fclose(ldat);
1606 }
1607 square_norm( &(sys->newton) );
1608 sys->newton.accurate = TRUE;
1609 #if DEBUG
1610 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Newton: ");
1611 debug_out_vector(LIF(sys),sys,&(sys->newton));
1612 #endif
1613 }
1614
1615
1616 /**
1617 Computes an update to the product B and newton.
1618 */
1619 static void calc_Bnewton( qrslv_system_t sys){
1620 if(sys->Bnewton.accurate)return;
1621
1622 if(!OPTIMIZING(sys)){
1623 zero_vector(&(sys->Bnewton));
1624 sys->Bnewton.norm2 = 0.0;
1625 }else{
1626 struct hessian_data *update;
1627 copy_vector(&(sys->newton),&(sys->Bnewton));
1628 for( update=sys->B; update != NULL; update = update->next ) {
1629 int32 col;
1630 real64 Yn = inner_product( &(update->y),&(sys->newton) );
1631 real64 sBn = inner_product( &(update->Bs),&(sys->newton) );
1632 col = sys->Bnewton.rng->low;
1633 for( ; col <= sys->Bnewton.rng->high; col++ ) {
1634 sys->Bnewton.vec[col] += update->ys > 0.0 ?
1635 (update->y.vec[col])*Yn/update->ys : 0.0;
1636 sys->Bnewton.vec[col] -= update->sBs > 0.0 ?
1637 (update->Bs.vec[col])*sBn/update->sBs : 0.0;
1638 }
1639 }
1640 square_norm( &(sys->Bnewton) );
1641 }
1642 sys->Bnewton.accurate = TRUE;
1643 #if DEBUG
1644 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Bnewton: ");
1645 debug_out_vector(LIF(sys),sys,&(sys->Bnewton));
1646 #endif
1647 }
1648
1649
1650 /**
1651 Calculate the nullspace move if OPTIMIZING.
1652 */
1653 static void calc_nullspace( qrslv_system_t sys){
1654 if(sys->nullspace.accurate)return;
1655
1656 if(!OPTIMIZING(sys)){
1657 zero_vector(&(sys->nullspace));
1658 sys->nullspace.norm2 = 0.0;
1659 }else{
1660 mtx_coord_t nz;
1661 zero_vector(&(sys->nullspace));
1662 for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1663 int32 col, dep, ndx;
1664 col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
1665 dep = mtx_col_to_org(sys->J.mtx,col);
1666 sys->nullspace.vec[col] = -sys->stationary.vec[col] -
1667 sys->Bnewton.vec[col];
1668 ndx = sys->J.reg.col.low;
1669 for( ; ndx <= sys->J.reg.col.high - sys->ZBZ.order; ndx++ ) {
1670 int32 ind = mtx_col_to_org(sys->J.mtx,ndx);
1671 if(set_is_member(sys->J.varpivots,ind)){
1672 sys->nullspace.vec[col] -=
1673 (sys->stationary.vec[ndx] + sys->Bnewton.vec[ndx])*
1674 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1675 }
1676 }
1677 }
1678 /*
1679 * Must invert ZBZ first. It's symmetric so
1680 * can find Cholesky factors. Essentially, find
1681 * the "square root" of the matrix such that
1682 *
1683 * T
1684 * L U = U U = ZBZ, where U is an upper triangular
1685 * matrix.
1686 */
1687 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1688 for( nz.col = nz.row; nz.col < sys->ZBZ.order; nz.col++ ) {
1689 int32 col;
1690 for( col = 0; col < nz.row; col++ )
1691 sys->ZBZ.mtx[nz.row][nz.col] -=
1692 sys->ZBZ.mtx[nz.row][col]*
1693 sys->ZBZ.mtx[col][nz.col];
1694 if(nz.row == nz.col )
1695 sys->ZBZ.mtx[nz.row][nz.col] =
1696 calc_sqrt_D0(sys->ZBZ.mtx[nz.row][nz.col]);
1697 else {
1698 sys->ZBZ.mtx[nz.row][nz.col] /=
1699 sys->ZBZ.mtx[nz.row][nz.row];
1700 sys->ZBZ.mtx[nz.col][nz.row] =
1701 sys->ZBZ.mtx[nz.row][nz.col];
1702 }
1703 }
1704 }
1705 #if DEBUG
1706 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nInverse Reduced Hessian: \n");
1707 debug_out_hessian(LIF(sys),sys);
1708 #endif
1709 /*
1710 * forward substitute
1711 */
1712 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1713 int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1714 for( nz.col = 0; nz.col < nz.row; nz.col++ ) {
1715 sys->nullspace.vec[nz.row+offset] -=
1716 sys->nullspace.vec[nz.col+offset]*
1717 sys->ZBZ.mtx[nz.row][nz.col];
1718 }
1719 sys->nullspace.vec[nz.row+offset] /=
1720 sys->ZBZ.mtx[nz.row][nz.row];
1721 }
1722
1723 /*
1724 * backward substitute
1725 */
1726 for( nz.row = sys->ZBZ.order-1; nz.row >= 0; nz.row-- ) {
1727 int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1728 for( nz.col = nz.row+1; nz.col < sys->ZBZ.order; nz.col++ ) {
1729 sys->nullspace.vec[nz.row+offset] -=
1730 sys->nullspace.vec[nz.col+offset]*
1731 sys->ZBZ.mtx[nz.row][nz.col];
1732 }
1733 sys->nullspace.vec[nz.row+offset] /=
1734 sys->ZBZ.mtx[nz.row][nz.row];
1735 }
1736 square_norm( &(sys->nullspace) );
1737 }
1738 sys->nullspace.accurate = TRUE;
1739 #if DEBUG
1740 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Nullspace: ");
1741 debug_out_vector(LIF(sys),sys,&(sys->nullspace));
1742 #endif
1743 }
1744
1745 /**
1746 Calculate the 1st order descent direction for phi
1747 in the variables.
1748 */
1749 static void calc_varstep1( qrslv_system_t sys){
1750 if(sys->varstep1.accurate )
1751 return;
1752
1753 if(!OPTIMIZING(sys)){
1754 copy_vector(&(sys->gamma),&(sys->varstep1));
1755 sys->varstep1.norm2 = sys->gamma.norm2;
1756 }else{
1757 int32 col;
1758 col = sys->varstep1.rng->low;
1759 for( ; col <= sys->varstep1.rng->high; col++ )
1760 sys->varstep1.vec[col] = SLV_PARAM_REAL(&(sys->p),RHO)*sys->gamma.vec[col] -
1761 sys->stationary.vec[col];
1762 square_norm( &(sys->varstep1) );
1763 }
1764 sys->varstep1.accurate = TRUE;
1765 #if DEBUG
1766 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Varstep1: ");
1767 debug_out_vector(LIF(sys),sys,&(sys->varstep1));
1768 #endif
1769 }
1770
1771
1772 /**
1773 Computes an update to the product B and varstep1.
1774 */
1775 static void calc_Bvarstep1( qrslv_system_t sys){
1776 if(sys->Bvarstep1.accurate )
1777 return;
1778
1779 if(!OPTIMIZING(sys)){
1780 zero_vector(&(sys->Bvarstep1));
1781 sys->Bvarstep1.norm2 = 0.0;
1782 }else{
1783 struct hessian_data *update;
1784 copy_vector(&(sys->varstep1),&(sys->Bvarstep1));
1785 for( update=sys->B; update != NULL; update = update->next ) {
1786 int32 col;
1787 real64 yv = inner_product( &(update->y),&(sys->varstep1) );
1788 real64 sBv = inner_product( &(update->Bs),&(sys->varstep1) );
1789 col = sys->Bvarstep1.rng->low;
1790 for( ; col <= sys->Bvarstep1.rng->high; col++ ) {
1791 sys->Bvarstep1.vec[col] += update->ys > 0.0 ?
1792 (update->y.vec[col])*yv/update->ys : 0.0;
1793 sys->Bvarstep1.vec[col] -= update->sBs > 0.0 ?
1794 (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1795 }
1796 }
1797 square_norm( &(sys->Bvarstep1) );
1798 }
1799 sys->Bvarstep1.accurate = TRUE;
1800 #if DEBUG
1801 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Bvarstep1: ");
1802 debug_out_vector(LIF(sys),sys,&(sys->Bvarstep1));
1803 #endif
1804 }
1805
1806
1807 /**
1808 Calculate the 2nd order descent direction for phi
1809 in the variables.
1810 */
1811 static void calc_varstep2( qrslv_system_t sys){
1812 if(sys->varstep2.accurate )
1813 return;
1814
1815 if(!OPTIMIZING(sys)){
1816 copy_vector(&(sys->newton),&(sys->varstep2));
1817 sys->varstep2.norm2 = sys->newton.norm2;
1818 }else{
1819 int32 col;
1820 col = sys->varstep2.rng->low;
1821 for( ; col <= sys->varstep2.rng->high - sys->ZBZ.order ; ++col ) {
1822 int32 dep;
1823 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1824 sys->varstep2.vec[col] = sys->newton.vec[col];
1825 if(set_is_member(sys->J.varpivots,ind) ) {
1826 dep = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1827 for( ; dep <= sys->varstep2.rng->high; dep++ )
1828 sys->varstep2.vec[col] += sys->nullspace.vec[dep]*
1829 (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1830 }
1831 }
1832 col = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1833 for( ; col <= sys->varstep2.rng->high; ++col )
1834 sys->varstep2.vec[col] = sys->nullspace.vec[col] +
1835 sys->newton.vec[col];
1836 square_norm( &(sys->varstep2) );
1837 }
1838 sys->varstep2.accurate = TRUE;
1839 #if DEBUG
1840 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Varstep2: ");
1841 debug_out_vector(LIF(sys),sys,&(sys->varstep2));
1842 #endif
1843 }
1844
1845
1846 /**
1847 Computes an update to the product B and varstep2.
1848 */
1849 static void calc_Bvarstep2( qrslv_system_t sys){
1850 if(sys->Bvarstep2.accurate )
1851 return;
1852
1853 if(!OPTIMIZING(sys)){
1854 zero_vector(&(sys->Bvarstep2));
1855 sys->Bvarstep2.norm2 = 0.0;
1856 }else{
1857 struct hessian_data *update;
1858 copy_vector(&(sys->varstep2),&(sys->Bvarstep2));
1859 for( update=sys->B; update != NULL; update = update->next ) {
1860 int32 col;
1861 real64 yv = inner_product( &(update->y),&(sys->varstep2) );
1862 real64 sBv = inner_product( &(update->Bs),&(sys->varstep2) );
1863 col = sys->Bvarstep2.rng->low;
1864 for( ; col <= sys->Bvarstep2.rng->high; col++ ) {
1865 sys->Bvarstep2.vec[col] += update->ys > 0.0 ?
1866 (update->y.vec[col])*yv/update->ys : 0.0;
1867 sys->Bvarstep2.vec[col] -= update->sBs > 0.0 ?
1868 (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1869 }
1870 }
1871 square_norm( &(sys->Bvarstep2) );
1872 }
1873 sys->Bvarstep2.accurate = TRUE;
1874 #if DEBUG
1875 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Bvarstep2: ");
1876 debug_out_vector(LIF(sys),sys,&(sys->Bvarstep2));
1877 #endif
1878 }
1879
1880
1881 /**
1882 Calculate the negative gradient direction of phi in the
1883 multipliers.
1884 */
1885 static void calc_mulstep1( qrslv_system_t sys){
1886 if(sys->mulstep1.accurate )
1887 return;
1888
1889 if(!OPTIMIZING(sys)){
1890 zero_vector(&(sys->mulstep1));
1891 sys->mulstep1.norm2 = 0.0;
1892 }else{
1893 int32 row;
1894 row = sys->mulstep1.rng->low;
1895 for( ; row <= sys->mulstep1.rng->high; row++ )
1896 sys->mulstep1.vec[row] = -sys->residuals.vec[row];
1897 square_norm( &(sys->mulstep1) );
1898 }
1899 sys->mulstep1.accurate = TRUE;
1900 #if DEBUG
1901 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Mulstep1: ");
1902 debug_out_vector(LIF(sys),sys,&(sys->mulstep1));
1903 #endif
1904 }
1905
1906
1907 /**
1908 Calculate the mulstep2 direction of phi in the
1909 multipliers.
1910 */
1911 static void calc_mulstep2( qrslv_system_t sys){
1912 if(sys->mulstep2.accurate )
1913 return;
1914
1915 if(!OPTIMIZING(sys)){
1916 zero_vector(&(sys->mulstep2));
1917 sys->mulstep2.norm2 = 0.0;
1918 }else{
1919 linsolqr_system_t lsys = sys->J.sys;
1920 int32 row;
1921 sys->J.rhs = linsolqr_get_rhs(lsys,2);
1922 mtx_zero_real64(sys->J.rhs,sys->cap);
1923 calc_rhs(sys, &(sys->Bvarstep2), -1.0, TRUE);
1924 calc_rhs(sys, &(sys->stationary), -1.0, TRUE);
1925 linsolqr_solve(lsys,sys->J.rhs);
1926 row = sys->mulstep2.rng->low;
1927 for( ; row <= sys->mulstep2.rng->high; row++ )
1928 sys->mulstep2.vec[row] = linsolqr_var_value
1929 (lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row));
1930 if(SLV_PARAM_BOOL(&(sys->p),SAVLIN)) {
1931 FILE *ldat;
1932 int32 ov;
1933 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1934 ldat=fopen(savlinfilename,"w");
1935 FPRINTF(ldat,
1936 "================= mulstep2rhs (orgcoled) itn %d =======\n",
1937 sys->s.iteration);
1938 debug_write_array(ldat,sys->J.rhs,sys->cap);
1939 FPRINTF(ldat,
1940 "================= mulstep2vars (orgrowed) ============\n");
1941 for(ov=0 ; ov < sys->cap; ov++ )
1942 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1943 fclose(ldat);
1944 }
1945 square_norm( &(sys->mulstep2) );
1946 }
1947 sys->mulstep2.accurate = TRUE;
1948 #if DEBUG
1949 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Mulstep2: ");
1950 debug_out_vector(LIF(sys),sys,&(sys->mulstep2));
1951 #endif
1952 }
1953
1954
1955 /**
1956 Computes the global minimizing function Phi.
1957 */
1958 static void calc_phi( qrslv_system_t sys){
1959 if(!OPTIMIZING(sys)){
1960 sys->phi = 0.5*sys->residuals.norm2;
1961 }else{
1962 sys->phi = sys->objective;
1963 sys->phi += inner_product( &(sys->multipliers),&(sys->residuals) );
1964 sys->phi += 0.5*SLV_PARAM_REAL(&(sys->p),RHO)*sys->residuals.norm2;
1965 }
1966 }
1967
1968 /*------------------------------------------------------------------------------
1969 STEP CALCULATION STUFF
1970
1971 * OK. Here's where we compute the actual step to be taken. It will
1972 * be some linear combination of the 1st order and 2nd order steps.
1973 */
1974
1975 typedef real64 sym_2x2_t[3]; /* Stores symmetric 2x2 matrices */
1976
1977 struct parms_t {
1978 real64 low,high,guess; /* Used to search for parameter */
1979 };
1980
1981 struct calc_step_vars {
1982 sym_2x2_t coef1, coef2;
1983 real64 rhs[2]; /* RHS for 2x2 system */
1984 struct parms_t parms;
1985 real64 alpha1, alpha2;
1986 real64 error; /* Error between step norm and sys->maxstep */
1987 };
1988
1989 /**
1990 Calculates 2x2 system (coef1,coef2,rhs).
1991 */
1992 static void calc_2x2_system(qrslv_system_t sys, struct calc_step_vars *vars){
1993 vars->coef1[0] = (2.0*sys->phi/sys->newton.norm2)*
1994 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1995 vars->coef1[1] = 1.0;
1996 vars->coef1[2] = (sys->Jgamma.norm2/sys->gamma.norm2)*
1997 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1998
1999 vars->coef2[0] = 1.0;
2000 vars->coef2[1] = 2.0*sys->phi/
2001 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
2002 vars->coef2[2] = 1.0;
2003
2004 vars->rhs[0] = 2.0*sys->phi/
2005 sys->maxstep/calc_sqrt_D0(sys->gamma.norm2);
2006 vars->rhs[1] = calc_sqrt_D0(sys->newton.norm2)/sys->maxstep;
2007 }
2008
2009 /**
2010 Determines alpha1 and alpha2 from the parameter (guess).
2011 */
2012 static void coefs_from_parm( qrslv_system_t sys, struct calc_step_vars *vars){
2013
2014 sym_2x2_t coef; /* Actual coefficient matrix */
2015 real64 det; /* Determinant of coefficient matrix */
2016 int i;
2017
2018 for(i=0; i<3; ++i) coef[i]= vars->coef1[i] + vars->parms.guess * vars->coef2[i];
2019
2020 det = coef[0]*coef[2] - coef[1]*coef[1];
2021 if(det < 0.0){
2022 ERROR_REPORTER_HERE(ASC_PROG_ERROR,"Unexpected negative determinant %f.", det);
2023 }
2024
2025 if(det <= SLV_PARAM_REAL(&(sys->p),DETZERO) ) {
2026 /*
2027 varstep2 and varstep1 are essentially parallel:
2028 adjust length of either
2029 */
2030 vars->alpha2 = 0.0;
2031 vars->alpha1 = 1.0;
2032 }else{
2033 vars->alpha2 = (vars->rhs[0]*coef[2] - vars->rhs[1]*coef[1])/det;
2034 vars->alpha1 = (vars->rhs[1]*coef[0] - vars->rhs[0]*coef[1])/det;
2035 }
2036 }
2037
2038 /**
2039 Computes step vector length based on 1st order and 2nd order
2040 vectors and their coefficients.
2041 */
2042 static real64 step_norm2( qrslv_system_t sys, struct calc_step_vars *vars){
2043 return sys->maxstep*sys->maxstep*
2044 (vars->alpha2 * vars->alpha2 +
2045 vars->alpha2 * vars->alpha1 * sys->phi/
2046 calc_sqrt_D0(sys->varstep2.norm2 + sys->mulstep2.norm2)/
2047 calc_sqrt_D0(sys->varstep1.norm2 + sys->mulstep1.norm2) +
2048 vars->alpha1 * vars->alpha1);
2049 }
2050
2051 /**
2052 Re-guesses the parameters based on step size vs. target value.
2053 */
2054 static void adjust_parms( qrslv_system_t sys, struct calc_step_vars *vars){
2055 vars->error = (calc_sqrt_D0(step_norm2(sys,vars))/sys->maxstep) - 1.0;
2056 if(vars->error > 0.0 ) {
2057 /* Increase parameter (to decrease step length) */
2058 vars->parms.low = vars->parms.guess;
2059 vars->parms.guess = (vars->parms.high>3.0*vars->parms.guess)
2060 ? 2.0*vars->parms.guess
2061 : 0.5*(vars->parms.low + vars->parms.high);
2062 }else{
2063 /* Decrease parameter (to increase step norm) */
2064 vars->parms.high = vars->parms.guess;
2065 vars->parms.guess = 0.5*(vars->parms.low + vars->parms.high);
2066 }
2067 }
2068
2069 /**
2070 Computes the step based on the coefficients in vars.
2071 */
2072 static void compute_step( qrslv_system_t sys, struct calc_step_vars *vars){
2073 int32 row,col;
2074 real64 tot1_norm2, tot2_norm2;
2075
2076 tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
2077 tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
2078 if(!sys->varstep.accurate ) {
2079 for( col=sys->varstep.rng->low ; col<=sys->varstep.rng->high ; ++col )
2080 if((vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
2081 sys->varstep.vec[col] = sys->maxstep *
2082 sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2);
2083 }else if((vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
2084 sys->varstep.vec[col] = sys->maxstep *
2085 sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2);
2086 }else if((vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
2087 sys->varstep.vec[col] = sys->maxstep*
2088 (
2089 vars->alpha2*sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2) +
2090 vars->alpha1*sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2)
2091 );
2092 }
2093 sys->varstep.accurate = TRUE;
2094 }
2095 if(!sys->mulstep.accurate ) {
2096 for( row=sys->mulstep.rng->low ; row<=sys->mulstep.rng->high ; ++row )
2097 if((vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
2098 sys->mulstep.vec[row] = sys->maxstep *
2099 sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2);
2100 }else if((vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
2101 sys->mulstep.vec[row] = sys->maxstep *
2102 sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2);
2103 }else if((vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
2104 sys->mulstep.vec[row] = sys->maxstep*
2105 (
2106 vars->alpha2*sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2) +
2107 vars->alpha1*sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2)
2108 );
2109 }
2110 sys->mulstep.accurate = TRUE;
2111 }
2112 #if DEBUG
2113 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Varstep: ");
2114 debug_out_vector(LIF(sys),sys,&(sys->varstep));
2115 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Mulstep: ");
2116 debug_out_vector(LIF(sys),sys,&(sys->mulstep));
2117 #endif
2118 }
2119
2120
2121 /**
2122 Calculates step vector, based on sys->maxstep, and the varstep2/
2123 varstep1 and mulstep2/mulstep1 vectors. Nothing is assumed to be
2124 calculated, except the weights and the jacobian (scaled). Also,
2125 the step is not checked for legitimacy.
2126 NOTE: the step is scaled.
2127 */
2128 static void calc_step( qrslv_system_t sys, int minor){
2129
2130 struct calc_step_vars vars;
2131 real64 tot1_norm2, tot2_norm2;
2132
2133 if(sys->varstep.accurate && sys->mulstep.accurate )
2134 return;
2135 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2136 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\n%-40s ---> %d\n", " Step trial",minor);
2137 }
2138
2139 tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
2140 tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
2141 if((tot1_norm2 == 0.0) && (tot2_norm2 == 0.0) ) {
2142 /* Take no step at all */
2143 vars.alpha1 = 0.0;
2144 vars.alpha2 = 0.0;
2145 sys->maxstep = 0.0;
2146 sys->varstep.norm2 = 0.0;
2147 sys->mulstep.norm2 = 0.0;
2148
2149 }else if(tot2_norm2 > 0.0 && OPTIMIZING(sys)){
2150 /* Stay in varstep2 direction */
2151 vars.alpha1 = 0.0;
2152 vars.alpha2 = 1.0;
2153 sys->maxstep = MIN(sys->maxstep,calc_sqrt_D0(tot2_norm2));
2154 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2155 sys->varstep2.norm2/tot2_norm2;
2156 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2157 sys->mulstep2.norm2/tot2_norm2;
2158
2159 }else if((tot2_norm2>0.0)&&(calc_sqrt_D0(tot2_norm2)<=sys->maxstep) ) {
2160 /* Attempt step in varstep2 direction */
2161 vars.alpha1 = 0.0;
2162 vars.alpha2 = 1.0;
2163 sys->maxstep = calc_sqrt_D0(tot2_norm2);
2164 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2165 sys->varstep2.norm2/tot2_norm2;
2166 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2167 sys->mulstep2.norm2/tot2_norm2;
2168
2169 }else if((tot2_norm2==0.0 || sys->s.block.current_size==1) &&
2170 (tot1_norm2 > 0.0) ) {
2171 /* Attempt step in varstep1 direction */
2172 vars.alpha1 = 1.0;
2173 vars.alpha2 = 0.0;
2174 if( (sys->gamma.norm2/sys->Jgamma.norm2)*
2175 calc_sqrt_D0(sys->gamma.norm2) <= sys->maxstep )
2176 sys->maxstep = (sys->gamma.norm2/sys->Jgamma.norm2)*
2177 calc_sqrt_D0(sys->gamma.norm2);
2178 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2179 sys->varstep1.norm2/tot1_norm2;
2180 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2181 sys->mulstep1.norm2/tot1_norm2;
2182
2183 }else{
2184 /* Attempt step in varstep1-varstep2 direction */
2185 vars.parms.low = 0.0;
2186 vars.parms.high = MAXDOUBLE;
2187 vars.parms.guess = 1.0;
2188 calc_2x2_system(sys,&vars);
2189 do {
2190 coefs_from_parm(sys, &vars);
2191 adjust_parms(sys, &vars);
2192 } while( fabs(vars.error) > SLV_PARAM_REAL(&(sys->p),STEPSIZEERR_MAX) &&
2193 vars.parms.high - vars.parms.low > SLV_PARAM_REAL(&(sys->p),PARMRNG_MIN) );
2194 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2195 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2196 " parameter high", vars.parms.high);
2197 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2198 " parameter low", vars.parms.low);
2199 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2200 " Error in step length", vars.error);
2201 }
2202 sys->varstep.norm2 = step_norm2(sys, &vars);
2203 sys->mulstep.norm2 = 0.0;
2204 }
2205
2206 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2207 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", " Alpha1 coefficient (normalized)",
2208 vars.alpha1);
2209 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", " Alpha2 coefficient (normalized)",
2210 vars.alpha2);
2211 }
2212 compute_step(sys,&vars);
2213 return;
2214
2215 }
2216
2217 /*------------------------------------------------------------------------------
2218 VARIABLE VALUES MAINTENANCE
2219 */
2220
2221 /**
2222 Restores the values of the variables before applying
2223 a step.
2224 */
2225 static void restore_variables( qrslv_system_t sys){
2226 int32 col;
2227 real64 *vec;
2228 vec = (sys->nominals.vec);
2229 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
2230 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2231 var_set_value(var,sys->variables.vec[col]*vec[col]);
2232 }
2233 }
2234
2235
2236 /**
2237 Calculates the maximum fraction of the step which can be
2238 taken without going out of bounds. If the entire step can be
2239 taken, 1.0 is returned. Otherwise a value less than 1 is
2240 returned. It is assumed that the current variable values
2241 are within their bounds. The step must be calculated.
2242 */
2243 static real64 required_coef_to_stay_inbounds( qrslv_system_t sys){
2244 real64 mincoef;
2245 int32 col;
2246 real64 *vec;
2247 vec = (sys->nominals.vec);
2248
2249 if(sys->p.ignore_bounds )
2250 return(1.0);
2251
2252 mincoef = 1.0;
2253 for( col=sys->varstep.rng->low; col <= sys->varstep.rng->high; col++ ) {
2254 struct var_variable *var;
2255 real64 coef,dx,val,bnd;
2256 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2257 coef = 1.0;
2258 dx = sys->varstep.vec[col] * vec[col];
2259 bnd = var_upper_bound(var);
2260 if((val=var_value(var)) + dx > bnd )
2261 coef = MIN((bnd-val)/dx, 1.0);
2262 bnd = var_lower_bound(var);
2263 if(val + dx < bnd )
2264 coef = MIN((bnd-val)/dx, 1.0);
2265 if(coef < mincoef )
2266 mincoef = coef;
2267 }
2268 return(mincoef);
2269 }
2270
2271
2272 /**
2273 Adds sys->varstep to the variable values in block: projecting
2274 near bounds.
2275 */
2276 static void apply_step( qrslv_system_t sys){
2277 FILE *lif = LIF(sys);
2278 int nproj = 0;
2279 real64 bounds_coef = 1.0;
2280 int32 col;
2281 real64 *vec;
2282 vec = (sys->nominals.vec);
2283
2284 if(SLV_PARAM_BOOL(&(sys->p),TRUNCATE) && (!sys->p.ignore_bounds))
2285 bounds_coef = required_coef_to_stay_inbounds(sys);
2286
2287 for( col=sys->varstep.rng->low; col <= sys->varstep.rng->high; col++ ) {
2288 struct var_variable *var;
2289 real64 dx,val,bnd;
2290 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2291 dx = vec[col]*sys->varstep.vec[col];
2292 val = var_value(var);
2293 if(bounds_coef < 1.0) {
2294 dx = dx*SLV_PARAM_REAL(&(sys->p),TOWARD_BOUNDS)*bounds_coef;
2295 sys->varstep.vec[col] = dx/vec[col];
2296 }else{
2297 if(!sys->p.ignore_bounds ) {
2298 if(val + dx > (bnd=var_upper_bound(var)) ) {
2299 dx = SLV_PARAM_REAL(&(sys->p),TOWARD_BOUNDS)*(bnd-val);
2300 sys->varstep.vec[col] = dx/vec[col];
2301 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2302 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> ",
2303 " Variable projected to upper bound");
2304 print_var_name(lif,sys,var); PUTC('\n',lif);
2305 }
2306 ++nproj;
2307 }else if(val + dx < (bnd=var_lower_bound(var)) ) {
2308 dx = SLV_PARAM_REAL(&(sys->p),TOWARD_BOUNDS)*(bnd-val);
2309 sys->varstep.vec[col] = dx/vec[col];
2310 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2311 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> ",
2312 " Variable projected to lower bound");
2313 print_var_name(lif,sys,var); PUTC('\n',lif);
2314 }
2315 ++nproj;
2316 }
2317 }
2318 }
2319 var_set_value(var,val+dx);
2320 }
2321
2322 if(!sys->p.ignore_bounds ) {
2323 if(nproj > 0) {
2324 square_norm(&(sys->varstep));
2325 sys->progress = calc_sqrt_D0
2326 (calc_sqrt_D0((sys->varstep.norm2 + sys->mulstep.norm2)*
2327 (sys->varstep1.norm2 + sys->mulstep1.norm2)));
2328 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2329 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", " Projected step length (scaled)",
2330 calc_sqrt_D0(sys->varstep.norm2 + sys->mulstep.norm2));
2331 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2332 " Projected progress", sys->progress);
2333 }
2334 }
2335 if(bounds_coef < 1.0) {
2336 square_norm(&(sys->varstep));
2337 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2338 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2339 " Truncated step length (scaled)",
2340 calc_sqrt_D0(sys->varstep.norm2 + sys->mulstep.norm2));
2341 }
2342 sys->progress = calc_sqrt_D0
2343 (calc_sqrt_D0((sys->varstep.norm2 + sys->mulstep.norm2)*
2344 (sys->varstep1.norm2 + sys->mulstep1.norm2)));
2345 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2346 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2347 " Truncated progress", sys->progress);
2348 }
2349 }
2350 }
2351
2352 /* Allow weighted residuals to be recalculated at new point */
2353 sys->residuals.accurate = FALSE;
2354
2355 return;
2356 }
2357
2358 /**
2359 This function should be called when the step is accepted.
2360 */
2361 static void step_accepted( qrslv_system_t sys){
2362 /* Maintain update status on jacobian and weights */
2363 if(--(sys->update.jacobian) <= 0)
2364 sys->J.accurate = FALSE;
2365
2366 sys->ZBZ.accurate = FALSE;
2367 sys->variables.accurate = FALSE;
2368 sys->gradient.accurate = FALSE;
2369 sys->multipliers.accurate = FALSE;
2370 sys->stationary.accurate = FALSE;
2371 sys->newton.accurate = FALSE;
2372 sys->Bnewton.accurate = FALSE;
2373 sys->nullspace.accurate = FALSE;
2374 sys->gamma.accurate = FALSE;
2375 sys->Jgamma.accurate = FALSE;
2376 sys->varstep1.accurate = FALSE;
2377 sys->Bvarstep1.accurate = FALSE;
2378 sys->varstep2.accurate = FALSE;
2379 sys->Bvarstep2.accurate = FALSE;
2380 sys->mulstep1.accurate = FALSE;
2381 sys->mulstep2.accurate = FALSE;
2382 sys->varstep.accurate = FALSE;
2383 sys->mulstep.accurate = FALSE;
2384
2385 if(!OPTIMIZING(sys)){
2386 sys->ZBZ.accurate = TRUE;
2387 sys->gradient.accurate = TRUE;
2388 sys->multipliers.accurate = TRUE;
2389 sys->stationary.accurate = TRUE;
2390 sys->Bnewton.accurate = TRUE;
2391 sys->nullspace.accurate = TRUE;
2392 sys->Bvarstep1.accurate = TRUE;
2393 sys->Bvarstep2.accurate = TRUE;
2394 }
2395 }
2396
2397 /**
2398 This function changes sys->maxstep to the given number and should be
2399 called whenever sys->maxstep is to be changed.
2400 */
2401 static void change_maxstep( qrslv_system_t sys, real64 maxstep){
2402 sys->maxstep = maxstep;
2403 sys->varstep.accurate = FALSE;
2404 if(OPTIMIZING(sys))sys->mulstep.accurate = FALSE;
2405 }
2406
2407
2408 /*------------------------------------------------------------------------------
2409 BLOCK ROUTINES
2410 */
2411
2412 /**
2413 Returns TRUE if the current block is feasible, FALSE otherwise.
2414 It is assumed that the residuals have been computed.
2415 */
2416 static boolean block_feasible( qrslv_system_t sys){
2417 int32 row;
2418
2419 if(!sys->s.calc_ok )
2420 return(FALSE);
2421
2422 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2423 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2424 if(!rel_satisfied(rel) ) return FALSE;
2425 }
2426 return TRUE;
2427 }
2428
2429 /**
2430 Moves on to the next block, updating all of the solver information.
2431 To move to the first block, set sys->s.block.current_block to -1 before
2432 calling. If already at the last block, then sys->s.block.current_block
2433 will equal the number of blocks and the system will be declared
2434 converged. Otherwise, the residuals for the new block will be computed
2435 and sys->s.calc_ok set according.
2436 */
2437 static void move_to_next_block( qrslv_system_t sys){
2438 struct var_variable *var;
2439 struct rel_relation *rel;
2440 int32 row;
2441 int32 col;
2442 int32 ci;
2443 boolean ok;
2444
2445 if(sys->s.block.current_block >= 0 ) {
2446
2447
2448 /* Record cost accounting info here. */
2449 ci=sys->s.block.current_block;
2450 sys->s.cost[ci].size = sys->s.block.current_size;
2451 sys->s.cost[ci].iterations = sys->s.block.iteration;
2452 sys->s.cost[ci].funcs = sys->s.block.funcs;
2453 sys->s.cost[ci].jacs = sys->s.block.jacs;
2454 sys->s.cost[ci].functime = sys->s.block.functime;
2455 sys->s.cost[ci].jactime = sys->s.block.jactime;
2456 sys->s.cost[ci].time = sys->s.block.cpu_elapsed;
2457 sys->s.cost[ci].resid = sys->s.block.residual;
2458
2459 /* De-initialize previous block */
2460 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT) && (sys->s.block.current_size >1 ||
2461 SLV_PARAM_BOOL(&(sys->p),LIFDS))) {
2462 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"Block %d converged.\n",
2463 sys->s.block.current_block);
2464 }
2465 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
2466 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2467 var_set_in_block(var,FALSE);
2468 }
2469 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2470 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2471 rel_set_in_block(rel,FALSE);
2472 }
2473 sys->s.block.previous_total_size += sys->s.block.current_size;
2474 }
2475
2476 sys->s.block.current_block++;
2477 if(sys->s.block.current_block < sys->s.block.number_of ) {
2478
2479 /* Initialize next block */
2480 if(OPTIMIZING(sys)){
2481 mtx_region(&(sys->J.reg), 0, sys->rank-1, 0, sys->vused-1 );
2482 }else{
2483 sys->J.reg =
2484 (slv_get_solvers_blocks(SERVER))->block[sys->s.block.current_block];
2485 }
2486
2487 row = sys->J.reg.row.high - sys->J.reg.row.low + 1;
2488 col = sys->J.reg.col.high - sys->J.reg.col.low + 1;
2489 sys->s.block.current_size = MAX(row,col);
2490
2491 sys->s.block.iteration = 0;
2492 sys->s.block.cpu_elapsed = 0.0;
2493 sys->s.block.functime = 0.0;
2494 sys->s.block.jactime = 0.0;
2495 sys->s.block.funcs = 0;
2496 sys->s.block.jacs = 0;
2497
2498 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT) && (SLV_PARAM_BOOL(&(sys->p),LIFDS) ||
2499 sys->s.block.current_size > 1)) {
2500 debug_delimiter(LIF(sys));
2501 debug_delimiter(LIF(sys));
2502 }
2503 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT) && SLV_PARAM_BOOL(&(sys->p),LIFDS)) {
2504 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\n%-40s ---> %d in [%d..%d]\n",
2505 "Current block number", sys->s.block.current_block,
2506 0, sys->s.block.number_of-1);
2507 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d\n", "Current block size",
2508 sys->s.block.current_size);
2509 }
2510 sys->s.calc_ok = TRUE;
2511
2512 if(!(ok = calc_objective(sys)) ) {
2513 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Objective calculation errors detected");
2514 }
2515 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT) && sys->obj) {
2516 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", "Objective", sys->objective);
2517 }
2518 sys->s.calc_ok = sys->s.calc_ok && ok;
2519
2520 if(!(sys->p.ignore_bounds) ) {
2521 slv_ensure_bounds(SERVER, sys->J.reg.col.low,
2522 sys->J.reg.col.high,MIF(sys));
2523 }
2524
2525 sys->residuals.accurate = FALSE;
2526 if(!(ok = calc_residuals(sys)) ) {
2527 /* error_reporter will have been called somewhere else already */
2528 CONSOLE_DEBUG("Residual calculation errors detected in move_to_next_block.");
2529 }
2530 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT) &&
2531 (sys->s.block.current_size >1 ||
2532 SLV_PARAM_BOOL(&(sys->p),LIFDS)) ) {
2533 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", "Residual norm (unscaled)",
2534 sys->s.block.residual);
2535 }
2536 sys->s.calc_ok = sys->s.calc_ok && ok;
2537
2538 /* Must be updated as soon as required */
2539 sys->J.accurate = FALSE;
2540 sys->update.weights = 0;
2541 sys->update.nominals = 0;
2542 sys->update.relnoms = 0;
2543 sys->update.iterative = 0;
2544 sys->ZBZ.accurate = FALSE;
2545 sys->variables.accurate = FALSE;
2546 sys->gradient.accurate = FALSE;
2547 sys->multipliers.accurate = FALSE;
2548 sys->stationary.accurate = FALSE;
2549 sys->newton.accurate = FALSE;
2550 sys->Bnewton.accurate = FALSE;
2551 sys->nullspace.accurate = FALSE;
2552 sys->gamma.accurate = FALSE;
2553 sys->Jgamma.accurate = FALSE;
2554 sys->varstep1.accurate = FALSE;
2555 sys->Bvarstep1.accurate = FALSE;
2556 sys->varstep2.accurate = FALSE;
2557 sys->Bvarstep2.accurate = FALSE;
2558 sys->mulstep1.accurate = FALSE;
2559 sys->mulstep2.accurate = FALSE;
2560 sys->varstep.accurate = FALSE;
2561 sys->mulstep.accurate = FALSE;
2562
2563 if(!OPTIMIZING(sys)){
2564 sys->ZBZ.accurate = TRUE;
2565 sys->gradient.accurate = TRUE;
2566 sys->multipliers.accurate = TRUE;
2567 sys->stationary.accurate = TRUE;
2568 sys->Bnewton.accurate = TRUE;
2569 sys->nullspace.accurate = TRUE;
2570 sys->Bvarstep1.accurate = TRUE;
2571 sys->Bvarstep2.accurate = TRUE;
2572 }
2573
2574 }else{
2575 /*
2576 * Before we claim convergence, we must check if we left behind
2577 * some unassigned relations. If and only if they happen to be
2578 * satisfied at the current point, convergence has been obtained.
2579 *
2580 * Also insures that all included relations have valid residuals.
2581 * Included inequalities will have correct residuals.
2582 * Unsatisfied included inequalities cause inconsistency.
2583 *
2584 * This of course ignores that fact an objective function might
2585 * be present. Then, feasibility isn't enough, is it now.
2586 */
2587 if(sys->s.struct_singular ) {
2588 /* black box w/singletons provoking bug here, maybe */
2589 sys->s.block.current_size = sys->rused - sys->rank;
2590 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2591 debug_delimiter(LIF(sys));
2592 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d\n", "Unassigned Relations",
2593 sys->s.block.current_size);
2594 }
2595 sys->J.reg.row.low = sys->J.reg.col.low = sys->rank;
2596 sys->J.reg.row.high = sys->J.reg.col.high = sys->rused - 1;
2597 sys->residuals.accurate = FALSE;
2598 if(!(ok=calc_residuals(sys)) ) {
2599 FPRINTF(MIF(sys),
2600 "Residual calculation errors detected in leftover equations.\n");
2601 }
2602
2603 /** @TODO does this 'ok' needed to be ANDed with sys->s.calc_ok? */
2604
2605 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2606 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n", "Residual norm (unscaled)",
2607 sys->s.block.residual);
2608 }
2609 if(block_feasible(sys) ) {
2610 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2611 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nUnassigned relations ok. Lucky you.\n");
2612 }
2613 sys->s.converged = TRUE;
2614 }else{
2615 ERROR_REPORTER_HERE(ASC_PROG_WARNING,"Problem inconsistent: unassigned relations not satisfied");
2616 /* if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2617 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\nProblem inconsistent: %s.\n",
2618 "Unassigned relations not satisfied");
2619 }
2620 */
2621 sys->s.inconsistent = TRUE;
2622 }
2623 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)) {
2624 debug_delimiter(LIF(sys));
2625 }
2626 }else{
2627 sys->s.converged = TRUE;
2628 }
2629 /* nearly done checking. Must verify included inequalities if
2630 we think equalities are ok. */
2631 if(sys->s.converged) {
2632 ok = calc_inequalities(sys);
2633 if(!ok && sys->s.inconsistent){
2634 sys->s.inconsistent = TRUE;
2635 ERROR_REPORTER_HERE(ASC_PROG_ERR,"System marked inconsistent after inspecting inequalities");
2636 }
2637 }
2638 }
2639 }
2640
2641 /**
2642 Calls the appropriate reorder function on a block
2643 */
2644 static void reorder_new_block(qrslv_system_t sys){
2645 int32 method;
2646 if(sys->s.block.current_block < sys->s.block.number_of ) {
2647 if(strcmp(SLV_PARAM_CHAR(&(sys->p),REORDER_OPTION),"SPK1") == 0) {
2648 method = 2;
2649 }else{
2650 method = 1;
2651 }
2652
2653 if(sys->s.block.current_block <= sys->s.block.current_reordered_block &&
2654 sys->s.cost[sys->s.block.current_block].reorder_method == method &&
2655 sys->s.block.current_block >= 0 ) {
2656 #if DEBUG
2657 FPRINTF(stderr,"YOU JUST AVOIDED A REORDERING\n");
2658 #endif
2659 slv_set_up_block(SERVER,sys->s.block.current_block);
2660 /* tell linsol to bless it and get on with things */
2661 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
2662 return; /*must have been reordered since last system build*/
2663 }
2664
2665 /* Let the slv client function take care of reordering things
2666 * and setting in block flags.
2667 */
2668 if(strcmp(SLV_PARAM_CHAR(&(sys->p),REORDER_OPTION),"SPK1") == 0) {
2669 sys->s.cost[sys->s.block.current_block].reorder_method = 2;
2670 slv_spk1_reorder_block(SERVER,sys->s.block.current_block,1);
2671 }else if(strcmp(SLV_PARAM_CHAR(&(sys->p),REORDER_OPTION),"TEAR_DROP") == 0) {
2672 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2673 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block
2674 ,SLV_PARAM_INT(&(sys->p),CUTOFF), 0,mtx_SPK1
2675 );
2676 /* khack: try tspk1 for transpose case */
2677 }else if(strcmp(SLV_PARAM_CHAR(&(sys->p),REORDER_OPTION),"OVER_TEAR") == 0) {
2678 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2679 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block
2680 ,SLV_PARAM_INT(&(sys->p),CUTOFF), 1,mtx_SPK1
2681 );
2682 }else{
2683 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2684 ERROR_REPORTER_START_NOLINE(ASC_PROG_ERROR);
2685 FPRINTF(MIF(sys),"QRSlv called with unknown reorder option\n");
2686 FPRINTF(MIF(sys),"QRSlv using single edge tear drop (TEAR_DROP).\n");
2687 error_reporter_end_flush();
2688
2689 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2690 SLV_PARAM_INT(&(sys->p),CUTOFF),0,mtx_SPK1);
2691 }
2692 /* tell linsol to bless it and get on with things */
2693 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
2694 if(sys->s.block.current_block > sys->s.block.current_reordered_block) {
2695 sys->s.block.current_reordered_block = sys->s.block.current_block;
2696 }
2697 }
2698 }
2699
2700 /**
2701 Moves to next unconverged block, assuming that the current block has
2702 converged (or is -1, to start).
2703 */
2704 static void find_next_unconverged_block( qrslv_system_t sys){
2705
2706 do{
2707 move_to_next_block(sys);
2708 #if DEBUG
2709 debug_out_var_values(stderr,sys);
2710 debug_out_rel_residuals(stderr,sys);
2711 #endif
2712 }while(!sys->s.converged && block_feasible(sys) && !OPTIMIZING(sys));
2713
2714 reorder_new_block(sys);
2715 }
2716
2717 /*------------------------------------------------------------------------------
2718 ITERATION BEGIN/END ROUTINES
2719 */
2720
2721 /**
2722 Prepares sys for entering an iteration, increasing the iteration counts
2723 and starting the clock.
2724 */
2725 static void iteration_begins( qrslv_system_t sys){
2726 sys->clock = tm_cpu_time();
2727 ++(sys->s.block.iteration);
2728 ++(sys->s.iteration);
2729 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT)&& (sys->s.block.current_size >1 ||
2730 SLV_PARAM_BOOL(&(sys->p),LIFDS))) {
2731 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"\n%-40s ---> %d\n",
2732 "Iteration", sys->s.block.iteration);
2733 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %d\n",
2734 "Total iteration", sys->s.iteration);
2735 }
2736 }
2737
2738 /**
2739 Prepares sys for exiting an iteration, stopping the clock and recording
2740 the cpu time.
2741 */
2742 static void iteration_ends( qrslv_system_t sys){
2743 double cpu_elapsed; /* elapsed this iteration */
2744
2745 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
2746 sys->s.block.cpu_elapsed += cpu_elapsed;
2747 sys->s.cpu_elapsed += cpu_elapsed;
2748 if(SLV_PARAM_BOOL(&(sys->p),SHOW_LESS_IMPT) && (sys->s.block.current_size >1 ||
2749 SLV_PARAM_BOOL(&(sys->p),LIFDS))) {
2750 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2751 "Elapsed time", sys->s.block.cpu_elapsed);
2752 ERROR_REPORTER_HERE(ASC_PROG_NOTE,"%-40s ---> %g\n",
2753 "Total elapsed time", sys->s.cpu_elapsed);
2754 }
2755 }
2756
2757 /**
2758 Updates the solver status.
2759 */
2760 static void update_status( qrslv_system_t sys){
2761 boolean unsuccessful;
2762
2763 if(!sys->s.converged ) {
2764 sys->s.time_limit_exceeded =
2765 (sys->s.block.cpu_elapsed >= SLV_PARAM_INT(&(sys->p),TIME_LIMIT));
2766 sys->s.iteration_limit_exceeded =
2767 (sys->s.block.iteration >= SLV_PARAM_INT(&(sys->p),ITER_LIMIT));
2768 }
2769