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

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

Parent Directory Parent Directory | Revision Log Revision Log


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