/[ascend]/trunk/ipslv-temp/slv5.mehrotraquad.c
ViewVC logotype

Contents of /trunk/ipslv-temp/slv5.mehrotraquad.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 916 - (show annotations) (download) (as text)
Tue Oct 31 02:39:16 2006 UTC (13 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 93782 byte(s)
Adding Vicente's IPSlv code (needs revising and testing against the new code)\
1 /*
2 * IPSlv ASCEND Interior Point Method Solver
3 * by Vicente Rico-Ramirez based on QRSlv
4 * Created:
5 * Version: $ $
6 * Version control file: $ $
7 * Date last modified: $ $
8 * Last modified by: $Author: $
9 *
10 * This file is part of the SLV solver.
11 *
12 * The SLV solver is free software; you can redistribute
13 * it and/or modify it under the terms of the GNU General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
16 *
17 * The SLV solver is distributed in hope that it will be
18 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with the program; if not, write to the Free Software Foundation,
24 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
25 * COPYING. COPYING is found in ../compiler.
26 *
27 */
28
29 #include <math.h>
30 #include <stdarg.h>
31 #include "utilities/ascConfig.h"
32 #include "utilities/ascSignal.h"
33 #include "utilities/ascMalloc.h"
34 #include "utilities/set.h"
35 #include "general/time.h"
36 #include "utilities/mem.h"
37 #include "utilities/ascPanic.h"
38 #include "general/list.h"
39 #include "compiler/fractions.h"
40 #include "compiler/dimen.h"
41 #include "compiler/functype.h"
42 #include "compiler/func.h"
43 #include "solver/mtx.h"
44 #include "solver/linsol.h"
45 #include "solver/linsolqr.h"
46 #include "solver/slv_types.h"
47 #include "solver/var.h"
48 #include "solver/rel.h"
49 #include "solver/discrete.h"
50 #include "solver/conditional.h"
51 #include "solver/logrel.h"
52 #include "solver/bnd.h"
53 #include "solver/calc.h"
54 #include "solver/relman.h"
55 #include "solver/slv_common.h"
56 #include "solver/slv_client.h"
57 #include "solver/slv5.h"
58 #include "solver/slv_stdcalls.h"
59
60 #if !defined(STATIC_IPSLV) && !defined(DYNAMIC_IPSLV)
61 int slv5_register(SlvFunctionsT *f)
62 {
63 (void)f; /* stop gcc whine about unused parameter */
64
65 FPRINTF(ASCERR,"IPSlv not compiled in this ASCEND IV.\n");
66 return 1;
67 }
68 #else /* either STATIC_IPSLV or DYNAMIC_IPSLV is defined */
69 #ifdef DYNAMIC_IPSLV
70 /* do dynamic loading stuff. yeah, right */
71 #else /* following is used if STATIC_IPSLV is defined */
72
73 #define DEBUG TRUE
74 #define DEBUG_COMPLEMENTARY_VAR TRUE
75 #define DEBUG_OBJ_VALUES TRUE
76 #define DEBUG_ITERATION TRUE
77 #define DEBUG_CENTERING TRUE
78
79 #define SLV5(s) ((slv5_system_t)(s))
80 #define SERVER (sys->slv)
81 #define slv5_PA_SIZE 29 /* MUST INCREMENT WHEN ADDING PARAMETERS */
82 #define slv5_RA_SIZE 4
83
84 /* do not delete (or extend) this array definition.
85 */
86 #define IEX(n) slv5_iaexpln[(n)]
87 #define slv5_IA_SIZE 10
88 static char *slv5_iaexpln[slv5_IA_SIZE] = {
89 "If lifds != 0 and showlessimportant is TRUE, show direct solve details",
90 "If savlin != 0, write out matrix data file at each iteration to SlvLinsol.dat",
91 "Cutoff is the block size cutoff for MODEL-based reordering of partitions",
92 "Update jacobian every this many major iterations",
93 "Update row scalings every this many major iterations",
94 "Update column scalings every this many major iterations",
95 "Check jacobian for poorly scaled columns and whine if found",
96 "Reorder option. 0 = MODEL based, 1 = MODEL based2, 2 = simple spk1",
97 "Use safe calculation routines",
98 "scaleopt = 0: 2norm"
99 };
100
101 /* change slv5_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
102 #define UPDATE_JACOBIAN_PTR (sys->parm_array[0])
103 #define UPDATE_JACOBIAN ((*(int *)UPDATE_JACOBIAN_PTR))
104 #define UPDATE_WEIGHTS_PTR (sys->parm_array[1])
105 #define UPDATE_WEIGHTS ((*(int *)UPDATE_WEIGHTS_PTR))
106 #define UPDATE_NOMINALS_PTR (sys->parm_array[2])
107 #define UPDATE_NOMINALS ((*(int *)UPDATE_NOMINALS_PTR))
108 #define DUMPCNORM_PTR (sys->parm_array[3])
109 #define DUMPCNORM ((*(int *)DUMPCNORM_PTR))
110 #define SAFE_CALC_PTR (sys->parm_array[4])
111 #define SAFE_CALC ((*(int *)SAFE_CALC_PTR))
112 #define SCALEOPT_PTR (sys->parm_array[5])
113 #define SCALEOPT ((*(char **)SCALEOPT_PTR))
114 #define TOO_SMALL_PTR (sys->parm_array[6])
115 #define TOO_SMALL ((*(real64 *)TOO_SMALL_PTR))
116 #define CNLOW_PTR (sys->parm_array[7])
117 #define CNLOW ((*(real64 *)CNLOW_PTR))
118 #define CNHIGH_PTR (sys->parm_array[8])
119 #define CNHIGH ((*(real64 *)CNHIGH_PTR))
120 #define TOWARD_BOUNDS_PTR (sys->parm_array[9])
121 #define TOWARD_BOUNDS ((*(real64 *)TOWARD_BOUNDS_PTR))
122 #define IGNORE_BOUNDS_PTR (sys->parm_array[10])
123 #define IGNORE_BOUNDS ((*(int32 *)IGNORE_BOUNDS_PTR))
124 #define RHO_PTR (sys->parm_array[11])
125 #define RHO ((*(real64 *)RHO_PTR))
126 #define PARTITION_PTR (sys->parm_array[12])
127 #define PARTITION ((*(int32 *)PARTITION_PTR))
128 #define SHOW_LESS_IMPT_PTR (sys->parm_array[13])
129 #define SHOW_LESS_IMPT ((*(int32 *)SHOW_LESS_IMPT_PTR))
130 #define TIME_LIMIT_PTR (sys->parm_array[14])
131 #define TIME_LIMIT ((*(int32 *)TIME_LIMIT_PTR))
132 #define ITER_LIMIT_PTR (sys->parm_array[15])
133 #define ITER_LIMIT ((*(int32 *)ITER_LIMIT_PTR))
134 #define ALPHA_PTR (sys->parm_array[16])
135 #define ALPHA ((*(real64 *)ALPHA_PTR))
136 #define SING_TOL_PTR (sys->parm_array[17])
137 #define SING_TOL ((*(real64 *)SING_TOL_PTR))
138 #define PIVOT_TOL_PTR (sys->parm_array[18])
139 #define PIVOT_TOL ((*(real64 *)PIVOT_TOL_PTR))
140 #define FEAS_TOL_PTR (sys->parm_array[19])
141 #define FEAS_TOL ((*(real64 *)FEAS_TOL_PTR))
142 #define LIFDS_PTR (sys->parm_array[20])
143 #define LIFDS ((*(int32 *)LIFDS_PTR))
144 #define SAVLIN_PTR (sys->parm_array[21])
145 #define SAVLIN ((*(int32 *)SAVLIN_PTR))
146 #define REORDER_OPTION_PTR (sys->parm_array[22])
147 #define REORDER_OPTION ((*(char **)REORDER_OPTION_PTR))
148 #define CUTOFF_PTR (sys->parm_array[23])
149 #define CUTOFF ((*(int32 *)CUTOFF_PTR))
150 #define FACTOR_OPTION_PTR (sys->parm_array[24])
151 #define FACTOR_OPTION ((*(char **)FACTOR_OPTION_PTR))
152 #define CONVOPT_PTR (sys->parm_array[25])
153 #define CONVOPT ((*(char **)CONVOPT_PTR))
154 #define LINTIME_PTR (sys->parm_array[26])
155 #define LINTIME ((*(int32 *)LINTIME_PTR))
156 #define ITER_BIS_LIMIT_PTR (sys->parm_array[27])
157 #define ITER_BIS_LIMIT ((*(int32 *)ITER_BIS_LIMIT_PTR))
158 #define BIS_TOL_PTR (sys->parm_array[28])
159 #define BIS_TOL ((*(real64 *)BIS_TOL_PTR))
160 /* change slv5_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
161
162
163 #define REX(n) slv5_raexpln[(n)]
164 static
165 char *slv5_raexpln[slv5_RA_SIZE] = {
166 "Var nominal to use if user specifies 0.0",
167 "Smallest column norm we won't complain about if checking",
168 "Largest column norm we won't complain about if checking",
169 "If bound is in the way, we go this fraction toward it"};
170
171 struct update_data {
172 int jacobian; /* Countdown on jacobian updating */
173 int weights; /* Countdown on weights updating */
174 int nominals; /* Countdown on nominals updating */
175 };
176
177
178 struct jacobian_data {
179 linsolqr_system_t sys; /* Linear system */
180 mtx_matrix_t mtx; /* Transpose gradient of residuals */
181 real64 *rhs; /* RHS from linear system */
182 dof_t *dofdata; /* dof data pointer from server */
183 mtx_region_t reg; /* Current block region */
184 int32 rank; /* Numerical rank of the jacobian */
185 enum factor_method fm; /* Linear factorization method */
186 boolean accurate; /* ? Recalculate matrix */
187 boolean singular; /* ? Can matrix be inverted */
188 boolean old_partition;/* old value of partition flag */
189 };
190
191 struct slv5_system_structure {
192
193 /*
194 * Problem definition
195 */
196 slv_system_t slv; /* slv_system_t back-link */
197 struct rel_relation *obj; /* Objective function: NULL = none */
198 struct var_variable **vlist; /* Variable list (NULL terminated) */
199 struct rel_relation **rlist; /* Relation list (NULL terminated) */
200
201 /*
202 * Solver information
203 */
204 int integrity; /* ? Has the system been created */
205 int32 presolved; /* ? Has the system been presolved */
206 slv_parameters_t p; /* Parameters */
207 slv_status_t s; /* Status (as of iteration end) */
208 struct update_data update; /* Jacobian frequency counters */
209 int32 cap; /* Order of matrix/vectors */
210 int32 rank; /* Symbolic rank of problem */
211 int32 vused; /* Free and incident variables */
212 int32 vtot; /* length of varlist */
213 int32 rused; /* Included relations */
214 int32 rtot; /* length of rellist */
215 double clock; /* CPU time */
216 void *parm_array[slv5_PA_SIZE]; /* array of pointers to param values */
217 struct slv_parameter pa[slv5_PA_SIZE];/* &pa[0] => sys->p.parms */
218
219 /*
220 * Calculated data (scaled)
221 */
222 struct jacobian_data J; /* linearized system */
223 struct vector_data nominals; /* Variable nominals */
224 struct vector_data weights; /* Relation weights */
225 struct vector_data variables; /* Variable values */
226 struct vector_data residuals; /* Relation residuals */
227 struct vector_data newton_residuals; /* Newton Relation residuals */
228 struct vector_data perturbed_residuals; /* Perturbed residuals */
229 struct vector_data correction; /* 2nd order correction */
230 struct vector_data newton; /* Dependent variables */
231 struct vector_data perturbed_newton; /* Perturbed Newton direction */
232 struct vector_data varnewstep; /* newton step in variables */
233 struct vector_data varstep; /* Step in variables */
234
235 real64 progress; /* expected progress */
236 real64 sigma; /* penalty parameter */
237 real64 mu; /* Complementary gap */
238 real64 muaff; /* Complementary gap after newton */
239 real64 sigmamu; /* Complementary gap times penalty */
240 real64 eta; /* penalty of objective */
241 real64 nu; /* objective residuals */
242 real64 psi; /* Objective */
243 real64 comp; /* no. of complementary eqns. */
244 };
245
246
247 /*
248 * Integrity checks
249 * ----------------
250 * check_system(sys)
251 */
252
253 #define OK ((int)813029392)
254 #define DESTROYED ((int)103289182)
255 static int check_system(slv5_system_t sys)
256 /*
257 * Checks sys for NULL and for integrity.
258 */
259 {
260 if( sys == NULL ) {
261 FPRINTF(ASCERR,"ERROR: (slv5) check_system\n");
262 FPRINTF(ASCERR," NULL system handle.\n");
263 return 1;
264 }
265
266 switch( sys->integrity ) {
267 case OK:
268 return 0;
269 case DESTROYED:
270 FPRINTF(ASCERR,"ERROR: (slv5) check_system\n");
271 FPRINTF(ASCERR," System was recently destroyed.\n");
272 return 1;
273 default:
274 FPRINTF(ASCERR,"ERROR: (slv5) check_system\n");
275 FPRINTF(ASCERR," System reused or never allocated.\n");
276 return 1;
277 }
278 }
279
280 /*
281 * General input/output routines
282 * -----------------------------
283 * print_var_name(out,sys,var)
284 * print_rel_name(out,sys,rel)
285 */
286
287 #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
288 #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
289
290 /*
291 * Debug output routines
292 * ---------------------
293 * debug_delimiter(fp)
294 * debug_out_vector(fp,vec)
295 * debug_out_var_values(fp,sys)
296 * debug_out_rel_residuals(fp,sys)
297 * debug_out_jacobian(fp,sys)
298 * debug_write_array(fp,real64 *,length)
299 * debug_out_parameters(fp)
300 */
301
302 /*
303 * Outputs a hyphenated line.
304 */
305 static void debug_delimiter( FILE *fp)
306 {
307 int i;
308 for( i=0; i<60; i++ ) PUTC('-',fp);
309 PUTC('\n',fp);
310 }
311
312 #if DEBUG
313 /*
314 * Outputs a vector.
315 */
316 static void debug_out_vector( FILE *fp, struct vector_data *vec)
317 {
318 int32 ndx;
319 FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n",
320 calc_sqrt_D0(vec->norm2), vec->accurate?"TRUE":"FALSE",
321 vec->rng->low,vec->rng->high);
322 FPRINTF(fp,"Vector --> ");
323 for( ndx=vec->rng->low ; ndx<=vec->rng->high ; ++ndx )
324 FPRINTF(fp, "%g ", vec->vec[ndx]);
325 PUTC('\n',fp);
326 }
327
328 /*
329 * Outputs all variable values in current block.
330 */
331 static void debug_out_var_values( FILE *fp, slv5_system_t sys)
332 {
333 int32 col;
334 struct var_variable *var;
335
336 FPRINTF(fp,"Var values --> \n");
337 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high ; col++ ) {
338 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
339 print_var_name(fp,sys,var);
340 FPRINTF(fp, "\nI Lb Value Ub Scale Col INom\n");
341 FPRINTF(fp,"%d\t%.4g\t%.4g\t%.4g\t%.4g\t%d\t%.4g\n",
342 var_sindex(var),var_lower_bound(var),var_value(var),
343 var_upper_bound(var),var_nominal(var),
344 col,sys->nominals.vec[col]);
345 }
346 }
347
348 /*
349 * Outputs all relation residuals in current block.
350 */
351 static void debug_out_rel_residuals( FILE *fp, slv5_system_t sys)
352 {
353 int32 row;
354
355 FPRINTF(fp,"Rel residuals --> \n");
356 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high ; row++ ) {
357 struct rel_relation *rel;
358 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
359 FPRINTF(fp," %g : ",rel_residual(rel));
360 print_rel_name(fp,sys,rel);
361 PUTC('\n',fp);
362 }
363 PUTC('\n',fp);
364 }
365
366 /*
367 * Outputs permutation and values of the nonzero elements in the
368 * the jacobian matrix.
369 */
370 static void debug_out_jacobian( FILE *fp, slv5_system_t sys)
371 {
372 mtx_coord_t nz;
373 real64 value;
374
375 nz.row = sys->J.reg.row.low;
376 for( ; nz.row <= sys->J.reg.row.high; ++(nz.row) ) {
377 FPRINTF(fp," Row %d (rel %d)\n", nz.row,
378 mtx_row_to_org(sys->J.mtx,nz.row));
379 nz.col = mtx_FIRST;
380 while( value = mtx_next_in_row(sys->J.mtx,&nz,&(sys->J.reg.col)),
381 nz.col != mtx_LAST ) {
382 FPRINTF(fp," Col %d (var %d) has value %g\n", nz.col,
383 mtx_col_to_org(sys->J.mtx,nz.col), value);
384 }
385 }
386 }
387
388 #endif
389
390 static void debug_write_array(FILE *fp,real64 *vec, int32 length)
391 {
392 int32 i;
393 for (i=0; i< length;i++)
394 FPRINTF(fp,"%.20g\n",vec[i]);
395 }
396
397 static char savlinfilename[]="SlvLinsol.dat.\0";
398 static char savlinfilebase[]="SlvLinsol.dat.\0";
399 static int savlinnum=0;
400
401 /* The number to postfix to savlinfilebase. increases with file accesses. */
402
403 /*
404 * Array/vector operations
405 * ----------------------------
406 * destroy_array(p)
407 * create_array(len,type)
408 * zero_vector(vec)
409 * copy_vector(vec1,vec2)
410 * prod = inner_product(vec1,vec2)
411 * norm2 = square_norm(vec)
412 * matrix_product(mtx,vec,prod,scale,transpose)
413 */
414
415 #define destroy_array(p) \
416 if( (p) != NULL ) ascfree((p))
417 #define create_array(len,type) \
418 ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
419 #define create_zero_array(len,type) \
420 ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
421
422 #define zero_vector(v) slv_zero_vector(v)
423 #define copy_vector(v,t) slv_copy_vector((v),(t))
424 #define inner_product(v,u) slv_inner_product((v),(u))
425 #define square_norm(v) slv_square_norm(v)
426 #define matrix_product(m,v,p,s,t) slv_matrix_product((m),(v),(p),(s),(t))
427
428 /*
429 * Calculation routines
430 * --------------------
431 * calc_eta(server)
432 * ok = calc_residuals(sys)
433 * ok = calc_mu(sys)
434 * ok = calc_newton_residuals(sys)
435 * ok = calc_muaff(sys)
436 * ok = calc_sigma(sys)
437 * ok = calc_sigmamu(sys)
438 * ok = calc_perturbed_residuals(sys)
439 * ok = calc_J(sys)
440 * calc_nominals(sys)
441 * calc_weights(sys)
442 * scale_J(sys)
443 * scale_variables(sys)
444 * scale_residuals(sys)
445 * scale_perturbed_residuals(sys)
446 * calc_pivots(sys)
447 * calc_rhs(sys)
448 * calc_newton(sys)
449 * calc_perturbed_newton(sys)
450 * calc_varnewstep(sys)
451 * calc_varstep(sys)
452 * calc_nu(sys)
453 * calc_psi(sys)
454 */
455
456
457 /*
458 * Calculates the penalty of objective function
459 */
460 static void calc_comp( slv5_system_t sys )
461 {
462 int32 row;
463 struct rel_relation *rel;
464 real64 comp;
465
466 comp = 0.0;
467 row = sys->residuals.rng->low;
468 for( ; row <= sys->residuals.rng->high; row++ ) {
469 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
470 #if DEBUG
471 if (!rel) {
472 int r;
473 r=mtx_row_to_org(sys->J.mtx,row);
474 FPRINTF(ASCERR,"NULL relation found !!\n");
475 FPRINTF(ASCERR,"at row %d rel %d in calc_comp\n",(int)row,r);
476 FFLUSH(ASCERR);
477 }
478 #endif
479 if (rel_active(rel) && rel_included (rel) && rel_complementary(rel)) {
480 #if DEBUG
481 FPRINTF(ASCERR,"Complementary equation in slv5 \n");
482 #endif /* DEBUG */
483 comp = comp + 1.0;
484 }
485 }
486
487 #if DEBUG_ITERATION
488 FPRINTF(ASCERR," No. of complementary eqns = %g \n",comp);
489 #endif /* DEBUG_ITERATION */
490 sys->comp = comp;
491 }
492
493 /*
494 * Calculates the penalty of objective function
495 */
496 static void calc_eta( slv5_system_t sys)
497 {
498
499 real64 eta;
500
501 if (sys->comp == 0.0) {
502 eta = 0.0;
503 } else {
504 eta = sys->comp * sqrt(sys->comp);
505 }
506 #if DEBUG_ITERATION
507 FPRINTF(ASCERR," eta = %g \n",eta);
508 #endif /* DEBUG_ITERATION */
509 sys->eta = eta;
510 }
511
512
513 /*
514 * Calculates all of the residuals in the current block and computes
515 * the residual norm for block status. Returns true iff calculations
516 * preceded without error.
517 */
518 static boolean calc_residuals( slv5_system_t sys)
519 {
520 int32 row;
521 struct rel_relation *rel;
522 double time0;
523
524 if( sys->residuals.accurate ) return TRUE;
525
526 calc_ok = TRUE;
527 row = sys->residuals.rng->low;
528 time0=tm_cpu_time();
529 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
530 for( ; row <= sys->residuals.rng->high; row++ ) {
531 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
532 #if DEBUG
533 if (!rel) {
534 int r;
535 r=mtx_row_to_org(sys->J.mtx,row);
536 FPRINTF(ASCERR,"NULL relation found !!\n");
537 FPRINTF(ASCERR,"at row %d rel %d in calc_residuals\n",(int)row,r);
538 FFLUSH(ASCERR);
539 }
540 #endif
541 sys->residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
542
543 if (strcmp(CONVOPT,"ABSOLUTE") == 0) {
544 relman_calc_satisfied(rel,FEAS_TOL);
545 } else if (strcmp(CONVOPT,"RELNOM_SCALE") == 0) {
546 relman_calc_satisfied_scaled(rel,FEAS_TOL);
547 }
548 }
549 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
550 sys->s.block.functime += (tm_cpu_time() -time0);
551 sys->s.block.funcs++;
552 square_norm( &(sys->residuals) );
553 sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
554 return(calc_ok);
555 }
556
557 /*
558 * Calculates the complementary gap in the current point
559 */
560 static boolean calc_mu( slv5_system_t sys)
561 {
562
563 int32 row;
564 struct rel_relation *rel;
565 real64 muk;
566
567 muk = 0.0;
568 row = sys->residuals.rng->low;
569
570 #if DEBUG_CENTERING
571 FPRINTF(ASCERR,"row low is = %d \n",row);
572 FPRINTF(ASCERR,"row high is = %d \n",sys->residuals.rng->high);
573 #endif /* DEBUG_CENTERING */
574
575 for(row ; row <= sys->residuals.rng->high; row++ ) {
576 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
577 #if DEBUG
578 if (!rel) {
579 int r;
580 r = mtx_row_to_org(sys->J.mtx,row);
581 FPRINTF(ASCERR,"NULL relation found !!\n");
582 FPRINTF(ASCERR,"at row %d rel %d in calc_mu \n",(int)row,r);
583 FFLUSH(ASCERR);
584 }
585 #endif
586 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
587 muk = muk + rel_residual(rel);
588 #if DEBUG
589 FPRINTF(ASCERR,"Complementary equation in calc_mu \n");
590 FPRINTF(ASCERR,"row = %d\n",row);
591 FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
592 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
593 FPRINTF(ASCERR,"Partial muk is = %g \n",muk);
594 #endif /* DEBUG */
595 }
596 }
597
598 if (sys->comp > 0.0) {
599 muk = muk / sys->comp;
600 } else {
601 muk = 0.0;
602 }
603
604 sys->mu = muk;
605
606 #if DEBUG_CENTERING
607 FPRINTF(ASCERR,"muk is = %g \n",sys->mu);
608 #endif /* DEBUG_CENTERING */
609
610 return(calc_ok);
611 }
612
613
614 /*
615 * Calculates all of the residuals in the current block and computes
616 * the residual norm for block status after a hypothetical Newton
617 * step. Returns true iff calculations
618 * preceded without error.
619 */
620 static boolean calc_newton_residuals( slv5_system_t sys)
621 {
622 int32 row;
623 struct rel_relation *rel;
624
625 if( sys->newton_residuals.accurate ) return TRUE;
626
627 calc_ok = TRUE;
628 row = sys->newton_residuals.rng->low;
629 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
630 for( ; row <= sys->newton_residuals.rng->high; row++ ) {
631 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
632 #if DEBUG
633 if (!rel) {
634 int r;
635 r=mtx_row_to_org(sys->J.mtx,row);
636 FPRINTF(ASCERR,"NULL relation found !!\n");
637 FPRINTF(ASCERR,"at row %d rel %d in calc_newton_residuals\n",(int)row,r);
638 FFLUSH(ASCERR);
639 }
640 #endif
641 sys->newton_residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
642 }
643 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
644 square_norm( &(sys->newton_residuals) );
645 return(calc_ok);
646 }
647
648 /*
649 * Calculates the complementary gap after a hypothetical Newton Step
650 */
651 static boolean calc_muaff( slv5_system_t sys)
652 {
653 int32 row;
654 struct rel_relation *rel;
655 real64 muaff;
656
657 muaff = 0.0;
658 row = sys->newton_residuals.rng->low;
659 for(row ; row <= sys->newton_residuals.rng->high; row++ ) {
660 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
661 #if DEBUG
662 if (!rel) {
663 int r;
664 r=mtx_row_to_org(sys->J.mtx,row);
665 FPRINTF(ASCERR,"NULL relation found !!\n");
666 FPRINTF(ASCERR,"at row %d rel %d in calc_muaff\n",(int)row,r);
667 FFLUSH(ASCERR);
668 }
669 #endif
670 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
671 muaff = muaff + sys->newton_residuals.vec[row];
672 #if DEBUG
673 FPRINTF(ASCERR,"Complementary equation in calc_muaff \n");
674 FPRINTF(ASCERR,"row = %d\n",row);
675 FPRINTF(ASCERR,"residual vector = %g\n",sys->newton_residuals.vec[row]);
676 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
677 FPRINTF(ASCERR,"Partial muaff is = %g \n",muaff);
678 #endif /* DEBUG */
679 }
680 }
681
682 if (sys->comp > 0.0) {
683 muaff = muaff / sys->comp;
684 } else {
685 muaff = 0.0;
686 }
687
688 sys->muaff = muaff;
689
690 #if DEBUG_CENTERING
691 FPRINTF(ASCERR,"muaff is = %g \n",sys->muaff);
692 #endif /* DEBUG_CENTERING */
693
694 return(calc_ok);
695 }
696
697
698 /*
699 * Calculates the penalty parameter sigma
700 */
701 static boolean calc_sigma( slv5_system_t sys)
702 {
703 real64 sigma, frac;
704
705 if ((sys->mu) > 0.0) {
706 frac = (sys->muaff) / (sys->mu);
707 sigma = (frac) * (frac) * (frac);
708 } else {
709 frac = 0.0;
710 sigma = 0.0;
711 }
712
713 sys->sigma = sigma;
714
715 #if DEBUG_CENTERING
716 FPRINTF(ASCERR,"sigma is = %g \n",sys->sigma);
717 #endif /* DEBUG_CENTERING */
718
719 return(calc_ok);
720 }
721
722 /*
723 * Calculates the penalty parameter time the complementary gap
724 */
725 static boolean calc_sigmamu( slv5_system_t sys)
726 {
727 real64 sigmamu;
728
729 if ((sys->mu) > 0.0 && (sys->sigma) > 0.0) {
730 sigmamu = (sys->mu) * (sys->sigma);
731 } else {
732 sigmamu = 0.0;
733 }
734
735 sys->sigmamu = sigmamu;
736
737 #if DEBUG_CENTERING
738 FPRINTF(ASCERR,"sigma mu is = %g \n",sys->sigmamu);
739 #endif /* DEBUG_CENTERING */
740
741 return(calc_ok);
742 }
743
744
745 /*
746 * Calculates the 2nd order correction for Mehrotra's corrector
747 * step. Returns true iff calculations preceded without error.
748 */
749 static boolean calc_correction( slv5_system_t sys)
750 {
751 int32 row;
752 struct rel_relation *rel;
753
754 if( sys->correction.accurate ) return TRUE;
755
756 calc_ok = TRUE;
757 row = sys->correction.rng->low;
758 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
759 for( ; row <= sys->correction.rng->high; row++ ) {
760 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
761 #if DEBUG
762 if (!rel) {
763 int r;
764 r = mtx_row_to_org(sys->J.mtx,row);
765 FPRINTF(ASCERR,"NULL relation found !!\n");
766 FPRINTF(ASCERR,"at row %d rel %d in calc_correction\n",
767 (int)row,r);
768 FFLUSH(ASCERR);
769 }
770 #endif
771 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
772 sys->correction.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
773 #if DEBUG
774 FPRINTF(ASCERR,"calc_correction \n");
775 FPRINTF(ASCERR,"row = %d\n",row);
776 FPRINTF(ASCERR,"correction = %g\n",sys->correction.vec[row]);
777 #endif /* DEBUG */
778 } else {
779 sys->correction.vec[row] = 0.0;
780 }
781 }
782 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
783 square_norm( &(sys->correction) );
784 sys->correction.accurate = TRUE;
785 return(calc_ok);
786 }
787
788
789 /*
790 * Calculates the perturbed residuals in the current block and computes
791 * the perturbed residual norm. Returns true iff calculations
792 * preceded without error.
793 */
794 static boolean calc_perturbed_residuals( slv5_system_t sys)
795 {
796 int32 row;
797 struct rel_relation *rel;
798
799 if( sys->perturbed_residuals.accurate ) return TRUE;
800
801 calc_ok = TRUE;
802 row = sys->perturbed_residuals.rng->low;
803 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
804 for( ; row <= sys->perturbed_residuals.rng->high; row++ ) {
805 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
806 #if DEBUG
807 if (!rel) {
808 int r;
809 r=mtx_row_to_org(sys->J.mtx,row);
810 FPRINTF(ASCERR,"NULL relation found !!\n");
811 FPRINTF(ASCERR,"at row %d rel %d in calc_perturbed _residuals\n",
812 (int)row,r);
813 FFLUSH(ASCERR);
814 }
815 #endif
816 #if DEBUG
817 FPRINTF(ASCERR,"calc_perturbed_residuals \n");
818 FPRINTF(ASCERR,"row = %d\n",row);
819 #endif /* DEBUG */
820
821 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
822 sys->perturbed_residuals.vec[row] = sys->correction.vec[row]
823 - sys->sigmamu ;
824 #if DEBUG
825 FPRINTF(ASCERR,"Complementary equation \n");
826 FPRINTF(ASCERR,"sigmamu = %g \n",sys->sigmamu);
827 FPRINTF(ASCERR,"correction = %g \n",sys->correction.vec[row]);
828 FPRINTF(ASCERR,"- sigmamu + correction = %g\n",
829 sys->perturbed_residuals.vec[row]);
830 #endif /* DEBUG */
831 } else {
832 sys->perturbed_residuals.vec[row] = 0.0;
833 }
834 }
835 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
836 square_norm( &(sys->perturbed_residuals) );
837 return(calc_ok);
838 }
839
840 /*
841 * Calculates the current block of the jacobian.
842 * It is initially unscaled.
843 */
844 static boolean calc_J( slv5_system_t sys)
845 {
846 int32 row;
847 var_filter_t vfilter;
848 double time0;
849 real64 resid;
850
851 if( sys->J.accurate )
852 return TRUE;
853
854 calc_ok = TRUE;
855 vfilter.matchbits = (VAR_INBLOCK | VAR_ACTIVE);
856 vfilter.matchvalue = (VAR_INBLOCK | VAR_ACTIVE);
857 time0=tm_cpu_time();
858 mtx_clear_region(sys->J.mtx,&(sys->J.reg));
859 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
860 struct rel_relation *rel;
861 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
862 relman_diffs(rel,&vfilter,sys->J.mtx,&resid,SAFE_CALC);
863 }
864 sys->s.block.jactime += (tm_cpu_time() - time0);
865 sys->s.block.jacs++;
866
867 if( --(sys->update.nominals) <= 0 ) sys->nominals.accurate = FALSE;
868 if( --(sys->update.weights) <= 0 ) sys->weights.accurate = FALSE;
869
870 linsolqr_matrix_was_changed(sys->J.sys);
871 return(calc_ok);
872 }
873
874 /*
875 * Retrieves the nominal values of all of the block variables,
876 * insuring that they are all strictly positive.
877 */
878 static void calc_nominals( slv5_system_t sys)
879 {
880 int32 col;
881 FILE *fp = MIF(sys);
882 if( sys->nominals.accurate ) return;
883 fp = MIF(sys);
884 col = sys->nominals.rng->low;
885 if(strcmp(SCALEOPT,"NONE") == 0){
886 for( ; col <= sys->nominals.rng->high; col++ ) {
887 sys->nominals.vec[col] = 1;
888 }
889 } else {
890 for( ; col <= sys->nominals.rng->high; col++ ) {
891 struct var_variable *var;
892 real64 n;
893 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
894 n = var_nominal(var);
895 if( n <= 0.0 ) {
896 if( n == 0.0 ) {
897 n = TOO_SMALL;
898 FPRINTF(fp,"ERROR: (slv5) calc_nominals\n");
899 FPRINTF(fp," Variable ");
900 print_var_name(fp,sys,var);
901 FPRINTF(fp," \nhas nominal value of zero.\n");
902 FPRINTF(fp," Resetting to %g.\n",n);
903 var_set_nominal(var,n);
904 } else {
905 n = -n;
906 FPRINTF(fp,"ERROR: (slv5) calc_nominals\n");
907 FPRINTF(fp," Variable ");
908 print_var_name(fp,sys,var);
909 FPRINTF(fp," \nhas negative nominal value.\n");
910 FPRINTF(fp," Resetting to %g.\n",n);
911 var_set_nominal(var,n);
912 }
913 }
914 #if DEBUG
915 FPRINTF(fp,"Column %d is ",col);
916 print_var_name(fp,sys,var);
917 FPRINTF(fp,"\nScaling of column %d is %g\n",col,n);
918 #endif
919 sys->nominals.vec[col] = n;
920 }
921 }
922 square_norm( &(sys->nominals) );
923 sys->update.nominals = UPDATE_NOMINALS;
924 sys->nominals.accurate = TRUE;
925 }
926
927 /*
928 * Calculates the weights of all of the block relations
929 * to scale the rows of the Jacobian.
930 */
931 static void calc_weights( slv5_system_t sys)
932 {
933 mtx_coord_t nz;
934 real64 sum;
935
936 if( sys->weights.accurate )
937 return;
938
939 nz.row = sys->weights.rng->low;
940 if(strcmp(SCALEOPT,"NONE") == 0) {
941 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
942 sys->weights.vec[nz.row] = 1;
943 }
944 } else if (strcmp(SCALEOPT,"ROW_2NORM") == 0 ) {
945 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
946 sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
947 sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
948 }
949 }
950 square_norm( &(sys->weights) );
951 sys->update.weights = UPDATE_WEIGHTS;
952 sys->residuals.accurate = FALSE;
953 sys->weights.accurate = TRUE;
954 }
955
956 /*
957 * Scales the jacobian.
958 */
959 static void scale_J( slv5_system_t sys)
960 {
961 int32 row;
962 int32 col;
963
964 if( sys->J.accurate ) return;
965
966 calc_nominals(sys);
967 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ )
968 mtx_mult_col(sys->J.mtx,col,sys->nominals.vec[col],&(sys->J.reg.row));
969
970 calc_weights(sys);
971 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ )
972 mtx_mult_row(sys->J.mtx,row,sys->weights.vec[row],&(sys->J.reg.col));
973 }
974
975 static void jacobian_scaled(slv5_system_t sys)
976 {
977 int32 col;
978 if (DUMPCNORM) {
979 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
980 real64 cnorm;
981 cnorm =
982 calc_sqrt_D0(mtx_sum_sqrs_in_col(sys->J.mtx,col,&(sys->J.reg.row)));
983 if (cnorm >CNHIGH || cnorm <CNLOW) {
984 FPRINTF(ASCERR,"[col %d org %d] %g\n", col,
985 mtx_col_to_org(sys->J.mtx,col), cnorm);
986 }
987 }
988 }
989
990 sys->update.jacobian = UPDATE_JACOBIAN;
991 sys->J.accurate = TRUE;
992 sys->J.singular = FALSE; /* yet to be determined */
993 #if DEBUG
994 FPRINTF(LIF(sys),"\nJacobian: \n");
995 debug_out_jacobian(LIF(sys),sys);
996 #endif
997 }
998
999 static void scale_variables( slv5_system_t sys)
1000 {
1001 int32 col;
1002
1003 if( sys->variables.accurate ) return;
1004
1005 col = sys->variables.rng->low;
1006 for( ; col <= sys->variables.rng->high; col++ ) {
1007 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1008 sys->variables.vec[col] = var_value(var)/sys->nominals.vec[col];
1009 }
1010 square_norm( &(sys->variables) );
1011 sys->variables.accurate = TRUE;
1012 #if DEBUG
1013 FPRINTF(LIF(sys),"Variables: ");
1014 debug_out_vector(LIF(sys),&(sys->variables));
1015 #endif
1016 }
1017
1018 /*
1019 * Scales the previously calculated residuals.
1020 */
1021 static void scale_residuals( slv5_system_t sys)
1022 {
1023 int32 row;
1024
1025 if( sys->residuals.accurate ) return;
1026
1027 row = sys->residuals.rng->low;
1028 for( ; row <= sys->residuals.rng->high; row++ ) {
1029 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1030 sys->residuals.vec[row] = rel_residual(rel)*sys->weights.vec[row];
1031 }
1032 square_norm( &(sys->residuals) );
1033 sys->residuals.accurate = TRUE;
1034 #if DEBUG
1035 FPRINTF(LIF(sys),"Residuals: ");
1036 debug_out_vector(LIF(sys),&(sys->residuals));
1037 #endif
1038 }
1039
1040 /*
1041 * Scales the previously calculated residuals.
1042 */
1043 static void scale_perturbed_residuals( slv5_system_t sys)
1044 {
1045 int32 row;
1046
1047 if( sys->perturbed_residuals.accurate ) return;
1048
1049 row = sys->perturbed_residuals.rng->low;
1050 for( ; row <= sys->perturbed_residuals.rng->high; row++ ) {
1051 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1052 #if DEBUG
1053 FPRINTF(ASCERR,"scale_perturbed_residuals \n");
1054 FPRINTF(ASCERR,"row = %d\n",row);
1055 FPRINTF(ASCERR,"residual vector = %g\n",
1056 sys->perturbed_residuals.vec[row]);
1057 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1058 #endif /* DEBUG */
1059 sys->perturbed_residuals.vec[row] = sys->perturbed_residuals.vec[row]
1060 * sys->weights.vec[row];
1061 }
1062 square_norm( &(sys->perturbed_residuals) );
1063 sys->perturbed_residuals.accurate = TRUE;
1064 #if DEBUG
1065 FPRINTF(LIF(sys),"Perturbed Residuals: ");
1066 debug_out_vector(LIF(sys),&(sys->perturbed_residuals));
1067 #endif
1068 }
1069
1070
1071 /*
1072 * Scale system dependent on interface parameters
1073 */
1074 static void scale_perturbed_system( slv5_system_t sys )
1075 {
1076 if(strcmp(SCALEOPT,"NONE") == 0){
1077 scale_perturbed_residuals(sys);
1078 return;
1079 }
1080 if(strcmp(SCALEOPT,"ROW_2NORM") == 0){
1081 scale_perturbed_residuals(sys);
1082 return;
1083 }
1084 }
1085
1086 /*
1087 * Scale system dependent on interface parameters
1088 */
1089 static void scale_system( slv5_system_t sys )
1090 {
1091 if(strcmp(SCALEOPT,"NONE") == 0){
1092 if(sys->J.accurate == FALSE){
1093 calc_nominals(sys);
1094 calc_weights(sys);
1095 jacobian_scaled(sys);
1096 }
1097 scale_variables(sys);
1098 scale_residuals(sys);
1099 return;
1100 }
1101 if(strcmp(SCALEOPT,"ROW_2NORM") == 0 ){
1102 if(sys->J.accurate == FALSE){
1103 scale_J(sys);
1104 jacobian_scaled(sys);
1105 }
1106 scale_variables(sys);
1107 scale_residuals(sys);
1108 return;
1109 }
1110 }
1111
1112
1113 /*
1114 * Obtain the equations and variables which
1115 * are able to be pivoted.
1116 * return value is the row rank deficiency, which we hope is 0.
1117 */
1118 static int calc_pivots(slv5_system_t sys)
1119 {
1120 int row_rank_defect=0, oldtiming;
1121 linsolqr_system_t lsys = sys->J.sys;
1122 FILE *fp = LIF(sys);
1123
1124 oldtiming = g_linsolqr_timing;
1125 g_linsolqr_timing =LINTIME;
1126 linsolqr_factor(lsys,sys->J.fm); /* factor */
1127 g_linsolqr_timing = oldtiming;
1128
1129 sys->J.rank = linsolqr_rank(lsys);
1130 sys->J.singular = FALSE;
1131 row_rank_defect = sys->J.reg.row.high -
1132 sys->J.reg.row.low+1 - sys->J.rank;
1133 if( row_rank_defect > 0 ) {
1134 int32 row,krow;
1135 mtx_sparse_t *uprows=NULL;
1136 sys->J.singular = TRUE;
1137 uprows = linsolqr_unpivoted_rows(lsys);
1138 if (uprows !=NULL) {
1139 for( krow=0; krow < uprows->len ; krow++ ) {
1140 int32 org_row;
1141 struct rel_relation *rel;
1142
1143 org_row = uprows->idata[krow];
1144 row = mtx_org_to_row(sys->J.mtx,org_row);
1145 rel = sys->rlist[org_row];
1146 FPRINTF(fp,"%-40s ---> ","Relation not pivoted");
1147 print_rel_name(fp,sys,rel);
1148 PUTC('\n',fp);
1149
1150 /*
1151 * assign zeros to the corresponding weights
1152 * so that subsequent calls to "scale_residuals"
1153 * will only measure the pivoted equations.
1154 */
1155 sys->weights.vec[row] = 0.0;
1156 sys->residuals.vec[row] = 0.0;
1157 sys->residuals.accurate = FALSE;
1158 sys->correction.vec[row] = 0.0;
1159 sys->correction.accurate = FALSE;
1160 sys->perturbed_residuals.vec[row] = 0.0;
1161 sys->perturbed_residuals.accurate = FALSE;
1162 sys->newton_residuals.vec[row] = 0.0;
1163 sys->newton_residuals.accurate = FALSE;
1164 mtx_mult_row(sys->J.mtx,row,0.0,&(sys->J.reg.col));
1165 }
1166 mtx_destroy_sparse(uprows);
1167 }
1168 if( !sys->residuals.accurate ) {
1169 square_norm( &(sys->residuals) );
1170 sys->residuals.accurate = TRUE;
1171 sys->update.weights = 0; /* re-compute weights next iteration. */
1172 }
1173 }
1174 if( sys->J.rank < sys->J.reg.col.high-sys->J.reg.col.low+1 ) {
1175 int32 col,kcol;
1176 mtx_sparse_t *upcols=NULL;
1177 if (NOTNULL(upcols)) {
1178 for( kcol=0; upcols != NULL && kcol < upcols->len ; kcol++ ) {
1179 int32 org_col;
1180 struct var_variable *var;
1181
1182 org_col = upcols->idata[kcol];
1183 col = mtx_org_to_col(sys->J.mtx,org_col);
1184 var = sys->vlist[org_col];
1185 FPRINTF(fp,"%-40s ---> ","Variable not pivoted");
1186 print_var_name(fp,sys,var);
1187 PUTC('\n',fp);
1188 /*
1189 * If we're not optimizing (everything should be
1190 * pivotable) or this was one of the dependent variables,
1191 * consider this variable as if it were fixed.
1192 */
1193 if( col <= sys->J.reg.col.high ) {
1194 mtx_mult_col(sys->J.mtx,col,0.0,&(sys->J.reg.row));
1195 }
1196 }
1197 mtx_destroy_sparse(upcols);
1198 }
1199 }
1200 if (SHOW_LESS_IMPT) {
1201 FPRINTF(LIF(sys),"%-40s ---> %d (%s)\n","Jacobian rank", sys->J.rank,
1202 sys->J.singular ? "deficient":"full");
1203 FPRINTF(LIF(sys),"%-40s ---> %g\n","Smallest pivot",
1204 linsolqr_smallest_pivot(sys->J.sys));
1205 }
1206 return row_rank_defect;
1207 }
1208
1209
1210 /*
1211 * Calculates just the jacobian RHS. This function should be used to
1212 * supplement calculation of the jacobian. The vector vec must
1213 * already be calculated and scaled so as to simply be added to the
1214 * rhs. Caller is responsible for initially zeroing the rhs vector.
1215 */
1216 static void calc_rhs(slv5_system_t sys, struct vector_data *vec,
1217 real64 scalar, boolean transpose)
1218 {
1219 if( transpose ) { /* vec is indexed by col */
1220 int32 col;
1221 for( col=vec->rng->low; col<=vec->rng->high; col++ ) {
1222 sys->J.rhs[mtx_col_to_org(sys->J.mtx,col)] += scalar*vec->vec[col];
1223 }
1224 } else { /* vec is indexed by row */
1225 int32 row;
1226 for( row=vec->rng->low; row<=vec->rng->high; row++ ) {
1227 sys->J.rhs[mtx_row_to_org(sys->J.mtx,row)] += scalar*vec->vec[row];
1228 }
1229 }
1230 linsolqr_rhs_was_changed(sys->J.sys,sys->J.rhs);
1231 }
1232
1233 /*
1234 * Computes a step to solve the linearized equations.
1235 */
1236 static void calc_newton( slv5_system_t sys)
1237 {
1238 linsolqr_system_t lsys = sys->J.sys;
1239 int32 col;
1240
1241 if( sys->newton.accurate )
1242 return;
1243
1244 sys->J.rhs = linsolqr_get_rhs(lsys,0);
1245 mtx_zero_real64(sys->J.rhs,sys->cap);
1246 calc_rhs(sys, &(sys->residuals), -1.0, FALSE);
1247 linsolqr_solve(lsys,sys->J.rhs);
1248 col = sys->newton.rng->low;
1249 for( ; col <= sys->newton.rng->high; col++ ) {
1250 sys->newton.vec[col] =
1251 linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1252 }
1253 if (SAVLIN) {
1254 FILE *ldat;
1255 int32 ov;
1256 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1257 ldat=fopen(savlinfilename,"w");
1258 FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1259 sys->s.iteration);
1260 debug_write_array(ldat,sys->J.rhs,sys->cap);
1261 FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1262 for(ov=0 ; ov < sys->cap; ov++ )
1263 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1264 fclose(ldat);
1265 }
1266 square_norm( &(sys->newton) );
1267 sys->newton.accurate = TRUE;
1268 #if DEBUG
1269 FPRINTF(LIF(sys),"Newton: ");
1270 debug_out_vector(LIF(sys),&(sys->newton));
1271 #endif
1272 }
1273
1274 /*
1275 * Computes a step to solve the linearized equations.
1276 */
1277 static void calc_perturbed_newton( slv5_system_t sys)
1278 {
1279 linsolqr_system_t lsys = sys->J.sys;
1280 int32 col;
1281
1282 if( sys->perturbed_newton.accurate )
1283 return;
1284
1285 sys->J.rhs = linsolqr_get_rhs(lsys,1);
1286 mtx_zero_real64(sys->J.rhs,sys->cap);
1287 calc_rhs(sys, &(sys->perturbed_residuals), -1.0, FALSE);
1288 linsolqr_solve(lsys,sys->J.rhs);
1289 col = sys->perturbed_newton.rng->low;
1290 for( ; col <= sys->perturbed_newton.rng->high; col++ ) {
1291 sys->perturbed_newton.vec[col] =
1292 linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1293 }
1294 if (SAVLIN) {
1295 FILE *ldat;
1296 int32 ov;
1297 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1298 ldat=fopen(savlinfilename,"w");
1299 FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1300 sys->s.iteration);
1301 debug_write_array(ldat,sys->J.rhs,sys->cap);
1302 FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1303 for(ov=0 ; ov < sys->cap; ov++ )
1304 FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1305 fclose(ldat);
1306 }
1307 square_norm( &(sys->perturbed_newton) );
1308 sys->perturbed_newton.accurate = TRUE;
1309 #if DEBUG
1310 FPRINTF(LIF(sys),"Perturbed Newton: ");
1311 debug_out_vector(LIF(sys),&(sys->perturbed_newton));
1312 #endif
1313 }
1314
1315
1316
1317 /*
1318 * Calculate the perturbed descent direction
1319 * in the variables.
1320 */
1321 static void calc_varstep( slv5_system_t sys)
1322 {
1323 if( sys->varstep.accurate )
1324 return;
1325 copy_vector(&(sys->perturbed_newton),&(sys->varstep));
1326 sys->varstep.norm2 = sys->perturbed_newton.norm2;
1327
1328 sys->varstep.accurate = TRUE;
1329 #if DEBUG
1330 FPRINTF(LIF(sys),"Varstep: ");
1331 debug_out_vector(LIF(sys),&(sys->varstep));
1332 #endif
1333 }
1334
1335 /*
1336 * Calculate the hypothetical netown direction
1337 * in the variables.
1338 */
1339 static void calc_varnewstep( slv5_system_t sys)
1340 {
1341 if( sys->varnewstep.accurate )
1342 return;
1343 copy_vector(&(sys->newton),&(sys->varnewstep));
1344 sys->varnewstep.norm2 = sys->newton.norm2;
1345
1346 sys->varnewstep.accurate = TRUE;
1347 #if DEBUG
1348 FPRINTF(LIF(sys),"Varnewstep: ");
1349 debug_out_vector(LIF(sys),&(sys->varnewstep));
1350 #endif
1351 }
1352
1353 /*
1354 * Computes the error nu.
1355 */
1356 static void calc_nu( slv5_system_t sys)
1357 {
1358 struct rel_relation *rel;
1359 int32 row, r;
1360 real64 error;
1361
1362 error = 0.0;
1363 row = sys->residuals.rng->low;
1364 for(row = sys->residuals.rng->low; row <= sys->residuals.rng->high; row++) {
1365 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1366 if (!rel) {
1367 r=mtx_row_to_org(sys->J.mtx,row);
1368 FPRINTF(ASCERR,"NULL relation found !!\n");
1369 FPRINTF(ASCERR,"at row %d rel %d in calc_nu\n",(int)row,r);
1370 FFLUSH(ASCERR);
1371 }
1372 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
1373 #if DEBUG
1374 FPRINTF(ASCERR,"Complementary equation in calc_nu \n");
1375 FPRINTF(ASCERR,"row = %d\n",row);
1376 FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
1377 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1378 #endif /* DEBUG */
1379 error = error + rel_residual(rel);
1380 } else {
1381 if (rel_active(rel) && rel_included(rel)) {
1382 error = error + (rel_residual(rel) *rel_residual(rel) );
1383 #if DEBUG
1384 FPRINTF(ASCERR,"Non complementary equation calc_nu \n");
1385 FPRINTF(ASCERR,"row = %d\n",row);
1386 FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
1387 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1388 #endif /* DEBUG */
1389 }
1390 }
1391 }
1392
1393 sys->nu = error;
1394
1395 #if DEBUG_OBJ_VALUES
1396 FPRINTF(ASCERR," Error nu = %g \n",sys->nu);
1397 #endif /* DEBUG_OBJ_VALUES */
1398 }
1399
1400
1401 /*
1402 * Computes the objective function Psi.
1403 */
1404 static void calc_psi( slv5_system_t sys)
1405 {
1406
1407 struct rel_relation *rel;
1408 int32 row, r;
1409 real64 sum, sumt;
1410
1411 calc_nu( sys);
1412
1413 sum = 0.0;
1414 for(row = sys->residuals.rng->low; row <= sys->residuals.rng->high; row++) {
1415 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1416 if (!rel) {
1417 r = mtx_row_to_org(sys->J.mtx,row);
1418 FPRINTF(ASCERR,"NULL relation found !!\n");
1419 FPRINTF(ASCERR,"at row %d rel %d in calc_psi\n",(int)row,r);
1420 FFLUSH(ASCERR);
1421 }
1422 if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
1423 #if DEBUG
1424 FPRINTF(ASCERR,"Complementary equation in calc_psi\n");
1425 FPRINTF(ASCERR,"row = %d\n",row);
1426 FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
1427 FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1428 #endif /* DEBUG */
1429 sumt = log10(rel_residual(rel));
1430 sum = sum + sumt;
1431 }
1432 }
1433
1434 if (sys->eta > 0.0) {
1435 sys->psi = ( sys->eta * log10(sys->nu) ) - sum;
1436 } else {
1437 sys->psi = 0.5*sys->residuals.norm2;
1438 }
1439
1440 #if DEBUG_OBJ_VALUES
1441 FPRINTF(ASCERR," Psi = %g \n",sys->psi);
1442 #endif /* DEBUG_OBJ_VALUES */
1443 }
1444
1445
1446
1447 /*
1448 * Variable values maintenance
1449 * ---------------------------
1450 * restore_variables(sys)
1451 * coef = required_coef_to_stay_inbounds(sys)
1452 * apply_step(sys,coef)
1453 * step_accepted(sys)
1454 */
1455
1456 /*
1457 * Restores the values of the variables before applying
1458 * a step.
1459 */
1460 static void restore_variables( slv5_system_t sys)
1461 {
1462 int32 col;
1463 real64 *vec;
1464 vec = (sys->nominals.vec);
1465 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
1466 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1467 var_set_value(var,sys->variables.vec[col]*vec[col]);
1468 }
1469 }
1470
1471
1472 /*
1473 * Calculates the maximum fraction of the step which can be
1474 * taken without making negative the complementary vars.
1475 * It is assumed that the current complementary variable
1476 * is positive. The step must be calculated.
1477 */
1478 static real64 factor_for_complementary_vars( slv5_system_t sys, int32 v)
1479 {
1480 real64 factor, minfactor;
1481 struct var_variable *var;
1482 real64 dx,val,bnd;
1483 int32 col;
1484 struct vector_data step;
1485 real64 *vec;
1486
1487 vec = (sys->nominals.vec);
1488
1489 if (v == 1) {
1490 step = sys->varstep;
1491 } else {
1492 step = sys->varnewstep;
1493 }
1494
1495 minfactor = 1.0;
1496 factor = 1.0;
1497 for( col=step.rng->low; col <= step.rng->high; col++ ) {
1498 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1499 val = var_value(var);
1500 if (var_complementary(var) && (!var_fixed(var))) {
1501 dx = step.vec[col] * vec[col];
1502 bnd = 0.0;
1503 if( val + dx < bnd ) {
1504 factor = MIN((bnd-val)/dx, 1.0);
1505 if (factor < 1.0) {
1506 #if DEBUG_COMPLEMENTARY_VAR
1507 FPRINTF(ASCERR,"Negative Complementary Variable : \n");
1508 print_var_name(ASCERR,sys,var);
1509 FPRINTF(ASCERR,"\n");
1510 FPRINTF(ASCERR,"factor = %g \n",factor);
1511 #endif /* DEBUG_COMPLEMENTARY_VAR */
1512 }
1513 if( factor < minfactor ) {
1514 minfactor = factor;
1515 }
1516 }
1517 }
1518 }
1519 #if DEBUG_COMPLEMENTARY_VAR
1520 FPRINTF(ASCERR,"Minimum complementary factor = %g \n",minfactor);
1521 #endif /* DEBUG_COMPLEMENTARY_VAR */
1522 return minfactor;
1523 }
1524
1525
1526 /*
1527 * Calculates the fraction of the mehrotra step which can be
1528 * taken without making negative the complementary vars.
1529 * It is assumed that the current complementary variable
1530 * is positive. The step must be calculated.
1531 */
1532 static real64 quadratic_factor_for_complementary_vars( slv5_system_t sys)
1533 {
1534 struct var_variable *var;
1535 struct vector_data predictor,corrector;
1536 real64 *vec;
1537 real64 dx, dxp, dxc, val, bnd, try;
1538 real64 factor, minfactor, fup, flow;
1539 real64 error;
1540 int32 col;
1541 int32 conv_flag, iter;
1542
1543 vec = (sys->nominals.vec);
1544 predictor = sys->varnewstep;
1545 corrector = sys->varstep;
1546
1547 minfactor = 1.0;
1548 for( col = corrector.rng->low; col <= corrector.rng->high; col++ ) {
1549 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1550 val = var_value(var);
1551 if (var_complementary(var) && (!var_fixed(var))) {
1552 factor = 1.0;
1553 fup = 1.0;
1554 flow = 0.0;
1555 dxp = predictor.vec[col] * vec[col];
1556 dxc = corrector.vec[col] * vec[col];
1557 error = BIS_TOL * vec[col];
1558 dx = factor * dxp + pow(factor,2) * dxc;
1559 bnd = 0.0;
1560 if( val + dx < bnd ) {
1561 conv_flag = 0;
1562 while (conv_flag == 0) {
1563 iter = 0;
1564 factor = ( fup + flow) / 2.0;
1565 dx = factor * dxp + pow(factor,2) * dxc;
1566 try = val + dx;
1567 if (try < bnd) {
1568 fup = factor;
1569 } else {
1570 flow = factor;
1571 if ( (try - bnd) <= error ) {
1572 conv_flag = 1;
1573 }
1574 }
1575 iter++;
1576 if (iter > ITER_BIS_LIMIT) {
1577 FPRINTF(ASCERR,"quadratic_factor_for_complementary_vars \n");
1578 FPRINTF(ASCERR,
1579 "Reached Maximum number of iterations for bisection \n");
1580 FPRINTF(ASCERR,"Returning zero \n");
1581 factor = 0.0;
1582 return factor;
1583 }
1584 }
1585 }
1586 if (factor < 1.0) {
1587 #if DEBUG_COMPLEMENTARY_VAR
1588 FPRINTF(ASCERR,"Negative Complementary Variable : \n");
1589 print_var_name(ASCERR,sys,var);
1590 FPRINTF(ASCERR,"\n");
1591 FPRINTF(ASCERR,"quadratic factor = %g \n",factor);
1592 #endif /* DEBUG_COMPLEMENTARY_VAR */
1593 }
1594 if( factor < minfactor ) {
1595 minfactor = factor;
1596 }
1597 }
1598 }
1599 #if DEBUG_COMPLEMENTARY_VAR
1600 FPRINTF(ASCERR,"Complementary minimum quadratic factor = %g \n",minfactor);
1601 #endif /* DEBUG_COMPLEMENTARY_VAR */
1602 return minfactor;
1603 }
1604
1605
1606 /*
1607 * Adds step to the variable values in block. It projects
1608 * non complementary varaibles near bounds.
1609 */
1610 static void apply_quadratic_step( slv5_system_t sys, real64 factor)
1611 {
1612 FILE *lif = LIF(sys);
1613 struct var_variable *var;
1614 real64 dx, dxp, dxc, val, bnd;
1615 struct vector_data predictor, corrector;
1616 int32 col;
1617 real64 *vec;
1618
1619 vec = (sys->nominals.vec);
1620 predictor = sys->varnewstep;
1621 corrector = sys->varstep;
1622
1623 for(col = corrector.rng->low; col <= corrector.rng->high;col++) {
1624 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1625 dxp = vec[col] * predictor.vec[col];
1626 dxc = vec[col] * corrector.vec[col];
1627 val = var_value(var);
1628 dx = factor * dxp + pow(factor,2)*dxc;
1629 if( val + dx > (bnd=var_upper_bound(var)) ) {
1630 dx = TOWARD_BOUNDS*(bnd-val);
1631 if (SHOW_LESS_IMPT) {
1632 FPRINTF(lif,"%-40s ---> ",
1633 " Variable projected to upper bound");
1634 print_var_name(lif,sys,var); PUTC('\n',lif);
1635 }
1636 } else if( val + dx < (bnd=var_lower_bound(var)) ) {
1637 dx = TOWARD_BOUNDS*(bnd-val);
1638 if (SHOW_LESS_IMPT) {
1639 FPRINTF(lif,"%-40s ---> ",
1640 " Variable projected to lower bound");
1641 print_var_name(lif,sys,var); PUTC('\n',lif);
1642 }
1643 }
1644 var_set_value(var,val+dx);
1645 }
1646
1647 /* Allow weighted residuals to be recalculated at new point */
1648 sys->residuals.accurate = FALSE;
1649 sys->newton_residuals.accurate = FALSE;
1650 sys->perturbed_residuals.accurate = FALSE;
1651
1652 return;
1653 }
1654
1655
1656 /*
1657 * Adds step to the variable values in block. It projects
1658 * non complementary varaibles near bounds.
1659 */
1660 static void apply_step( slv5_system_t sys, int32 v, real64 factor)
1661 {
1662 FILE *lif = LIF(sys);
1663 struct var_variable *var;
1664 real64 dx,val,bnd;
1665 struct vector_data step;
1666 int32 col;
1667 real64 *vec;
1668
1669 vec = (sys->nominals.vec);
1670
1671 if (v == 1) {
1672 step = sys->varstep;
1673 } else {
1674 step = sys->varnewstep;
1675 }
1676
1677 for(col=step.rng->low; col<=step.rng->high;col++) {
1678 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1679 dx = vec[col]*step.vec[col];
1680 val = var_value(var);
1681 dx = dx*factor;
1682 if( val + dx > (bnd=var_upper_bound(var)) ) {
1683 dx = TOWARD_BOUNDS*(bnd-val);
1684 if (SHOW_LESS_IMPT) {
1685 FPRINTF(lif,"%-40s ---> ",
1686 " Variable projected to upper bound");
1687 print_var_name(lif,sys,var); PUTC('\n',lif);
1688 }
1689 } else if( val + dx < (bnd=var_lower_bound(var)) ) {
1690 dx = TOWARD_BOUNDS*(bnd-val);
1691 if (SHOW_LESS_IMPT) {
1692 FPRINTF(lif,"%-40s ---> ",
1693 " Variable projected to lower bound");
1694 print_var_name(lif,sys,var); PUTC('\n',lif);
1695 }
1696 }
1697 var_set_value(var,val+dx);
1698 }
1699
1700 /* Allow weighted residuals to be recalculated at new point */
1701 sys->residuals.accurate = FALSE;
1702 sys->newton_residuals.accurate = FALSE;
1703 sys->perturbed_residuals.accurate = FALSE;
1704
1705 return;
1706 }
1707
1708
1709 /*
1710 * Adds step to the variable values in block. It projects
1711 * non complementary varaibles near bounds.
1712 */
1713 static void apply_2nd_order_correction( slv5_system_t sys)
1714 {
1715 struct var_variable *var;
1716 real64 dx,val;
1717 struct vector_data step;
1718 int32 col;
1719 real64 *vec;
1720
1721 vec = (sys->nominals.vec);
1722 step = sys->varnewstep;
1723
1724 for(col=step.rng->low; col<=step.rng->high;col++) {
1725 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1726 if (var_active(var) && var_incident(var) && var_complementary(var)
1727 && (!var_fixed(var))) {
1728 dx = vec[col]*step.vec[col];
1729 val = dx;
1730 var_set_value(var,val);
1731 }
1732 }
1733
1734 /* Allow 2nd order correction to be recalculated */
1735 sys->correction.accurate = FALSE;
1736
1737 return;
1738 }
1739
1740
1741 /*
1742 * This function should be called when the step is accepted.
1743 */
1744 static void step_accepted( slv5_system_t sys)
1745 {
1746 /* Maintain update status on jacobian and weights */
1747 if (--(sys->update.jacobian) <= 0) {
1748 sys->J.accurate = FALSE;
1749 }
1750
1751 sys->variables.accurate = FALSE;
1752 sys->newton_residuals.accurate = FALSE;
1753 sys->perturbed_residuals.accurate = FALSE;
1754 sys->newton.accurate = FALSE;
1755 sys->correction.accurate = FALSE;
1756 sys->perturbed_newton.accurate = FALSE;
1757 sys->varnewstep.accurate = FALSE;
1758 sys->varstep.accurate = FALSE;
1759 }
1760
1761
1762 /*
1763 * Block routines
1764 * --------------
1765 * feas = block_feasible(sys)
1766 * move_to_next_block(sys)
1767 * find_next_unconverged_block(sys)
1768 */
1769
1770 /*
1771 * Returns TRUE if the current block is feasible, FALSE otherwise.
1772 * It is assumed that the residuals have been computed.
1773 */
1774 static boolean block_feasible( slv5_system_t sys)
1775 {
1776 int32 row;
1777
1778 if( !sys->s.calc_ok )
1779 return(FALSE);
1780
1781 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
1782 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1783 if( !rel_satisfied(rel) ) return FALSE;
1784 }
1785 return TRUE;
1786 }
1787
1788 /*
1789 * Moves on to the next block, updating all of the solver information.
1790 * To move to the first block, set sys->s.block.current_block to -1 before
1791 * calling. If already at the last block, then sys->s.block.current_block
1792 * will equal the number of blocks and the system will be declared
1793 * converged. Otherwise, the residuals for the new block will be computed
1794 * and sys->s.calc_ok set according.
1795 */
1796 static void move_to_next_block( slv5_system_t sys)
1797 {
1798 struct var_variable *var;
1799 struct rel_relation *rel;
1800 int32 row;
1801 int32 col;
1802 int32 ci;
1803
1804 if( sys->s.block.current_block >= 0 ) {
1805
1806 /* Record cost accounting info here. */
1807 ci=sys->s.block.current_block;
1808 sys->s.cost[ci].size = sys->s.block.current_size;
1809 sys->s.cost[ci].iterations = sys->s.block.iteration;
1810 sys->s.cost[ci].funcs = sys->s.block.funcs;
1811 sys->s.cost[ci].jacs = sys->s.block.jacs;
1812 sys->s.cost[ci].functime = sys->s.block.functime;
1813 sys->s.cost[ci].jactime = sys->s.block.jactime;
1814 sys->s.cost[ci].time = sys->s.block.cpu_elapsed;
1815 sys->s.cost[ci].resid = sys->s.block.residual;
1816
1817 /* De-initialize previous block */
1818 if (SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
1819 LIFDS)) {
1820 FPRINTF(LIF(sys),"Block %d converged.\n",
1821 sys->s.block.current_block);
1822 }
1823 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
1824 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1825 var_set_in_block(var,FALSE);
1826 }
1827 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
1828 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1829 rel_set_in_block(rel,FALSE);
1830 }
1831 sys->s.block.previous_total_size += sys->s.block.current_size;
1832 }
1833
1834 sys->s.block.current_block++;
1835 if( sys->s.block.current_block < sys->s.block.number_of ) {
1836 boolean ok;
1837
1838 /* Initialize next block */
1839
1840 sys->J.reg =
1841 (slv_get_solvers_blocks(SERVER))->block[sys->s.block.current_block];
1842
1843
1844 row = sys->J.reg.row.high - sys->J.reg.row.low + 1;
1845 col = sys->J.reg.col.high - sys->J.reg.col.low + 1;
1846 sys->s.block.current_size = MAX(row,col);
1847
1848 sys->s.block.iteration = 0;
1849 sys->s.block.cpu_elapsed = 0.0;
1850 sys->s.block.functime = 0.0;
1851 sys->s.block.jactime = 0.0;
1852 sys->s.block.funcs = 0;
1853 sys->s.block.jacs = 0;
1854
1855 if(SHOW_LESS_IMPT && (LIFDS ||
1856 sys->s.block.current_size > 1)) {
1857 debug_delimiter(LIF(sys));
1858 debug_delimiter(LIF(sys));
1859 }
1860 if(SHOW_LESS_IMPT && LIFDS) {
1861 FPRINTF(LIF(sys),"\n%-40s ---> %d in [%d..%d]\n",
1862 "Current block number", sys->s.block.current_block,
1863 0, sys->s.block.number_of-1);
1864 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Current block size",
1865 sys->s.block.current_size);
1866 }
1867 sys->s.calc_ok = TRUE;
1868
1869 if (!(sys->p.ignore_bounds) ) {
1870 slv_insure_bounds(SERVER, sys->J.reg.col.low,
1871 sys->J.reg.col.high,MIF(sys));
1872 }
1873
1874 sys->residuals.accurate = FALSE;
1875 if( !(ok = calc_residuals(sys)) ) {
1876 FPRINTF(MIF(sys),
1877 "Residual calculation errors detected in move_to_next_block.\n");
1878 }
1879
1880 /*
1881 * Update number of complementary equations
1882 */
1883 calc_comp(sys);
1884 calc_eta(sys);
1885
1886 if( SHOW_LESS_IMPT &&
1887 (sys->s.block.current_size >1 ||
1888 LIFDS) ) {
1889 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
1890 sys->s.block.residual);
1891 }
1892 sys->s.calc_ok = sys->s.calc_ok && ok;
1893
1894 /* Must be updated as soon as required */
1895 sys->J.accurate = FALSE;
1896 sys->update.weights = 0;
1897 sys->update.nominals = 0;
1898 sys->variables.accurate = FALSE;
1899 sys->newton_residuals.accurate = FALSE;
1900 sys->perturbed_residuals.accurate = FALSE;
1901 sys->newton.accurate = FALSE;
1902 sys->correction.accurate = FALSE;
1903 sys->perturbed_newton.accurate = FALSE;
1904 sys->varnewstep.accurate = FALSE;
1905 sys->varstep.accurate = FALSE;
1906
1907 } else {
1908 boolean ok;
1909 /*
1910 * Before we claim convergence, we must check if we left behind
1911 * some unassigned relations. If and only if they happen to be
1912 * satisfied at the current point, convergence has been obtained.
1913 *
1914 * Also insures that all included relations have valid residuals.
1915 * Included inequalities will have correct residuals.
1916 * Unsatisfied included inequalities cause inconsistency.
1917 *
1918 * This of course ignores that fact an objective function might
1919 * be present. Then, feasibility isn't enough, is it now.
1920 */
1921 if( sys->s.struct_singular ) {
1922 /* black box w/singletons provoking bug here, maybe */
1923 sys->s.block.current_size = sys->rused - sys->rank;
1924 if(SHOW_LESS_IMPT) {
1925 debug_delimiter(LIF(sys));
1926 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Unassigned Relations",
1927 sys->s.block.current_size);
1928 }
1929 sys->J.reg.row.low = sys->J.reg.col.low = sys->rank;
1930 sys->J.reg.row.high = sys->J.reg.col.high = sys->rused - 1;
1931 sys->residuals.accurate = FALSE;
1932 if( !(ok=calc_residuals(sys)) ) {
1933 FPRINTF(MIF(sys),
1934 "Residual calculation errors detected in leftover equations.\n");
1935 }
1936 if(SHOW_LESS_IMPT) {
1937 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
1938 sys->s.block.residual);
1939 }
1940 if( block_feasible(sys) ) {
1941 if(SHOW_LESS_IMPT) {
1942 FPRINTF(LIF(sys),"\nUnassigned relations ok. You lucked out.\n");
1943 }
1944 sys->s.converged = TRUE;
1945 } else {
1946 if(SHOW_LESS_IMPT) {
1947 FPRINTF(LIF(sys),"\nProblem inconsistent: %s.\n",
1948 "Unassigned relations not satisfied");
1949 }
1950 sys->s.inconsistent = TRUE;
1951 }
1952 if(SHOW_LESS_IMPT) {
1953 debug_delimiter(LIF(sys));
1954 }
1955 } else {
1956 sys->s.converged = TRUE;
1957 }
1958 }
1959 }
1960
1961
1962 /*
1963 * Calls the appropriate reorder function on a block
1964 */
1965 static void reorder_new_block(slv5_system_t sys)
1966 {
1967 int32 method;
1968 if( sys->s.block.current_block < sys->s.block.number_of ) {
1969 if (strcmp(REORDER_OPTION,"SPK1") == 0) {
1970 method = 2;
1971 } else {
1972 method = 1;
1973 }
1974
1975 if( sys->s.block.current_block <= sys->s.block.current_reordered_block &&
1976 sys->s.cost[sys->s.block.current_block].reorder_method == method &&
1977 sys->s.block.current_block >= 0 ) {
1978 #if DEBUG
1979 FPRINTF(ASCERR,"YOU JUST AVOIDED A REORDERING\n");
1980 #endif
1981 slv_set_up_block(SERVER,sys->s.block.current_block);
1982 /* tell linsol to bless it and get on with things */
1983 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
1984 return; /*must have been reordered since last system build*/
1985 }
1986
1987 /* Let the slv client function take care of reordering things
1988 * and setting in block flags.
1989 */
1990 if (strcmp(REORDER_OPTION,"SPK1") == 0) {
1991 sys->s.cost[sys->s.block.current_block].reorder_method = 2;
1992 slv_spk1_reorder_block(SERVER,sys->s.block.current_block,1);
1993 } else if (strcmp(REORDER_OPTION,"TEAR_DROP") == 0) {
1994 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
1995 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
1996 CUTOFF,
1997 0,mtx_SPK1);
1998 /* khack: try tspk1 for transpose case */
1999 } else if (strcmp(REORDER_OPTION,"OVER_TEAR") == 0) {
2000 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2001 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2002 CUTOFF,
2003 1,mtx_SPK1);
2004 } else {
2005 sys->s.cost[sys->s.block.current_block].reorder_method = 1;
2006 FPRINTF(MIF(sys),"IPSlv called with unknown reorder option\n");
2007 FPRINTF(MIF(sys),"IPSlv using single edge tear drop (TEAR_DROP).\n");
2008 slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
2009 CUTOFF,0,mtx_SPK1);
2010 }
2011 /* tell linsol to bless it and get on with things */
2012 linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
2013 if (sys->s.block.current_block > sys->s.block.current_reordered_block) {
2014 sys->s.block.current_reordered_block = sys->s.block.current_block;
2015 }
2016 }
2017 }
2018
2019 /*
2020 * Moves to next unconverged block, assuming that the current block has
2021 * converged (or is -1, to start).
2022 */
2023 static void find_next_unconverged_block( slv5_system_t sys)
2024 {
2025 do {
2026 move_to_next_block(sys);
2027 #if DEBUG
2028 debug_out_var_values(ASCERR,sys);
2029 debug_out_rel_residuals(ASCERR,sys);
2030 #endif
2031 } while( !sys->s.converged && block_feasible(sys));
2032 reorder_new_block(sys);
2033 }
2034
2035
2036 /*
2037 * Iteration begin/end routines
2038 * ----------------------------
2039 * iteration_begins(sys)
2040 * iteration_ends(sys)
2041 */
2042
2043 /*
2044 * Prepares sys for entering an iteration, increasing the iteration counts
2045 * and starting the clock.
2046 */
2047 static void iteration_begins( slv5_system_t sys)
2048 {
2049 sys->clock = tm_cpu_time();
2050 ++(sys->s.block.iteration);
2051 ++(sys->s.iteration);
2052 if(SHOW_LESS_IMPT&& (sys->s.block.current_size >1 ||
2053 LIFDS)) {
2054 FPRINTF(LIF(sys),"\n%-40s ---> %d\n",
2055 "Iteration", sys->s.block.iteration);
2056 FPRINTF(LIF(sys),"%-40s ---> %d\n",
2057 "Total iteration", sys->s.iteration);
2058 }
2059 }
2060
2061 /*
2062 * Prepares sys for exiting an iteration, stopping the clock and recording
2063 * the cpu time.
2064 */
2065 static void iteration_ends( slv5_system_t sys)
2066 {
2067 double cpu_elapsed; /* elapsed this iteration */
2068
2069 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
2070 sys->s.block.cpu_elapsed += cpu_elapsed;
2071 sys->s.cpu_elapsed += cpu_elapsed;
2072 if(SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
2073 LIFDS)) {
2074 FPRINTF(LIF(sys),"%-40s ---> %g\n",
2075 "Elapsed time", sys->s.block.cpu_elapsed);
2076 FPRINTF(LIF(sys),"%-40s ---> %g\n",
2077 "Total elapsed time", sys->s.cpu_elapsed);
2078 }
2079 }
2080
2081 /*
2082 * Updates the solver status.
2083 */
2084 static void update_status( slv5_system_t sys)
2085 {
2086 boolean unsuccessful;
2087
2088 if( !sys->s.converged ) {
2089 sys->s.time_limit_exceeded =
2090 (sys->s.block.cpu_elapsed >= TIME_LIMIT);
2091 sys->s.iteration_limit_exceeded =
2092 (sys->s.block.iteration >= ITER_LIMIT);
2093 }
2094
2095 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
2096 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
2097
2098 sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
2099 sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
2100 }
2101
2102 static
2103 int32 slv5_get_default_parameters(slv_system_t server, SlvClientToken asys,
2104 slv_parameters_t *parameters)
2105 {
2106 slv5_system_t sys;
2107 union parm_arg lo,hi,val;
2108 struct slv_parameter *new_parms = NULL;
2109 int32 make_macros = 0;
2110
2111 static char *factor_names[] = {
2112 "SPK1/RANKI","SPK1/RANKI+ROW",
2113 "Fast-SPK1/RANKI","Fast-SPK1/RANKI+ROW",
2114 "Fastest-SPK1/MR-RANKI","CondQR","CPQR"
2115 /* ,"GAUSS","GAUSS_EASY" currently only works for ken */
2116 };
2117 static char *reorder_names[] = {
2118 "SPK1","TEAR_DROP","OVER_TEAR"
2119 };
2120 static char *converge_names[] = {
2121 "ABSOLUTE","RELNOM_SCALE"
2122 };
2123 static char *scaling_names[] = {
2124 "NONE","ROW_2NORM"
2125 };
2126
2127 if (server != NULL && asys != NULL) {
2128 sys = SLV5(asys);
2129 make_macros = 1;
2130 }
2131
2132 #ifndef NDEBUG /* keep purify from whining on UMR */
2133 lo.argr = hi.argr = val.argr = 0.0;
2134 #endif
2135
2136 if (parameters->parms == NULL) {
2137 /* an external client wants our parameter list.
2138 * an instance of slv5_system_structure has this pointer
2139 * already set in slv5_create
2140 */
2141 new_parms = (struct slv_parameter *)
2142 ascmalloc((slv5_PA_SIZE)*sizeof(struct slv_parameter));
2143 if (new_parms == NULL) {
2144 return -1;
2145 }
2146 parameters->parms = new_parms;
2147 parameters->dynamic_parms = 1;
2148 }
2149 parameters->num_parms = 0;
2150
2151 /* begin defining parameters */
2152
2153 slv_define_parm(parameters, bool_parm,
2154 "ignorebounds","ignore bounds?","ignore bounds?",
2155 U_p_bool(val,0),U_p_bool(lo,0),U_p_bool(hi,1),-1);
2156 SLV_BPARM_MACRO(IGNORE_BOUNDS_PTR,parameters);
2157
2158 slv_define_parm(parameters, real_parm,
2159 "rho", "search parameter", "search parameter",
2160 U_p_real(val,0.5),U_p_real(lo, 0),U_p_real(hi,1.0), 1);
2161 SLV_RPARM_MACRO(RHO_PTR,parameters);
2162
2163 slv_define_parm(parameters, bool_parm,
2164 "partition", "partitioning enabled", "partitioning enabled",
2165 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2166 SLV_BPARM_MACRO(PARTITION_PTR,parameters);
2167
2168 slv_define_parm(parameters, bool_parm,
2169 "showlessimportant", "detailed solving info",
2170 "detailed solving info",
2171 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2172 SLV_BPARM_MACRO(SHOW_LESS_IMPT_PTR,parameters);
2173
2174 slv_define_parm(parameters, int_parm,
2175 "timelimit", "time limit (CPU sec/block)",
2176 "time limit (CPU sec/block)",
2177 U_p_int(val,1500),U_p_int(lo, 1),U_p_int(hi,20000),1);
2178 SLV_IPARM_MACRO(TIME_LIMIT_PTR,parameters);
2179
2180 slv_define_parm(parameters, int_parm,
2181 "iterationlimit", "max iterations/block",
2182 "max iterations/block",
2183 U_p_int(val, 30),U_p_int(lo, 1),U_p_int(hi,20000),1);
2184 SLV_IPARM_MACRO(ITER_LIMIT_PTR,parameters);
2185
2186 slv_define_parm(parameters, int_parm,
2187 "bisiterationlimit", "max iterations for bisection",
2188 "max iterations for bisection",
2189 U_p_int(val, 50),U_p_int(lo, 1),U_p_int(hi,20000),1);
2190 SLV_IPARM_MACRO(ITER_BIS_LIMIT_PTR,parameters);
2191
2192 slv_define_parm(parameters,real_parm,
2193 "bistol","tolerance for bisection" ,
2194 "tolerance for bisection" ,
2195 U_p_real(val,1e-08),U_p_real(lo,0),U_p_real(hi,1.0), 1);
2196 SLV_RPARM_MACRO(BIS_TOL_PTR,parameters);
2197
2198 slv_define_parm(parameters,real_parm,
2199 "alpha","search coefficient" ,"search coefficient" ,
2200 U_p_real(val,0.5),U_p_real(lo,0),U_p_real(hi,1.0), 1);
2201 SLV_RPARM_MACRO(ALPHA_PTR,parameters);
2202
2203 slv_define_parm(parameters, real_parm,
2204 "singtol", "epsilon (min pivot)", "epsilon (min pivot)",
2205 U_p_real(val, 1e-12),U_p_real(lo, 1e-12),U_p_real(hi,1.0),1);
2206 SLV_RPARM_MACRO(SING_TOL_PTR,parameters);
2207
2208 slv_define_parm(parameters, real_parm,
2209 "pivottol", "condition tolerance", "condition tolerance",
2210 U_p_real(val, 0.5),U_p_real(lo, 0),U_p_real(hi, 1),1);
2211 SLV_RPARM_MACRO(PIVOT_TOL_PTR,parameters);
2212
2213 slv_define_parm(parameters, real_parm,
2214 "feastol", "max residual (absolute)",
2215 "max residual (absolute)",
2216 U_p_real(val, 1e-8),U_p_real(lo, 1e-13),
2217 U_p_real(hi,100000.5),1);
2218 SLV_RPARM_MACRO(FEAS_TOL_PTR,parameters);
2219
2220 slv_define_parm(parameters, bool_parm,
2221 "lifds", "show singletons details", IEX(0),
2222 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2223 SLV_BPARM_MACRO(LIFDS_PTR,parameters);
2224
2225 slv_define_parm(parameters, bool_parm,
2226 "savlin", "write to file SlvLinsol.dat", IEX(1),
2227 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2228 SLV_BPARM_MACRO(SAVLIN_PTR,parameters);
2229
2230 slv_define_parm(parameters, bool_parm,
2231 "safe_calc", "safe calculations", IEX(8),
2232 U_p_bool(val, 1),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2233
2234 SLV_BPARM_MACRO(SAFE_CALC_PTR,parameters);
2235
2236
2237 slv_define_parm(parameters, int_parm,
2238 "cutoff", "block size cutoff (MODEL-based)", IEX(2),
2239 U_p_int(val, 500),U_p_int(lo,0),U_p_int(hi,20000), 2);
2240 SLV_IPARM_MACRO(CUTOFF_PTR,parameters);
2241
2242
2243 slv_define_parm(parameters, int_parm,
2244 "upjac", "Jacobian update frequency", IEX(3),
2245 U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
2246 SLV_IPARM_MACRO(UPDATE_JACOBIAN_PTR,parameters);
2247
2248 slv_define_parm(parameters, int_parm,
2249 "upwts", "Row scaling update frequency", IEX(4),
2250 U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
2251 SLV_IPARM_MACRO(UPDATE_WEIGHTS_PTR,parameters);
2252
2253 slv_define_parm(parameters, int_parm,
2254 "upnom", "Column scaling update frequency", IEX(5),
2255 U_p_int(val, 1000),U_p_int(lo,0),U_p_int(hi,20000), 3);
2256 SLV_IPARM_MACRO(UPDATE_NOMINALS_PTR,parameters);
2257
2258 slv_define_parm(parameters, char_parm,
2259 "convopt", "convergence test", "convergence test",
2260 U_p_string(val,converge_names[0]),
2261 U_p_strings(lo,converge_names),
2262 U_p_int(hi,sizeof(converge_names)/sizeof(char *)),1);
2263 SLV_CPARM_MACRO(CONVOPT_PTR,parameters);
2264
2265 slv_define_parm(parameters, char_parm,
2266 "scaleopt", "jacobian scaling option", IEX(9),
2267 U_p_string(val,scaling_names[1]),
2268 U_p_strings(lo,scaling_names),
2269 U_p_int(hi,sizeof(scaling_names)/sizeof(char *)),1);
2270 SLV_CPARM_MACRO(SCALEOPT_PTR,parameters);
2271
2272 slv_define_parm(parameters, bool_parm,
2273 "cncols", "Check poorly scaled columns", IEX(6),
2274 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2275 SLV_BPARM_MACRO(DUMPCNORM_PTR,parameters);
2276
2277 slv_define_parm(parameters, bool_parm,
2278 "lintime", "Enable linsolqr timing",
2279 "Enable linsolqr factor timing",
2280 U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2281 SLV_BPARM_MACRO(LINTIME_PTR,parameters);
2282
2283
2284 slv_define_parm(parameters, char_parm,
2285 "reorder", "reorder method", IEX(7),
2286 U_p_string(val,reorder_names[0]),
2287 U_p_strings(lo,reorder_names),
2288 U_p_int(hi,sizeof(reorder_names)/sizeof(char *)),1);
2289 SLV_CPARM_MACRO(REORDER_OPTION_PTR,parameters);
2290
2291
2292 slv_define_parm(parameters, real_parm,
2293 "toosmall", "default for zero nominal", REX(0),
2294 U_p_real(val, 1e-8),U_p_real(lo, 1e-12),U_p_real(hi,1.5), 3);
2295 SLV_RPARM_MACRO(TOO_SMALL_PTR,parameters);
2296
2297 slv_define_parm(parameters, real_parm,
2298 "cnlow", "smallest allowable column norm", REX(1),
2299 U_p_real(val, 0.01),U_p_real(lo, 0),
2300 U_p_real(hi,100000000.5), 3);
2301 SLV_RPARM_MACRO(CNLOW_PTR,parameters);
2302
2303 slv_define_parm(parameters, real_parm,
2304 "cnhigh", "largest allowable column norm", REX(2),
2305 U_p_real(val, 100.0),U_p_real(lo,0),
2306 U_p_real(hi,100000000.5), 3);
2307 SLV_RPARM_MACRO(CNHIGH_PTR,parameters);
2308
2309 slv_define_parm(parameters, real_parm,
2310 "tobnds", "fraction move to bounds", REX(3),
2311 U_p_real(val, 0.995),U_p_real(lo, 0),U_p_real(hi,1.0), 3);
2312 SLV_RPARM_MACRO(TOWARD_BOUNDS_PTR,parameters);
2313
2314 slv_define_parm(parameters, char_parm,
2315 "bppivoting","linear method","linear method choice",
2316 U_p_string(val,factor_names[4]),
2317 U_p_strings(lo,factor_names),
2318 U_p_int(hi,sizeof(factor_names)/sizeof(char *)),1);
2319 SLV_CPARM_MACRO(FACTOR_OPTION_PTR,parameters);
2320
2321 return 1;
2322 }
2323
2324 /*
2325 * External routines
2326 * -----------------
2327 * See slv_client.h
2328 */
2329
2330 static SlvClientToken slv5_create(slv_system_t server, int *statusindex)
2331 {
2332 slv5_system_t sys;
2333
2334 sys = (slv5_system_t)asccalloc(1, sizeof(struct slv5_system_structure) );
2335 if (sys==NULL) {
2336 *statusindex = 1;
2337 return sys;
2338 }
2339 SERVER = server;
2340 sys->p.parms = sys->pa;
2341 sys->p.dynamic_parms = 0;
2342 slv5_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
2343 sys->integrity = OK;
2344 sys->presolved = 0;
2345 sys->p.output.more_important = stdout;
2346 sys->p.output.less_important = stdout;
2347 sys->J.old_partition = TRUE;
2348 sys->p.whose = (*statusindex);
2349 sys->s.ok = TRUE;
2350 sys->s.calc_ok = TRUE;
2351 sys->s.costsize = 0;
2352 sys->s.cost = NULL; /*redundant, but sanity preserving */
2353 sys->vlist = slv_get_solvers_var_list(server);
2354 sys->rlist = slv_get_solvers_rel_list(server);
2355 sys->obj = slv_get_obj_relation(server);
2356 if (sys->vlist == NULL) {
2357 ascfree(sys);
2358 FPRINTF(ASCERR,"IPSlv called with no variables.\n");
2359 *statusindex = -2;
2360 return NULL; /*_prolly leak here */
2361 }
2362 if (sys->rlist == NULL && sys->obj == NULL) {
2363 ascfree(sys);
2364 FPRINTF(ASCERR,"IPSlv called with no relations or objective.\n");
2365 *statusindex = -1;
2366 return NULL; /*_prolly leak here */
2367 }
2368 /* we don't give a damn about the objective list or the pars or
2369 * bounds or extrels or any of the other crap.
2370 */
2371 slv_check_var_initialization(server);
2372 *statusindex = 0;
2373 return((SlvClientToken)sys);
2374
2375 }
2376
2377 static void destroy_matrices( slv5_system_t sys)
2378 {
2379 if( sys->J.sys ) {
2380 int count = linsolqr_number_of_rhs(sys->J.sys)-1;
2381 for( ; count >= 0; count-- ) {
2382 destroy_array(linsolqr_get_rhs(sys->J.sys,count));
2383 }
2384 mtx_destroy(linsolqr_get_matrix(sys->J.sys));
2385 linsolqr_set_matrix(sys->J.sys,NULL);
2386 linsolqr_destroy(sys->J.sys);
2387 sys->J.sys = NULL;
2388 }
2389 }
2390
2391 static void destroy_vectors( slv5_system_t sys)
2392 {
2393 destroy_array(sys->nominals.vec);
2394 destroy_array(sys->weights.vec);
2395 destroy_array(sys->variables.vec);
2396 destroy_array(sys->residuals.vec);
2397 destroy_array(sys->newton_residuals.vec);
2398 destroy_array(sys->perturbed_residuals.vec);
2399 destroy_array(sys->newton.vec);
2400 destroy_array(sys->correction.vec);
2401 destroy_array(sys->perturbed_newton.vec);
2402 destroy_array(sys->varnewstep.vec);
2403 destroy_array(sys->varstep.vec);
2404 }
2405
2406 static int slv5_eligible_solver(slv_system_t server)
2407 {
2408 struct rel_relation **rp;
2409 rel_filter_t rfilter;
2410
2411 rfilter.matchbits = (REL_INCLUDED | REL_ACTIVE);
2412 rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
2413
2414 if (!slv_count_solvers_rels(server,&rfilter)) {
2415 return (FALSE);
2416 }
2417
2418 for( rp=slv_get_solvers_rel_list(server); *rp != NULL ; ++rp ) {
2419 if( rel_less(*rp) || rel_greater(*rp) ) return(FALSE);
2420 }
2421 return(TRUE);
2422 }
2423
2424 static
2425 void slv5_get_parameters(slv_system_t server, SlvClientToken asys,
2426 slv_parameters_t *parameters)
2427 {
2428 slv5_system_t sys;
2429 (void) server;
2430 sys = SLV5(asys);
2431 if (check_system(sys)) return;
2432 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
2433 }
2434
2435 static void slv5_set_parameters(slv_system_t server, SlvClientToken asys,
2436 slv_parameters_t *parameters)
2437 {
2438 slv5_system_t sys;
2439 (void) server;
2440 sys = SLV5(asys);
2441 if (check_system(sys)) return;
2442 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
2443 }
2444
2445 static void slv5_get_status(slv_system_t server, SlvClientToken asys,
2446 slv_status_t *status)
2447 {
2448 slv5_system_t sys;
2449 (void) server;
2450 sys = SLV5(asys);
2451 if (check_system(sys)) return;
2452 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
2453 }
2454
2455 static linsolqr_system_t slv5_get_linsolqr_sys(slv_system_t server,
2456 SlvClientToken asys)
2457 {
2458 slv5_system_t sys;
2459 (void) server;
2460 sys = SLV5(asys);
2461 if (check_system(sys)) return NULL;
2462 return(sys->J.sys);
2463 }
2464
2465 static linsol_system_t slv5_get_linsol_sys(slv_system_t server,
2466 SlvClientToken asys)
2467 {
2468 (void) server;
2469 (void) asys;
2470 return( NULL );
2471 }
2472
2473 /*
2474 * Performs structural analysis on the system, setting the flags in
2475 * status. The problem must be set up, the relation/variable list
2476 * must be non-NULL. The
2477 * jacobian (linear) system must be created and have the correct order
2478 * (stored in sys->cap). Everything else will be determined here.
2479 * On entry there isn't yet a correspondence between var_sindex and
2480 * jacobian column. Here we establish that.
2481 */
2482 static void structural_analysis(slv_system_t server, slv5_system_t sys)
2483 {
2484 var_filter_t vfilter;
2485 rel_filter_t rfilter;
2486
2487 /*
2488 * The server has marked incidence flags already.
2489 */
2490 /* count included equalities */
2491 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
2492 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
2493 sys->rused = slv_count_solvers_rels(server,&rfilter);
2494
2495 /* count free and incident vars */
2496 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
2497 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
2498 sys->vused = slv_count_solvers_vars(server,&vfilter);
2499
2500 /* Symbolic analysis */
2501 sys->rtot = slv_get_num_solvers_rels(server);
2502 sys->vtot = slv_get_num_solvers_vars(server);
2503 if (sys->rtot) {
2504 slv_block_partition(server);
2505 }
2506 sys->J.dofdata = slv_get_dofdata(server);
2507 sys->rank = sys->J.dofdata->structural_rank;
2508
2509 if( !(PARTITION) ) {
2510 /* maybe we should reorder blocks here? maybe not */
2511 slv_block_unify(server);
2512 }
2513
2514 slv_check_bounds(SERVER,sys->vused,sys->vtot-1,MIF(sys),"fixed");
2515
2516 /* Initialize Status */
2517 sys->s.over_defined = (sys->rused > sys->vused);
2518 sys->s.under_defined = (sys->rused < sys->vused);
2519 sys->s.struct_singular = (sys->rank < sys->rused);
2520 sys->s.block.number_of = (slv_get_solvers_blocks(SERVER))->nblocks;
2521 }
2522
2523
2524
2525 static void set_factor_options (slv5_system_t sys)
2526 {
2527 if (strcmp(FACTOR_OPTION,"SPK1/RANKI") == 0) {
2528 sys->J.fm = ranki_kw;
2529 } else if (strcmp(FACTOR_OPTION,"SPK1/RANKI+ROW") == 0) {
2530 sys->J.fm = ranki_jz;
2531 } else if (strcmp(FACTOR_OPTION,"Fast-SPK1/RANKI") == 0) {
2532 sys->J.fm = ranki_kw2;
2533 } else if (strcmp(FACTOR_OPTION,"Fast-SPK1/RANKI+ROW") == 0) {
2534 sys->J.fm = ranki_jz2;
2535 } else if (strcmp(FACTOR_OPTION,"Fastest-SPK1/MR-RANKI") == 0) {
2536 sys->J.fm = ranki_ba2;
2537 } else {
2538 sys->J.fm = ranki_ba2;
2539 }
2540 mtx_set_order(sys->J.mtx,sys->cap);
2541
2542 linsolqr_set_matrix(sys->J.sys,sys->J.mtx);
2543 linsolqr_prep(sys->J.sys,linsolqr_fmethod_to_fclass(sys->J.fm));
2544 linsolqr_set_pivot_zero(sys->J.sys, SING_TOL);
2545
2546 /* KHACK NOTE: looks like drop tol never set on interface */
2547 linsolqr_set_drop_tolerance(sys->J.sys, sys->p.tolerance.drop);
2548 linsolqr_set_pivot_tolerance(sys->J.sys, PIVOT_TOL);
2549 /* this next one is fishy, but we don't use qr so not panicking */
2550 linsolqr_set_condition_tolerance(sys->J.sys, PIVOT_TOL);
2551 }
2552
2553
2554 /*
2555 configures linsolqr system, among other things.
2556 sets type to be ranki_kw or ranki_jz as determined by
2557 parameters.
2558 */
2559 static void create_matrices(slv_system_t server, slv5_system_t sys)
2560 {
2561 sys->J.sys = linsolqr_create();
2562 sys->J.mtx = mtx_create();
2563
2564 set_factor_options(sys);
2565
2566 /* rhs 0 for sys->newton */
2567 sys->J.rhs = create_array(sys->cap,real64);
2568 linsolqr_add_rhs(sys->J.sys,sys->J.rhs,FALSE);
2569 /* rhs 1 for sys->perturbed_newton */
2570 sys->J.rhs = create_array(sys->cap,real64);
2571 linsolqr_add_rhs(sys->J.sys,sys->J.rhs,FALSE);
2572 structural_analysis(server,sys);
2573
2574 }
2575
2576 static void create_vectors(sys)
2577 slv5_system_t sys;
2578 {
2579 sys->nominals.vec = create_array(sys->cap,real64);
2580 sys->nominals.rng = &(sys->J.reg.col);
2581 sys->weights.vec = create_array(sys->cap,real64);
2582 sys->weights.rng = &(sys->J.reg.row);
2583 sys->variables.vec = create_array(sys->cap,real64);
2584 sys->variables.rng = &(sys->J.reg.col);
2585 sys->residuals.vec = create_array(sys->cap,real64);
2586 sys->residuals.rng = &(sys->J.reg.row);
2587 sys->newton_residuals.vec = create_array(sys->cap,real64);
2588 sys->newton_residuals.rng = &(sys->J.reg.row);
2589 sys->perturbed_residuals.vec = create_array(sys->cap,real64);
2590 sys->perturbed_residuals.rng = &(sys->J.reg.row);
2591 sys->newton.vec = create_array(sys->cap,real64);
2592 sys->newton.rng = &(sys->J.reg.col);
2593 sys->correction.vec = create_array(sys->cap,real64);
2594 sys->correction.rng = &(sys->J.reg.col);
2595 sys->perturbed_newton.vec = create_array(sys->cap,real64);
2596 sys->perturbed_newton.rng = &(sys->J.reg.col);
2597 sys->varnewstep.vec = create_array(sys->cap,real64);
2598 sys->varnewstep.rng = &(sys->J.reg.col);
2599 sys->varstep.vec = create_array(sys->cap,real64);
2600 sys->varstep.rng = &(sys->J.reg.col);
2601
2602 }
2603
2604 static
2605 void slv5_dump_internals(slv_system_t server, SlvClientToken sys,int level)
2606 {
2607 (void) server;
2608 check_system(sys);
2609 if (level > 0) {
2610 FPRINTF(ASCERR,"ERROR: (slv5) slv5_dump_internals\n");
2611 FPRINTF(ASCERR," slv5 does not dump its internals.\n");
2612 }
2613 }
2614
2615 /*
2616 * Here we will check if any fixed or included flags have
2617 * changed since the last presolve.
2618 */
2619 static
2620 int32 slv5_dof_changed(slv5_system_t sys)
2621 {
2622 int32 ind, result = 0;
2623 /* Currently we have two copies of the fixed and included flags
2624 which must be kept in sync. The var_fixed and rel_included
2625 functions perform the syncronization and hence must be called
2626 over the whole var list and rel list respectively. When we move
2627 to using only one set of flags (bit flags) this function can
2628 be changed to return 1 at the first indication of a change
2629 in the dof. */
2630
2631 /* search for vars that were fixed and are now free */
2632 for( ind = sys->vused; ind < sys->vtot; ++ind ) {
2633 if( !var_fixed(sys->vlist[ind]) && var_active(sys->vlist[ind]) ) {
2634 ++result;
2635 }
2636 }
2637 /* search for rels that were unincluded and are now included */
2638 for( ind = sys->rused; ind < sys->rtot; ++ind ) {
2639 if( rel_included(sys->rlist[ind]) && rel_active(sys->rlist[ind])) {
2640 ++result;
2641 }
2642 }
2643 /* search for vars that were free and are now fixed */
2644 for( ind = sys->vused -1; ind >= 0; --ind ) {
2645 if( var_fixed(sys->vlist[ind]) || !var_active(sys->vlist[ind])) {
2646 ++result;
2647 }
2648 }
2649 /* search for rels that were included and are now unincluded */
2650 for( ind = sys->rused -1; ind >= 0; --ind ) {
2651 if( !rel_included(sys->rlist[ind]) || !rel_active(sys->rlist[ind]) ) {
2652 ++result;
2653 }
2654 }
2655 return result;
2656 }
2657
2658
2659 static void reset_cost(struct slv_block_cost *cost,int32 costsize)
2660 {
2661 int32 ind;
2662 for( ind = 0; ind < costsize; ++ind ) {
2663 cost[ind].size = 0;
2664 cost[ind].iterations = 0;
2665 cost[ind].funcs = 0;
2666 cost[ind].jacs = 0;
2667 cost[ind].functime = 0;
2668 cost[ind].jactime = 0;
2669 cost[ind].time = 0;
2670 cost[ind].resid = 0;
2671 }
2672 }
2673
2674 static void slv5_update_linsolqr(slv5_system_t sys)
2675 {
2676 if (strcmp(FACTOR_OPTION,"SPK1/RANKI") == 0) {
2677 sys->J.fm = ranki_kw;
2678 } else if (strcmp(FACTOR_OPTION,"SPK1/RANKI+ROW") == 0) {
2679 sys->J.fm = ranki_jz;
2680 } else if (strcmp(FACTOR_OPTION,"Fast-SPK1/RANKI") == 0) {
2681 sys->J.fm = ranki_kw2;
2682 } else if (strcmp(FACTOR_OPTION,"Fast-SPK1/RANKI+ROW") == 0) {
2683 sys->J.fm = ranki_jz2;
2684 } else if (strcmp(FACTOR_OPTION,"Fastest-SPK1/MR-RANKI") == 0) {
2685 sys->J.fm = ranki_ba2;
2686 /* } else if (strcmp(FACTOR_OPTION,"GAUSS_EASY") == 0) {
2687 sys->J.fm = gauss_easy;
2688 } else if (strcmp(FACTOR_OPTION,"NGSLV-2-LEVEL") == 0) {
2689 sys->J.fm = ranki_kt2;*/
2690 } else {
2691 sys->J.fm = ranki_ba2;
2692 }
2693 linsolqr_set_pivot_zero(sys->J.sys, SING_TOL);
2694 linsolqr_set_drop_tolerance(sys->J.sys, sys->p.tolerance.drop);
2695 linsolqr_set_pivot_tolerance(sys->J.sys, PIVOT_TOL);
2696 /* this next one is fishy, but we don't use qr so not panicking */
2697 linsolqr_set_condition_tolerance(sys->J.sys, PIVOT_TOL);
2698 }
2699
2700 static
2701 void slv5_presolve(slv_system_t server, SlvClientToken asys)
2702 {
2703 struct var_variable **vp;
2704 struct rel_relation **rp;
2705 int32 cap, ind;
2706 int32 matrix_creation_needed = 1;
2707 slv5_system_t sys;
2708
2709 sys = SLV5(asys);
2710 iteration_begins(sys);
2711 check_system(sys);
2712 if( sys->vlist == NULL ) {
2713 FPRINTF(ASCERR,"ERROR: (slv5) slv5_presolve\n");
2714 FPRINTF(ASCERR," Variable list was never set.\n");
2715 return;
2716 }
2717 if( sys->rlist == NULL && sys->obj == NULL ) {
2718 FPRINTF(ASCERR,"ERROR: (slv5) slv5_presolve\n");
2719 FPRINTF(ASCERR," Relation list and objective never set.\n");
2720 return;
2721 }
2722
2723 if(sys->presolved > 0) { /* system has been presolved before */
2724 if(!slv5_dof_changed(sys) /* no changes in fixed or included flags */
2725 && PARTITION == sys->J.old_partition) {
2726 #if DEBUG
2727 FPRINTF(ASCERR,"YOU JUST AVOIDED MATRIX DESTRUCTION/CREATION\n");
2728 #endif
2729 matrix_creation_needed = 0;
2730 }
2731 }
2732
2733 rp=sys->rlist;
2734 for( ind = 0; ind < sys->rtot; ++ind ) {
2735 rel_set_satisfied(rp[ind],FALSE);
2736 }
2737 if( matrix_creation_needed ) {
2738
2739 cap = slv_get_num_solvers_rels(SERVER);
2740 sys->cap = slv_get_num_solvers_vars(SERVER);
2741 sys->cap = MAX(sys->cap,cap);
2742 vp=sys->vlist;
2743 for( ind = 0; ind < sys->vtot; ++ind ) {
2744 var_set_in_block(vp[ind],FALSE);
2745 }
2746 rp=sys->rlist;
2747 for( ind = 0; ind < sys->rtot; ++ind ) {
2748 rel_set_in_block(rp[ind],FALSE);
2749 rel_set_satisfied(rp[ind],FALSE);
2750 }
2751
2752 sys->presolved = 1; /* full presolve recognized here */
2753 sys->J.old_partition = PARTITION;
2754 destroy_matrices(sys);
2755 destroy_vectors(sys);
2756 create_matrices(server,sys);
2757 create_vectors(sys);
2758
2759 sys->s.block.current_reordered_block = -2;
2760 } else {
2761 slv5_update_linsolqr(sys);
2762 }
2763
2764 /* Reset status */
2765 sys->s.iteration = 0;
2766 sys->s.cpu_elapsed = 0.0;
2767 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
2768 sys->s.block.previous_total_size = 0;
2769 sys->s.costsize = 1+sys->s.block.number_of;
2770
2771 if( matrix_creation_needed ) {
2772 destroy_array(sys->s.cost);
2773 sys->s.cost = create_zero_array(sys->s.costsize,struct slv_block_cost);
2774 for( ind = 0; ind < sys->s.costsize; ++ind ) {
2775 sys->s.cost[ind].reorder_method = -1;
2776 }
2777 } else {
2778 reset_cost(sys->s.cost,sys->s.costsize);
2779 }
2780
2781 /* set to go to first unconverged block */
2782 sys->s.block.current_block = -1;
2783 sys->s.block.current_size = 0;
2784 sys->s.calc_ok = TRUE;
2785 sys->s.block.iteration = 0;
2786
2787 update_status(sys);
2788 iteration_ends(sys);
2789 sys->s.cost[sys->s.block.number_of].time=sys->s.cpu_elapsed;
2790 }
2791
2792
2793 static void slv5_resolve(slv_system_t server, SlvClientToken asys)
2794 {
2795 struct var_variable **vp;
2796 struct rel_relation **rp;
2797 slv5_system_t sys;
2798 (void) server;
2799 sys = SLV5(asys);
2800
2801 check_system(sys);
2802 for( vp = sys->vlist ; *vp != NULL ; ++vp ) {
2803 var_set_in_block(*vp,FALSE);
2804 }
2805 for( rp = sys->rlist ; *rp != NULL ; ++rp ) {
2806 rel_set_in_block(*rp,FALSE);
2807 rel_set_satisfied(*rp,FALSE);
2808 }
2809
2810 /* Reset status */
2811 sys->s.iteration = 0;
2812 sys->s.cpu_elapsed = 0.0;
2813 sys->s.converged = sys->s.diverged = sys->s.inconsistent = FALSE;
2814 sys->s.block.previous_total_size = 0;
2815
2816 /* go to first unconverged block */
2817 sys->s.block.current_block = -1;
2818 sys->s.block.current_size = 0;
2819 sys->s.calc_ok = TRUE;
2820 sys->s.block.iteration = 0;
2821
2822 update_status(sys);
2823 }
2824
2825
2826 static void slv5_iterate(slv_system_t server, SlvClientToken asys)
2827 {
2828 slv5_system_t sys;
2829 FILE *mif;
2830 FILE *lif;
2831 boolean first, new_ok, descent_ok;
2832 int32 m, ds_status, rank_defect;
2833 int32 optimizing;
2834 real64 oldpsi;
2835 real64 time0;
2836 real64 factor;
2837
2838 sys = SLV5(asys);
2839 mif = MIF(sys);
2840 lif = LIF(sys);
2841 ds_status = 0;
2842 rank_defect = 0;
2843
2844 if (server == NULL || sys==NULL) return;
2845 if (check_system(SLV5(sys))) return;
2846 if( !sys->s.ready_to_solve ) {
2847 FPRINTF(ASCERR,"ERROR: (slv5) slv5_iterate\n");
2848 FPRINTF(ASCERR," Not ready to solve.\n");
2849 return;
2850 }
2851
2852 if (sys->s.block.current_block==-1) {
2853 find_next_unconverged_block(sys);
2854 update_status(sys);
2855 return;
2856 }
2857 if (SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
2858 LIFDS)) {
2859 debug_delimiter(lif);
2860 }
2861 iteration_begins(sys);
2862
2863 optimizing = sys->obj ? (sys->vused - sys->rank) : 0;
2864
2865 if( optimizing ) {
2866 FPRINTF(ASCERR,"ERROR: (slv5) slv5_iterate\n");
2867 FPRINTF(ASCERR," IPSlv cannot presently optimize.\n");
2868 sys->s.diverged = 1;
2869 iteration_ends(sys);
2870 update_status(sys);
2871 return;
2872 }
2873
2874 /*
2875 * Attempt direct solve if appropriate
2876 */
2877 if( sys->s.block.iteration == 1 && sys->s.block.current_size == 1 ) {
2878 struct var_variable *var;
2879 struct rel_relation *rel;
2880 var = sys->vlist[mtx_col_to_org(sys->J.mtx,sys->J.reg.col.low)];
2881 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,sys->J.reg.row.low)];
2882 if (SHOW_LESS_IMPT && LIFDS) {
2883 FPRINTF(lif,"%-40s ---> (%d)", "Singleton relation",
2884 mtx_row_to_org(sys->J.mtx,sys->J.reg.row.low));
2885 print_rel_name(lif,sys,rel); PUTC('\n',lif);
2886 FPRINTF(lif,"%-40s ---> (%d)", "Singleton variable",
2887 mtx_col_to_org(sys->J.mtx,sys->J.reg.col.low));
2888 print_var_name(lif,sys,var); PUTC('\n',lif);
2889 }
2890
2891 /* Attempt direct solve */
2892 time0=tm_cpu_time();
2893 ds_status=slv_direct_solve(SERVER,rel,var,mif,FEAS_TOL,IGNORE_BOUNDS,0);
2894 sys->s.block.functime += (tm_cpu_time()-time0);
2895
2896 switch( ds_status ) {
2897 case 0:
2898 if (SHOW_LESS_IMPT) {
2899 FPRINTF(lif,"Unable to directly solve symbolically.\n");
2900 }
2901 break;
2902 case 1:
2903 if (SHOW_LESS_IMPT && LIFDS) {
2904 FPRINTF(lif,"Directly solved.\n");
2905 }
2906 sys->s.calc_ok = calc_residuals(sys);
2907 if (sys->s.calc_ok) {
2908 iteration_ends(sys);
2909 find_next_unconverged_block(sys);
2910 update_status(sys);
2911 return;
2912 }
2913 FPRINTF(mif,"Direct solve found numerically impossible equation\n");
2914 FPRINTF(mif,"given variables solved in previous blocks.\n");
2915 case -1:
2916 sys->s.inconsistent = TRUE;
2917 FPRINTF(mif,"No solution exists within the bounds given for:\n");
2918 print_var_name(mif,sys,var); PUTC('\n',mif);
2919 FPRINTF(mif,"when inverting relation:\n");
2920 print_rel_name(mif,sys,rel); PUTC('\n',mif);
2921 iteration_ends(sys);
2922 update_status(sys);
2923 return;
2924 }
2925 } /* if direct solve fails with a 0, go on to newton a 1x1 */
2926 if( !calc_J(sys) ) {
2927 FPRINTF(MIF(sys),"Jacobian calculation errors detected.\n");
2928 }
2929
2930 calc_mu(sys); /* Complementary gap */
2931 calc_psi(sys);
2932
2933 scale_system(sys);
2934
2935 set_factor_options(sys);
2936 rank_defect = calc_pivots(sys);
2937
2938 calc_newton(sys);
2939 calc_varnewstep(sys); /* only call */
2940
2941 factor = factor_for_complementary_vars(sys,0);
2942 factor = TOWARD_BOUNDS * factor;
2943
2944 apply_step(sys,0,factor);
2945 calc_newton_residuals(sys);
2946 calc_muaff(sys); /* Complementary gap */
2947 calc_sigma(sys); /* Penalty Parameter */
2948 calc_sigmamu(sys); /* Penalty Parameter time complementary gap */
2949 restore_variables(sys);
2950 apply_2nd_order_correction(sys);
2951 calc_correction(sys);
2952 restore_variables(sys);
2953 calc_perturbed_residuals(sys);
2954 scale_perturbed_system(sys);
2955 calc_perturbed_newton(sys);
2956
2957 calc_varstep(sys);
2958
2959 factor = quadratic_factor_for_complementary_vars(sys);
2960 if (factor == 0.0) {
2961 FPRINTF(ASCERR,"ERROR: (slv5) slv5_iterate\n");
2962 FPRINTF(ASCERR,"Error in bisection method\n");
2963 sys->s.diverged = 1;
2964 iteration_ends(sys);
2965 update_status(sys);
2966 return;
2967 }
2968 factor = TOWARD_BOUNDS * factor;
2969
2970 first = TRUE;
2971 new_ok= FALSE;
2972 descent_ok= FALSE;
2973 oldpsi = sys->psi;
2974 m = 0;
2975 #if DEBUG_ITERATION
2976 FPRINTF(ASCERR,"m = %d \n",m);
2977 #endif /* DEBUG_ITERATION */
2978 while (first || !new_ok || !descent_ok){
2979 if (first) {
2980 first = FALSE;
2981 } else {
2982 if (!new_ok || !descent_ok ) {
2983 sys->psi = oldpsi;
2984 restore_variables(sys);
2985 m = m + 1;
2986 factor = factor * pow(RHO,m);
2987 #if DEBUG_ITERATION
2988 FPRINTF(ASCERR,"m = %d \n",m);
2989 FPRINTF(ASCERR,"factor = %g \n",factor);
2990 #endif /* DEBUG_ITERATION */
2991 if (m >= 10) {
2992 FPRINTF(ASCERR,"ERROR: (slv5) slv5_iterate\n");
2993 FPRINTF(ASCERR,"Coefficient m is too big \n");
2994 sys->s.diverged = 1;
2995 iteration_ends(sys);
2996 update_status(sys);
2997 return;
2998 }
2999 }
3000 }
3001 if (sys->eta > 0.0) {
3002 sys->progress = -1.0 * ALPHA * factor * (1.0 - sys->sigma) *
3003 (sys->eta - sys->comp);
3004 } else {
3005 sys->progress = 0.0 ;
3006 }
3007
3008 #if DEBUG_ITERATION
3009 FPRINTF(ASCERR,"Minimum expected progress = %g \n",sys->progress);
3010 #endif /* DEBUG_ITERATION */
3011
3012 /*
3013 * Apply step.
3014 */
3015 apply_quadratic_step(sys,factor);
3016
3017 /*
3018 * Check calculations at new point.
3019 */
3020 new_ok = (calc_residuals(sys));
3021
3022 if( !new_ok ) {
3023 continue;
3024 }
3025
3026 /*
3027 * Check for descent.
3028 */
3029 calc_psi(sys);
3030