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

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