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

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

Parent Directory Parent Directory | Revision Log Revision Log


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