/[ascend]/trunk/solvers/ipslv/slv5.mehrotra_no_quad.c
ViewVC logotype

Contents of /trunk/solvers/ipslv/slv5.mehrotra_no_quad.c

Parent Directory Parent Directory | Revision Log Revision Log


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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22