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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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