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

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