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