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

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

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