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

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

Parent Directory Parent Directory | Revision Log Revision Log


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