/[ascend]/trunk/ascend4/solver/slv0.c
ViewVC logotype

Contents of /trunk/ascend4/solver/slv0.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2 - (show annotations) (download) (as text)
Fri Nov 5 21:21:36 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 106392 byte(s)
corrected runaway minor loop problem exposed by Helder
1 /*
2 * SLV: Ascend Nonlinear Solver
3 * by Karl Michael Westerberg
4 * Created: 2/6/90
5 * Version: $Revision: 1.33 $
6 * Version control file: $RCSfile: slv0.c,v $
7 * Date last modified: $Date: 2000/01/25 02:27:16 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1990 Karl Michael Westerberg
13 * Copyright (C) 1993 Joseph Zaher
14 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15 *
16 * The SLV solver is free software; you can redistribute
17 * it and/or modify it under the terms of the GNU General Public License as
18 * published by the Free Software Foundation; either version 2 of the
19 * License, or (at your option) any later version.
20 *
21 * The SLV solver is distributed in hope that it will be
22 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24 * General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with the program; if not, write to the Free Software Foundation,
28 * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
29 * COPYING. COPYING is found in ../compiler.
30 *
31 */
32
33 #include <stdarg.h>
34 #include "utilities/ascConfig.h"
35 #include "utilities/ascMalloc.h"
36 #include "utilities/ascPanic.h"
37 #include "general/list.h"
38 #include "utilities/set.h"
39 #include "general/tm_time.h"
40 #include "utilities/mem.h"
41 #include "solver/mtx.h"
42 #include "solver/linsol.h"
43 #include "solver/linsolqr.h"
44 #include "solver/slv_types.h"
45 #include "solver/var.h"
46 #include "solver/rel.h"
47 #include "solver/discrete.h"
48 #include "solver/conditional.h"
49 #include "solver/logrel.h"
50 #include "solver/bnd.h"
51 #include "solver/calc.h"
52 #include "solver/relman.h"
53 #include "solver/slv_common.h"
54 #include "solver/slv_client.h"
55 #include "solver/slv0.h"
56 #if !defined(STATIC_SLV) && !defined(DYNAMIC_SLV)
57 /* do nothing */
58 int slv0_register(SlvFunctionsT *f)
59 {
60 (void)f; /* stop gcc whine about unused parameter */
61
62 FPRINTF(stderr,"Slv not compiled in this ASCEND IV.\n");
63 return 1;
64 }
65 #else /* either STATIC_SLV or DYNAMIC_SLV is defined */
66 #ifdef DYNAMIC_SLV
67 /* do dynamic loading stuff. yeah, right */
68 #else /* following is used if STATIC_SLV is defined */
69
70
71 #define KILL 0
72 /* code that needs to be deleted compiles only with kill = 1 */
73 #define DEBUG TRUE
74
75 #define slv0_solver_name "Slv" /* Solver's name. don't mess with the caps!*/
76 #define slv0_solver_number 0 /* Solver's number */
77
78 #define SLV0(s) ((slv0_system_t)(s))
79
80 #define slv0_IA_SIZE 2
81 #define slv0_RA_SIZE 0
82 #define slv0_CA_SIZE 0
83 #define slv0_VA_SIZE 0
84
85 /* subscripts for ia */
86 #define SP0_LIFDS 0
87 #define SP0_SAVLIN 1
88
89 /*********************************************************************\
90 Subparameters implemented: (value/meaning)
91 sp.ia[SP0_LIFDS] 0=>do not show full detail info for singletons
92 1=>do (this value ignored if detailed solve info not on.
93 sp.ia[SP0_SAVLIN] 0=>do not append linearizations arising in the newton
94 scheme to the file SlvLinsol.dat.
95 1=>do.
96 \*********************************************************************/
97 struct update_data {
98 int jacobian; /* Countdown on jacobian updating */
99 int weights; /* Countdown on weights updating */
100 int nominals; /* Countdown on nominals updating */
101 };
102
103 struct jacobian_data {
104 linsol_system_t sys; /* Linear system */
105 mtx_matrix_t mtx; /* Transpose gradient of residuals */
106 real64 *rhs; /* RHS from linear system */
107 unsigned *varpivots; /* Pivoted variables */
108 unsigned *relpivots; /* Pivoted relations */
109 unsigned *subregions; /* Set of subregion indeces */
110 mtx_region_t reg; /* Current block region */
111 int32 rank; /* Numerical rank of the jacobian */
112 boolean accurate; /* ? Recalculate matrix */
113 boolean singular; /* ? Can matrix be inverted */
114 };
115
116 struct hessian_data {
117 struct vector_data Bs; /* Product of B and s */
118 struct vector_data y; /* Difference in stationaries */
119 real64 ys; /* inner product of y and s */
120 real64 sBs; /* inner product of s and Bs */
121 struct hessian_data *next; /* previous iteration data */
122 };
123
124 struct reduced_data {
125 real64 **mtx; /* Dense matrix */
126 real64 *ZBs; /* Reduced Bs */
127 real64 *Zy; /* Reduced y */
128 int32 order; /* Degrees of freedom */
129 boolean accurate; /* Ready to re-compute ? */
130 };
131
132 struct slv0_system_structure {
133
134 /**
135 *** Problem definition
136 **/
137 slv_system_t slv; /* slv_system_t back-link */
138 struct rel_relation *obj; /* Objective function: NULL = none */
139 struct var_variable **vlist; /* Variable list (NULL terminated) */
140 struct var_variable **vlist_user; /* User vlist (NULL=determine) */
141 struct rel_relation **rlist; /* Relation list (NULL terminated) */
142 struct rel_relation **rlist_user; /* User rlist (NULL = none) */
143 struct ExtRelCache **erlist; /* External relations cache list */
144 struct ExtRelCache **erlist_user; /* User erlist (NULL = none) */
145
146 /**
147 *** Solver information
148 **/
149 int integrity; /* ? Has the system been created */
150 slv_parameters_t p; /* Parameters */
151 slv_status_t s; /* Status (as of iteration end) */
152 struct update_data update; /* Jacobian frequency counters */
153 int32 cap; /* Order of matrix/vectors */
154 int32 rank; /* Symbolic rank of problem */
155 int32 vused; /* Free and incident variables */
156 int32 vtot; /* length of varlist */
157 int32 bused; /* Included boundaries */
158 int32 rused; /* Included relations */
159 int32 rtot; /* length of varlist */
160 double clock; /* CPU time */
161 int iarray[slv0_IA_SIZE];
162
163 /**
164 *** Calculated data (scaled)
165 **/
166 struct jacobian_data J; /* linearized system */
167 struct hessian_data *B; /* Curvature information */
168 struct reduced_data ZBZ; /* Reduced hessian */
169
170 struct vector_data nominals; /* Variable nominals */
171 struct vector_data weights; /* Relation weights */
172 struct vector_data variables; /* Variable values */
173 struct vector_data residuals; /* Relation residuals */
174 struct vector_data gradient; /* Objective gradient */
175 struct vector_data multipliers; /* Relation multipliers */
176 struct vector_data stationary; /* Lagrange gradient */
177 struct vector_data gamma; /* Feasibility steepest descent */
178 struct vector_data Jgamma; /* Product of J and gamma */
179 struct vector_data newton; /* Dependent variables */
180 struct vector_data Bnewton; /* Product of B and newton */
181 struct vector_data nullspace; /* Independent variables */
182 struct vector_data varstep1; /* 1st order in variables */
183 struct vector_data Bvarstep1; /* Product of B and varstep1 */
184 struct vector_data varstep2; /* 2nd order in variables */
185 struct vector_data Bvarstep2; /* Product of B and varstep2 */
186 struct vector_data mulstep1; /* 1st order in multipliers */
187 struct vector_data mulstep2; /* 2nd order in multipliers */
188 struct vector_data varstep; /* Step in variables */
189 struct vector_data mulstep; /* Step in multipliers */
190
191 real64 objective; /* Objective function evaluation */
192 real64 phi; /* Unconstrained minimizer */
193 real64 maxstep; /* Maximum step size allowed */
194 real64 progress; /* Steepest directional derivative */
195 };
196
197
198 /**
199 *** Integrity checks
200 *** ----------------
201 *** check_system(sys)
202 **/
203
204 #define OK ((int)813028392)
205 #define DESTROYED ((int)103289182)
206 static int check_system(sys)
207 slv0_system_t sys;
208 /**
209 *** Checks sys for NULL and for integrity.
210 **/
211 {
212 if( sys == NULL ) {
213 FPRINTF(stderr,"ERROR: (slv0) check_system\n");
214 FPRINTF(stderr," NULL system handle.\n");
215 return 1;
216 }
217
218 switch( sys->integrity ) {
219 case OK:
220 return 0;
221 case DESTROYED:
222 FPRINTF(stderr,"ERROR: (slv0) check_system\n");
223 FPRINTF(stderr," System was recently destroyed.\n");
224 return 1;
225 default:
226 FPRINTF(stderr,"ERROR: (slv0) check_system\n");
227 FPRINTF(stderr," System reused or never allocated.\n");
228 return 1;
229 }
230 }
231
232
233 /**
234 *** General input/output routines
235 *** -----------------------------
236 *** print_var_name(out,sys,var)
237 *** print_bnd_name(out,sys,bnd)
238 *** print_rel_name(out,sys,rel)
239 **/
240
241 #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
242 #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
243 #define print_bnd_name(a,b,c) slv_print_bnd_name((a),(c))
244
245 /**
246 *** Debug output routines
247 *** ---------------------
248 *** debug_delimiter(fp)
249 *** debug_out_vector(fp,vec)
250 *** debug_out_var_values(fp,sys)
251 *** debug_out_rel_residuals(fp,sys)
252 *** debug_out_jacobian(fp,sys)
253 *** debug_out_hessian(fp,sys)
254 *** debug_write_array(fp,real64 *,length)
255 **/
256
257 static void debug_delimiter(fp)
258 FILE *fp;
259 /**
260 *** Outputs a hyphenated line.
261 **/
262 {
263 int i;
264 for( i=0; i<60; i++ ) PUTC('-',fp);
265 PUTC('\n',fp);
266 }
267
268
269 static void debug_out_vector(fp,sys,vec)
270 FILE *fp;
271 slv0_system_t sys;
272 struct vector_data *vec;
273 /**
274 *** Outputs a vector.
275 **/
276 {
277 int32 ndx;
278 FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n",
279 calc_sqrt_D0(vec->norm2), vec->accurate?"TRUE":"FALSE",
280 vec->rng->low,vec->rng->high);
281 FPRINTF(fp,"Vector --> ");
282 for( ndx=vec->rng->low ; ndx<=vec->rng->high ; ++ndx )
283 FPRINTF(fp, "%g ", vec->vec[ndx]);
284 PUTC('\n',fp);
285 }
286
287 static void debug_out_var_values(fp,sys)
288 FILE *fp;
289 slv0_system_t sys;
290 /**
291 *** Outputs all variable values in current block.
292 **/
293 {
294 int32 col;
295 int32 orig_col;
296
297 FPRINTF(fp,"Var values --> ");
298 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high ; col++ ) {
299 struct var_variable *var;
300 orig_col = mtx_col_to_org(sys->J.mtx,col);
301 var = sys->vlist[orig_col];
302 /* print_var_name(fp,sys,var); **** KAA_DEBUG */
303 FPRINTF(fp,"x[%d] = %12.8f\n", orig_col, var_value(var));
304 }
305 PUTC('\n',fp);
306 }
307
308 static void debug_out_rel_residuals(fp,sys)
309 FILE *fp;
310 slv0_system_t sys;
311 /**
312 *** Outputs all relation residuals in current block.
313 **/
314 {
315 int32 row;
316
317 FPRINTF(fp,"Rel residuals --> ");
318 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high ; row++ ) {
319 struct rel_relation *rel;
320 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
321 print_rel_name(fp,sys,rel);
322 FPRINTF(fp,"=%g ",rel_residual(rel));
323 }
324 PUTC('\n',fp);
325 }
326
327 static void debug_out_jacobian(fp,sys)
328 FILE *fp;
329 slv0_system_t sys;
330 /**
331 *** Outputs permutation and values of the nonzero elements in the
332 *** the jacobian matrix.
333 **/
334 {
335 mtx_coord_t nz;
336 real64 value;
337
338 nz.row = sys->J.reg.row.low;
339 for( ; nz.row <= sys->J.reg.row.high; ++(nz.row) ) {
340 FPRINTF(fp," Row %d (rel %d)\n",
341 nz.row,
342 mtx_row_to_org(sys->J.mtx,nz.row));
343 nz.col = mtx_FIRST;
344 while( value = mtx_next_in_row(sys->J.mtx,&nz,&(sys->J.reg.col)),
345 nz.col != mtx_LAST )
346 FPRINTF(fp," Col %d (var %d) has value %g\n",
347 nz.col, mtx_col_to_org(sys->J.mtx,nz.col), value);
348 }
349 }
350
351 static void debug_out_hessian(fp,sys)
352 FILE *fp;
353 slv0_system_t sys;
354 /**
355 *** Outputs permutation and values of the nonzero elements in the
356 *** reduced hessian matrix.
357 **/
358 {
359 mtx_coord_t nz;
360
361 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
362 nz.col = nz.row + sys->J.reg.col.high + 1 - sys->ZBZ.order;
363 FPRINTF(fp," ZBZ[%d (var %d)] = ",
364 nz.row, mtx_col_to_org(sys->J.mtx,nz.col));
365 for( nz.col = 0; nz.col < sys->ZBZ.order; nz.col++ )
366 FPRINTF(fp,"%10g ",sys->ZBZ.mtx[nz.row][nz.col]);
367 PUTC('\n',fp);
368 }
369 }
370
371
372
373 static void debug_write_array(FILE *fp,real64 *vec, int32 length)
374 {
375 int32 i;
376 for (i=0; i< length;i++)
377 FPRINTF(fp,"%.20g\n",vec[i]);
378 }
379
380
381 static char savlinfilename[]="SlvLinsol.dat. \0";
382 static char savlinfilebase[]="SlvLinsol.dat.\0";
383 static int savlinnum=0;
384 /** The number to postfix to savlinfilebase. increases with file accesses. **/
385
386 /**
387 *** Array/vector operations
388 *** ----------------------------
389 *** destroy_array(p)
390 *** create_array(len,type)
391 *** zero_array(arr,len,type)
392 ***
393 *** zero_vector(vec)
394 *** copy_vector(vec1,vec2)
395 *** prod = inner_product(vec1,vec2)
396 *** norm2 = square_norm(vec)
397 *** matrix_product(mtx,vec,prod,scale,transpose)
398 **/
399
400 #define destroy_array(p) \
401 if( (p) != NULL ) ascfree((p))
402 #define create_array(len,type) \
403 ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
404 #define create_zero_array(len,type) \
405 ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
406 #define zero_array(arr,nelts,type) \
407 mem_zero_byte_cast((arr),0,(nelts)*sizeof(type))
408 /* Zeros an array of nelts objects, each having given type. */
409
410 #define zero_vector(v) slv_zero_vector(v)
411 #define copy_vector(v,t) slv_copy_vector((v),(t))
412 #define inner_product(v,u) slv_inner_product((v),(u))
413 #define square_norm(v) slv_square_norm(v)
414 #define matrix_product(m,v,p,s,t) slv_matrix_product((m),(v),(p),(s),(t))
415
416 /**
417 *** Calculation routines
418 *** --------------------
419 *** ok = calc_objective(sys)
420 *** ok = calc_residuals(sys)
421 *** ok = calc_J(sys)
422 *** calc_nominals(sys)
423 *** calc_weights(sys)
424 *** scale_J(sys)
425 *** scale_variables(sys)
426 *** scale_residuals(sys)
427 *** ok = calc_gradient(sys)
428 *** calc_B(sys)
429 *** calc_ZBZ(sys)
430 *** calc_pivots(sys)
431 *** calc_rhs(sys)
432 *** calc_multipliers(sys)
433 *** calc_stationary(sys)
434 *** calc_newton(sys)
435 *** calc_Bnewton(sys)
436 *** calc_nullspace(sys)
437 *** calc_gamma(sys)
438 *** calc_Jgamma(sys)
439 *** calc_varstep1(sys)
440 *** calc_Bvarstep1(sys)
441 *** calc_varstep2(sys)
442 *** calc_Bvarstep2(sys)
443 *** calc_mulstep1(sys)
444 *** calc_mulstep2(sys)
445 *** calc_varstep(sys)
446 *** calc_mulstep(sys)
447 *** calc_phi(sys)
448 **/
449
450 #define UPDATE_JACOBIAN 1 /* Update jacobian this often */
451 #define UPDATE_WEIGHTS 1 /* Update weights this often */
452 #define UPDATE_NOMINALS 1000 /* Update nominals this often */
453 #define TOO_SMALL 1e-8
454 #define TOWARD_BOUNDS 0.95
455
456
457 #define OPTIMIZING(sys) ((sys)->ZBZ.order > 0)
458
459 static boolean calc_objective(sys)
460 slv0_system_t sys;
461 /**
462 *** Evaluate the objective function.
463 **/
464 {
465 calc_ok = TRUE;
466 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
467 sys->objective = (sys->obj ? relman_eval(sys->obj,&calc_ok) : 0.0);
468 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
469 return calc_ok;
470 }
471
472 static boolean calc_inequalities(sys)
473 slv0_system_t sys;
474 /**
475 *** Calculates all of the residuals of included inequalities.
476 *** Returns true iff (calculations preceded without error and
477 *** all inequalities are satisfied.)
478 **/
479 {
480 struct rel_relation **rp;
481 boolean satisfied=TRUE;
482 static rel_filter_t rfilter;
483 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
484 rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
485
486 calc_ok = TRUE;
487 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
488 for (rp=sys->rlist;*rp != NULL; rp++) {
489 if (rel_apply_filter(*rp,&rfilter)) {
490 relman_eval(*rp,&calc_ok);
491 satisfied=
492 satisfied && relman_calc_satisfied(*rp,sys->p.tolerance.feasible);
493 }
494 }
495 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
496 return (calc_ok && satisfied);
497 }
498
499 static boolean calc_residuals(sys)
500 slv0_system_t sys;
501 /**
502 *** Calculates all of the residuals in the current block and computes
503 *** the residual norm for block status. Returns true iff calculations
504 *** preceded without error.
505 **/
506 {
507 int32 row;
508 struct rel_relation *rel;
509 double time0;
510
511 if( sys->residuals.accurate )
512 return TRUE;
513
514 calc_ok = TRUE;
515 row = sys->residuals.rng->low;
516 time0=tm_cpu_time();
517 Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
518 for( ; row <= sys->residuals.rng->high; row++ ) {
519 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
520 #if DEBUG
521 if (!rel) {
522 int r;
523 r=mtx_row_to_org(sys->J.mtx,row);
524 FPRINTF(stderr,"NULL relation found !!\n");
525 FPRINTF(stderr,"at row %d rel %d in calc_residuals\n",(int)row,r);
526 FFLUSH(stderr);
527 }
528 #endif
529 sys->residuals.vec[row] = relman_eval(rel,&calc_ok);
530 /*
531 FPRINTF(stderr,"RESIDUAL[%d] = %12.8g\n",
532 mtx_row_to_org(sys->J.mtx,row), sys->residuals.vec[row]);
533 FFLUSH(stderr);
534 */
535
536 relman_calc_satisfied(rel,sys->p.tolerance.feasible);
537 }
538 Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
539 sys->s.block.functime += (tm_cpu_time() -time0);
540 sys->s.block.funcs++;
541 square_norm( &(sys->residuals) );
542 sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
543 return(calc_ok);
544 }
545
546
547 static boolean calc_J(sys)
548 slv0_system_t sys;
549 /**
550 *** Calculates the current block of the jacobian.
551 *** It is initially unscaled.
552 **/
553 {
554 int32 row;
555 var_filter_t vfilter;
556 double time0;
557
558 if( sys->J.accurate )
559 return TRUE;
560
561 calc_ok = TRUE;
562 vfilter.matchbits = (VAR_INBLOCK | VAR_ACTIVE);
563 vfilter.matchvalue = (VAR_INBLOCK | VAR_ACTIVE);
564 /* vfilter.fixed = var_ignore;
565 vfilter.incident = var_ignore;
566 vfilter.in_block = var_true; */
567 time0=tm_cpu_time();
568 mtx_clear_region(sys->J.mtx,&(sys->J.reg));
569 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
570 struct rel_relation *rel;
571 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
572 relman_diffs(rel,&vfilter,sys->J.mtx);
573 }
574 sys->s.block.jactime += (tm_cpu_time() - time0);
575 sys->s.block.jacs++;
576
577 if( --(sys->update.nominals) <= 0 )
578 sys->nominals.accurate = FALSE;
579 if( --(sys->update.weights) <= 0 )
580 sys->weights.accurate = FALSE;
581
582 linsol_matrix_was_changed(sys->J.sys);
583 return(calc_ok);
584 }
585
586
587 static void calc_nominals(sys)
588 slv0_system_t sys;
589 /**
590 *** Retrieves the nominal values of all of the block variables,
591 *** insuring that they are all strictly positive.
592 **/
593 {
594 int32 col;
595
596 if( sys->nominals.accurate )
597 return;
598
599 col = sys->nominals.rng->low;
600 for( ; col <= sys->nominals.rng->high; col++ ) {
601 struct var_variable *var;
602 real64 n;
603 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
604 n = var_nominal(var);
605 if( n <= 0.0 ) {
606 FILE *fp = MIF(sys);
607 if( n == 0.0 ) {
608 n = TOO_SMALL;
609 FPRINTF(fp,"ERROR: (slv0) calc_nominals\n");
610 FPRINTF(fp," Variable ");
611 print_var_name(fp,sys,var);
612 FPRINTF(fp," \nhas nominal value of zero.\n");
613 FPRINTF(fp," Resetting to %g.\n",n);
614 var_set_nominal(var,n);
615 } else {
616 n = -n;
617 FPRINTF(fp,"ERROR: (slv0) calc_nominals\n");
618 FPRINTF(fp," Variable ");
619 print_var_name(fp,sys,var);
620 FPRINTF(fp," \nhas negative nominal value.\n");
621 FPRINTF(fp," Resetting to %g.\n",n);
622 var_set_nominal(var,n);
623 }
624 }
625 sys->nominals.vec[col] = n;
626 }
627 square_norm( &(sys->nominals) );
628 sys->update.nominals = UPDATE_NOMINALS;
629 sys->nominals.accurate = TRUE;
630 }
631
632
633 static void calc_weights(sys)
634 slv0_system_t sys;
635 /**
636 *** Calculates the weights of all of the block relations
637 *** to scale the rows of the Jacobian.
638 **/
639 {
640 mtx_coord_t nz;
641
642 if( sys->weights.accurate )
643 return;
644
645 nz.row = sys->weights.rng->low;
646 for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
647 real64 sum;
648 sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
649 sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
650 }
651 square_norm( &(sys->weights) );
652 sys->update.weights = UPDATE_WEIGHTS;
653 sys->residuals.accurate = FALSE;
654 sys->weights.accurate = TRUE;
655 }
656
657
658 static void scale_J(sys)
659 slv0_system_t sys;
660 /**
661 *** Scales the jacobian.
662 **/
663 {
664 int32 row;
665 int32 col;
666
667 if( sys->J.accurate )
668 return;
669
670 calc_nominals(sys);
671
672 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ )
673 mtx_mult_col(sys->J.mtx,col,sys->nominals.vec[col],&(sys->J.reg.row));
674
675 calc_weights(sys);
676 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ )
677 mtx_mult_row(sys->J.mtx,row,sys->weights.vec[row],&(sys->J.reg.col));
678
679 sys->update.jacobian = UPDATE_JACOBIAN;
680 sys->J.accurate = TRUE;
681 sys->J.singular = FALSE; /* yet to be determined */
682 #if DEBUG
683 FPRINTF(LIF(sys),"\nJacobian: \n");
684 debug_out_jacobian(LIF(sys),sys);
685 #endif
686 }
687
688
689 static void scale_variables(sys)
690 slv0_system_t sys;
691 {
692 int32 col;
693
694 if( sys->variables.accurate )
695 return;
696
697 col = sys->variables.rng->low;
698 for( ; col <= sys->variables.rng->high; col++ ) {
699 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
700 sys->variables.vec[col] = var_value(var)/sys->nominals.vec[col];
701 }
702 square_norm( &(sys->variables) );
703 sys->variables.accurate = TRUE;
704 #if DEBUG
705 FPRINTF(LIF(sys),"Variables: ");
706 debug_out_vector(LIF(sys),sys,&(sys->variables));
707 #endif
708 }
709
710
711 static void scale_residuals(sys)
712 slv0_system_t sys;
713 /**
714 *** Scales the previously calculated residuals.
715 **/
716 {
717 int32 row;
718
719 if( sys->residuals.accurate )
720 return;
721
722 row = sys->residuals.rng->low;
723 for( ; row <= sys->residuals.rng->high; row++ ) {
724 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
725 sys->residuals.vec[row] = rel_residual(rel)*sys->weights.vec[row];
726 }
727 square_norm( &(sys->residuals) );
728 sys->residuals.accurate = TRUE;
729 #if DEBUG
730 FPRINTF(LIF(sys),"Residuals: ");
731 debug_out_vector(LIF(sys),sys,&(sys->residuals));
732 #endif
733 }
734
735 static boolean calc_gradient(sys)
736 slv0_system_t sys;
737 /**
738 *** Calculate scaled gradient of the objective function.
739 **/
740 {
741 if( sys->gradient.accurate )
742 return TRUE;
743
744 calc_ok = TRUE;
745 if ( !OPTIMIZING(sys) ) {
746 zero_vector(&(sys->gradient));
747 sys->gradient.norm2 = 0.0;
748 } else {
749 struct var_variable **vp,**tmp;
750 var_filter_t vfilter;
751 vfilter.matchbits = (VAR_INBLOCK | VAR_ACTIVE);
752 vfilter.matchvalue = (VAR_INBLOCK | VAR_ACTIVE);
753 /* vfilter.fixed = var_ignore;
754 vfilter.incident = var_ignore;
755 vfilter.in_block = var_true; */
756 zero_vector(&(sys->gradient));
757 #if 0 /* this needs to be reimplemented */
758 for( vp=tmp=expr_incidence_list(sys->obj,NULL) ; *vp != NULL ; ++vp ) {
759 int32 col;
760 col = mtx_org_to_col(sys->J.mtx,var_sindex(*vp));
761 if( var_apply_filter(*vp,&vfilter) )
762 sys->gradient.vec[col] = sys->nominals.vec[col]*
763 exprman_diff(sys->obj,*vp);
764 }
765 ascfree(tmp);
766 square_norm( &(sys->gradient) );
767 #else
768 Asc_Panic(2, "calc_gradient",
769 "need to reimplemnt objective eval slv0\n");
770 #endif
771 }
772 sys->gradient.accurate = TRUE;
773 #if DEBUG
774 FPRINTF(LIF(sys),"Gradient: ");
775 debug_out_vector(LIF(sys),sys,&(sys->gradient));
776 #endif
777 return calc_ok;
778 }
779
780
781 static void create_update(sys)
782 slv0_system_t sys;
783 /**
784 *** Create a new hessian_data structure for storing
785 *** latest update information.
786 **/
787 {
788 struct hessian_data *update;
789
790 if( !OPTIMIZING(sys) )
791 return;
792
793 update = (struct hessian_data *)ascmalloc(sizeof(struct hessian_data));
794 update->y.vec = create_array(sys->cap,real64);
795 update->y.rng = &(sys->J.reg.col);
796 update->y.accurate = FALSE;
797 update->Bs.vec = create_array(sys->cap,real64);
798 update->Bs.rng = &(sys->J.reg.col);
799 update->Bs.accurate = FALSE;
800 update->next = sys->B;
801 sys->B = update;
802 }
803
804
805 #define POSITIVE_DEFINITE 1.0e-2
806 static void calc_B(sys)
807 slv0_system_t sys;
808 /**
809 *** Computes a rank 2 BFGS update to the hessian matrix
810 *** B which accumulates curvature.
811 **/
812 {
813 if( sys->s.block.iteration > 1 )
814 create_update(sys);
815 else
816 if( sys->B ) {
817 struct hessian_data *update;
818 for( update=sys->B; update != NULL; ) {
819 struct hessian_data *handle;
820 handle = update;
821 update = update->next;
822 destroy_array(handle->y.vec);
823 destroy_array(handle->Bs.vec);
824 ascfree(handle);
825 }
826 sys->B = NULL;
827 }
828
829 if( sys->B ) {
830 real64 theta;
831 /**
832 *** The y vector
833 **/
834 if( !sys->B->y.accurate ) {
835 int32 col;
836 matrix_product(sys->J.mtx, &(sys->multipliers),
837 &(sys->B->y), 1.0, TRUE);
838 col = sys->B->y.rng->low;
839 for( ; col <= sys->B->y.rng->high; col++ )
840 sys->B->y.vec[col] += sys->gradient.vec[col] -
841 sys->stationary.vec[col];
842 square_norm( &(sys->B->y) );
843 sys->B->y.accurate = TRUE;
844 }
845
846 /**
847 *** The Bs vector
848 **/
849 if( !sys->B->Bs.accurate ) {
850 struct hessian_data *update;
851 copy_vector(&(sys->varstep),&(sys->B->Bs));
852 for( update=sys->B->next; update != NULL; update = update->next ) {
853 int32 col;
854 real64 ys = inner_product( &(update->y),&(sys->varstep) );
855 real64 sBs = inner_product( &(update->Bs),&(sys->varstep) );
856 col = sys->B->Bs.rng->low;
857 for( ; col<=sys->B->Bs.rng->high; col++) {
858 sys->B->Bs.vec[col] += update->ys > 0.0 ?
859 (update->y.vec[col])*ys/update->ys : 0.0;
860 sys->B->Bs.vec[col] -= update->sBs > 0.0 ?
861 (update->Bs.vec[col])*sBs/update->sBs : 0.0;
862 }
863 }
864 square_norm( &(sys->B->Bs) );
865 sys->B->Bs.accurate = TRUE;
866 }
867
868 sys->B->ys = inner_product( &(sys->B->y),&(sys->varstep) );
869 sys->B->sBs = inner_product( &(sys->B->Bs),&(sys->varstep) );
870
871 if( sys->B->ys == 0.0 && sys->B->sBs == 0.0 ) theta = 0.0;
872 else theta = sys->B->ys < POSITIVE_DEFINITE*sys->B->sBs ?
873 (1.0-POSITIVE_DEFINITE)*sys->B->sBs/(sys->B->sBs - sys->B->ys):1.0;
874 #if DEBUG
875 FPRINTF(LIF(sys),"ys, sBs, PD, theta = %g, %g, %g, %g\n",
876 sys->B->ys,
877 sys->B->sBs,
878 POSITIVE_DEFINITE,
879 theta);
880 #endif
881 if( theta < 1.0 ) {
882 int32 col;
883 col = sys->B->y.rng->low;
884 for( ; col <= sys->B->y.rng->high; col++ )
885 sys->B->y.vec[col] = theta*sys->B->y.vec[col] +
886 (1.0-theta)*sys->B->Bs.vec[col];
887 square_norm( &(sys->B->y) );
888 sys->B->ys = theta*sys->B->ys + (1.0-theta)*sys->B->sBs;
889 }
890 }
891 }
892 #undef POSITIVE_DEFINITE
893
894
895 static void calc_pivots(sys)
896 slv0_system_t sys;
897 /**
898 *** Obtain the equations and variables which
899 *** are able to be pivoted.
900 **/
901 {
902 linsol_system_t lsys = sys->J.sys;
903 FILE *fp = LIF(sys);
904 #undef KAA_DEBUG
905 #ifdef KAA_DEBUG
906 double time1 = tm_cpu_time();
907 linsol_invert(lsys,&(sys->J.reg));
908 time1 = tm_cpu_time() - time1;
909 FPRINTF(stdout,"Time to invert block = %g seconds\n",time1);
910 #else
911 linsol_invert(lsys,&(sys->J.reg));
912 #endif
913
914 set_null(sys->J.relpivots,sys->cap);
915 set_null(sys->J.varpivots,sys->cap);
916 linsol_get_pivot_sets(lsys,sys->J.relpivots,sys->J.varpivots);
917 sys->J.rank = linsol_rank(lsys);
918 sys->J.singular = FALSE;
919 if( sys->J.rank < sys->J.reg.row.high-sys->J.reg.row.low+1 ) {
920 int32 row;
921 sys->J.singular = TRUE;
922 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high ; row++ ) {
923 int32 org_row;
924 org_row = mtx_row_to_org(sys->J.mtx,row);
925 if( !set_is_member(sys->J.relpivots,org_row) ) {
926 struct rel_relation *rel;
927 rel = sys->rlist[org_row];
928 FPRINTF(fp,"%-40s ---> ","Relation not pivoted");
929 print_rel_name(fp,sys,rel);
930 PUTC('\n',fp);
931
932 /**
933 *** assign zeros to the corresponding weights
934 *** so that subsequent calls to "scale_residuals"
935 *** will only measure the pivoted equations.
936 **/
937 sys->weights.vec[row] = 0.0;
938 sys->residuals.vec[row] = 0.0;
939 sys->residuals.accurate = FALSE;
940 mtx_mult_row(sys->J.mtx,row,0.0,&(sys->J.reg.col));
941 }
942 }
943 if( !sys->residuals.accurate ) {
944 square_norm( &(sys->residuals) );
945 sys->residuals.accurate = TRUE;
946 sys->update.weights = 0; /* re-compute weights next iteration. */
947 }
948 }
949 if( sys->J.rank < sys->J.reg.col.high-sys->J.reg.col.low+1 ) {
950 int32 col;
951 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high ; col++ ) {
952 int32 org_col;
953 org_col = mtx_col_to_org(sys->J.mtx,col);
954 if( !set_is_member(sys->J.varpivots,org_col) ) {
955 struct var_variable *var;
956 var = sys->vlist[org_col];
957 FPRINTF(fp,"%-40s ---> ","Variable not pivoted");
958 print_var_name(fp,sys,var);
959 PUTC('\n',fp);
960 /**
961 *** If we're not optimizing (everything should be
962 *** pivotable) or this was one of the dependent variables,
963 *** consider this variable as if it were fixed.
964 **/
965 if( col <= sys->J.reg.col.high - sys->ZBZ.order )
966 mtx_mult_col(sys->J.mtx,col,0.0,&(sys->J.reg.row));
967 }
968 }
969 }
970 FPRINTF(fp,"%-40s ---> %d (%s)\n","Jacobian rank", sys->J.rank,
971 sys->J.singular ? "deficient":"full");
972 if (sys->p.output.less_important) {
973 FPRINTF(LIF(sys),"%-40s ---> %g\n","Smallest pivot",
974 linsol_smallest_pivot(sys->J.sys));
975 }
976 }
977
978
979 static void calc_ZBZ(sys)
980 slv0_system_t sys;
981 /**
982 *** Updates the reduced hessian matrix.
983 **/
984 {
985 mtx_coord_t nz;
986
987 if( sys->ZBZ.accurate ) return;
988
989 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
990 for( nz.col = 0; nz.col <= nz.row; nz.col++ ) {
991 int32 col, depr, depc;
992 col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
993 depr = mtx_col_to_org(sys->J.mtx,col);
994 col = nz.col+sys->J.reg.col.high+1-sys->ZBZ.order;
995 depc = mtx_col_to_org(sys->J.mtx,col);
996 sys->ZBZ.mtx[nz.row][nz.col] = (nz.row==nz.col ? 1.0 : 0.0);
997 col = sys->J.reg.col.low;
998 for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
999 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1000 if( set_is_member(sys->J.varpivots,ind) )
1001 sys->ZBZ.mtx[nz.row][nz.col] +=
1002 (-linsol_org_col_dependency(sys->J.sys,depr,ind))*
1003 (-linsol_org_col_dependency(sys->J.sys,depc,ind));
1004 }
1005 if( nz.row != nz.col )
1006 sys->ZBZ.mtx[nz.col][nz.row] =
1007 sys->ZBZ.mtx[nz.row][nz.col];
1008 }
1009 }
1010 if( OPTIMIZING(sys) ) {
1011 struct hessian_data *update;
1012 for( update=sys->B; update != NULL; update = update->next ) {
1013 mtx_coord_t nz;
1014 for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1015 int32 col, dep;
1016 col = nz.row + sys->J.reg.col.high + 1 - sys->ZBZ.order;
1017 dep = mtx_col_to_org(sys->J.mtx,col);
1018 sys->ZBZ.Zy[nz.row] = update->y.vec[col];
1019 sys->ZBZ.ZBs[nz.row] = update->Bs.vec[col];
1020 col = sys->J.reg.col.low;
1021 for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
1022 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1023 if( set_is_member(sys->J.varpivots,ind) ) {
1024 sys->ZBZ.Zy[nz.row] += update->y.vec[col]*
1025 (-linsol_org_col_dependency(sys->J.sys,dep,ind));
1026 sys->ZBZ.ZBs[nz.row] += update->Bs.vec[col]*
1027 (-linsol_org_col_dependency(sys->J.sys,dep,ind));
1028 }
1029 }
1030 for( nz.col=0; nz.col <= nz.row; nz.col++ ) {
1031 sys->ZBZ.mtx[nz.row][nz.col] += update->ys > 0.0 ?
1032 sys->ZBZ.Zy[nz.row]*sys->ZBZ.Zy[nz.col]/update->ys : 0.0;
1033 sys->ZBZ.mtx[nz.row][nz.col] -= update->sBs > 0.0 ?
1034 sys->ZBZ.ZBs[nz.row]*sys->ZBZ.ZBs[nz.col]/update->sBs : 0.0;
1035 if( nz.row != nz.col )
1036 sys->ZBZ.mtx[nz.col][nz.row] =
1037 sys->ZBZ.mtx[nz.row][nz.col];
1038 }
1039 }
1040 }
1041 }
1042 sys->ZBZ.accurate = TRUE;
1043 #if DEBUG
1044 FPRINTF(LIF(sys),"\nReduced Hessian: \n");
1045 debug_out_hessian(LIF(sys),sys);
1046 #endif
1047 }
1048
1049
1050 static void calc_rhs(sys,vec,scalar,transpose)
1051 slv0_system_t sys;
1052 struct vector_data *vec;
1053 real64 scalar;
1054 boolean transpose;
1055 /**
1056 *** Calculates just the jacobian RHS. This function should be used to
1057 *** supplement calculation of the jacobian. The vector vec must
1058 *** already be calculated and scaled so as to simply be added to the
1059 *** rhs. Caller is responsible for initially zeroing the rhs vector.
1060 **/
1061 {
1062 if( transpose ) { /* vec is indexed by col */
1063 int32 col;
1064 for( col=vec->rng->low; col<=vec->rng->high; col++ )
1065 sys->J.rhs[mtx_col_to_org(sys->J.mtx,col)] +=
1066 scalar*vec->vec[col];
1067
1068 } else { /* vec is indexed by row */
1069 int32 row;
1070 for( row=vec->rng->low; row<=vec->rng->high; row++ )
1071 sys->J.rhs[mtx_row_to_org(sys->J.mtx,row)] +=
1072 scalar*vec->vec[row];
1073 }
1074 linsol_rhs_was_changed(sys->J.sys,sys->J.rhs);
1075 }
1076
1077
1078 static void calc_multipliers(sys)
1079 slv0_system_t sys;
1080 /**
1081 *** Computes the lagrange multipliers for the equality constraints.
1082 **/
1083 {
1084 if( sys->multipliers.accurate )
1085 return;
1086
1087 if ( !OPTIMIZING(sys) ) {
1088 zero_vector(&(sys->multipliers));
1089 sys->multipliers.norm2 = 0.0;
1090 } else {
1091 linsol_system_t lsys = sys->J.sys;
1092 int32 row;
1093 sys->J.rhs = linsol_get_rhs(lsys,0);
1094 mem_zero_byte_cast(sys->J.rhs,0.0,sys->cap*sizeof(real64));
1095 calc_rhs(sys, &(sys->gradient), -1.0, TRUE );
1096 linsol_solve(lsys,sys->J.rhs);
1097 row = sys->multipliers.rng->low;
1098 for( ; row <= sys->multipliers.rng->high; row++ ) {
1099 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1100 sys->multipliers.vec[row] = linsol_var_value
1101 (lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row));
1102 rel_set_multiplier(rel,sys->multipliers.vec[row]*
1103 sys->weights.vec[row]);
1104
1105 }
1106 if (sys->iarray[SP0_SAVLIN]) {
1107 FILE *ldat;
1108 int32 ov;
1109 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1110 ldat=fopen(savlinfilename,"w");
1111 FPRINTF(ldat,
1112 "================= multipliersrhs (orgcoled) itn %d =====\n",
1113 sys->s.iteration);
1114 debug_write_array(ldat,sys->J.rhs,sys->cap);
1115 FPRINTF(ldat,
1116 "================= multipliers (orgrowed) ============\n");
1117 for(ov=0 ; ov < sys->cap; ov++ )
1118 FPRINTF(ldat,"%.20g\n",linsol_var_value(lsys,sys->J.rhs,ov));
1119 fclose(ldat);
1120 }
1121 square_norm( &(sys->multipliers) );
1122 }
1123 sys->multipliers.accurate = TRUE;
1124 #if DEBUG
1125 FPRINTF(LIF(sys),"Multipliers: ");
1126 debug_out_vector(LIF(sys),sys,&(sys->multipliers));
1127 #endif
1128 }
1129
1130
1131 static void calc_stationary(sys)
1132 slv0_system_t sys;
1133 /**
1134 *** Computes the gradient of the lagrangian which
1135 *** should be zero at the optimum solution.
1136 **/
1137 {
1138 if( sys->stationary.accurate )
1139 return;
1140
1141 if ( !OPTIMIZING(sys) ) {
1142 zero_vector(&(sys->stationary));
1143 sys->stationary.norm2 = 0.0;
1144 } else {
1145 int32 col;
1146 matrix_product(sys->J.mtx, &(sys->multipliers),
1147 &(sys->stationary), 1.0, TRUE);
1148 col = sys->stationary.rng->low;
1149 for( ; col <= sys->stationary.rng->high; col++ )
1150 sys->stationary.vec[col] += sys->gradient.vec[col];
1151 square_norm( &(sys->stationary) );
1152 }
1153 sys->stationary.accurate = TRUE;
1154 #if DEBUG
1155 FPRINTF(LIF(sys),"Stationary: ");
1156 debug_out_vector(LIF(sys),sys,&(sys->stationary));
1157 #endif
1158 }
1159
1160
1161 static void calc_gamma(sys)
1162 slv0_system_t sys;
1163 /**
1164 *** Calculate the gamma vector.
1165 **/
1166 {
1167 if( sys->gamma.accurate )
1168 return;
1169
1170 matrix_product(sys->J.mtx, &(sys->residuals),
1171 &(sys->gamma), -1.0, TRUE);
1172 square_norm( &(sys->gamma) );
1173 sys->gamma.accurate = TRUE;
1174
1175 #if DEBUG
1176 FPRINTF(LIF(sys),"Gamma: ");
1177 debug_out_vector(LIF(sys),sys,&(sys->gamma));
1178 #endif
1179 }
1180
1181 static void calc_Jgamma(sys)
1182 slv0_system_t sys;
1183 /**
1184 *** Calculate the Jgamma vector.
1185 **/
1186 {
1187 if( sys->Jgamma.accurate )
1188 return;
1189
1190 matrix_product(sys->J.mtx, &(sys->gamma),
1191 &(sys->Jgamma), 1.0, FALSE);
1192 square_norm( &(sys->Jgamma) );
1193 sys->Jgamma.accurate = TRUE;
1194 #if DEBUG
1195 FPRINTF(LIF(sys),"Jgamma: ");
1196 debug_out_vector(LIF(sys),sys,&(sys->Jgamma));
1197 #endif
1198 }
1199
1200
1201 static void calc_newton(sys)
1202 slv0_system_t sys;
1203 /**
1204 *** Computes a step to solve the linearized equations.
1205 **/
1206 {
1207 linsol_system_t lsys = sys->J.sys;
1208 int32 col;
1209
1210 if( sys->newton.accurate )
1211 return;
1212
1213 sys->J.rhs = linsol_get_rhs(lsys,1);
1214 mem_zero_byte_cast(sys->J.rhs,0.0,sys->cap*sizeof(real64));
1215 calc_rhs(sys, &(sys->residuals), -1.0, FALSE);
1216 linsol_solve(lsys,sys->J.rhs);
1217 col = sys->newton.rng->low;
1218 for( ; col <= sys->newton.rng->high; col++ )
1219 sys->newton.vec[col] = linsol_var_value
1220 (lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1221 if (sys->iarray[SP0_SAVLIN]) {
1222 FILE *ldat;
1223 int32 ov;
1224 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1225 ldat=fopen(savlinfilename,"w");
1226 FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1227 sys->s.iteration);
1228 debug_write_array(ldat,sys->J.rhs,sys->cap);
1229 FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1230 for(ov=0 ; ov < sys->cap; ov++ )
1231 FPRINTF(ldat,"%.20g\n",linsol_var_value(lsys,sys->J.rhs,ov));
1232 fclose(ldat);
1233 }
1234 square_norm( &(sys->newton) );
1235 sys->newton.accurate = TRUE;
1236 #if DEBUG
1237 FPRINTF(LIF(sys),"Newton: ");
1238 debug_out_vector(LIF(sys),sys,&(sys->newton));
1239 #endif
1240 }
1241
1242
1243 static void calc_Bnewton(sys)
1244 slv0_system_t sys;
1245 /**
1246 *** Computes an update to the product B and newton.
1247 **/
1248 {
1249 if( sys->Bnewton.accurate )
1250 return;
1251
1252 if ( !OPTIMIZING(sys) ) {
1253 zero_vector(&(sys->Bnewton));
1254 sys->Bnewton.norm2 = 0.0;
1255 } else {
1256 struct hessian_data *update;
1257 copy_vector(&(sys->newton),&(sys->Bnewton));
1258 for( update=sys->B; update != NULL; update = update->next ) {
1259 int32 col;
1260 real64 yn = inner_product( &(update->y),&(sys->newton) );
1261 real64 sBn = inner_product( &(update->Bs),&(sys->newton) );
1262 col = sys->Bnewton.rng->low;
1263 for( ; col <= sys->Bnewton.rng->high; col++ ) {
1264 sys->Bnewton.vec[col] += update->ys > 0.0 ?
1265 (update->y.vec[col])*yn/update->ys : 0.0;
1266 sys->Bnewton.vec[col] -= update->sBs > 0.0 ?
1267 (update->Bs.vec[col])*sBn/update->sBs : 0.0;
1268 }
1269 }
1270 square_norm( &(sys->Bnewton) );
1271 }
1272 sys->Bnewton.accurate = TRUE;
1273 #if DEBUG
1274 FPRINTF(LIF(sys),"Bnewton: ");
1275 debug_out_vector(LIF(sys),sys,&(sys->Bnewton));
1276 #endif
1277 }
1278
1279
1280 static void calc_nullspace(sys)
1281 slv0_system_t sys;
1282 /**
1283 *** Calculate the nullspace move.
1284 **/
1285 {
1286 if( sys->nullspace.accurate )
1287 return;
1288
1289 if( !OPTIMIZING(sys) ) {
1290 zero_vector(&(sys->nullspace));
1291 sys->nullspace.norm2 = 0.0;
1292 } else {
1293 mtx_coord_t nz;
1294 zero_vector(&(sys->nullspace));
1295 for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1296 int32 col, dep, ndx;
1297 col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
1298 dep = mtx_col_to_org(sys->J.mtx,col);
1299 sys->nullspace.vec[col] = -sys->stationary.vec[col] -
1300 sys->Bnewton.vec[col];
1301 ndx = sys->J.reg.col.low;
1302 for( ; ndx <= sys->J.reg.col.high - sys->ZBZ.order; ndx++ ) {
1303 int32 ind = mtx_col_to_org(sys->J.mtx,ndx);
1304 if( set_is_member(sys->J.varpivots,ind) )
1305 sys->nullspace.vec[col] -=
1306 (sys->stationary.vec[ndx] + sys->Bnewton.vec[ndx])*
1307 (-linsol_org_col_dependency(sys->J.sys,dep,ind));
1308 }
1309 }
1310 /**
1311 *** Must invert ZBZ first. It's symmetric so
1312 *** can find Cholesky factors. Essentially, find
1313 *** the "square root" of the matrix such that
1314 ***
1315 *** T
1316 *** L U = U U = ZBZ, where U is an upper triangular
1317 *** matrix.
1318 **/
1319 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1320 for( nz.col = nz.row; nz.col < sys->ZBZ.order; nz.col++ ) {
1321 int32 col;
1322 for( col = 0; col < nz.row; col++ )
1323 sys->ZBZ.mtx[nz.row][nz.col] -=
1324 sys->ZBZ.mtx[nz.row][col]*
1325 sys->ZBZ.mtx[col][nz.col];
1326 if( nz.row == nz.col )
1327 sys->ZBZ.mtx[nz.row][nz.col] =
1328 calc_sqrt_D0(sys->ZBZ.mtx[nz.row][nz.col]);
1329 else {
1330 sys->ZBZ.mtx[nz.row][nz.col] /=
1331 sys->ZBZ.mtx[nz.row][nz.row];
1332 sys->ZBZ.mtx[nz.col][nz.row] =
1333 sys->ZBZ.mtx[nz.row][nz.col];
1334 }
1335 }
1336 }
1337 #if DEBUG
1338 FPRINTF(LIF(sys),"\nInverse Reduced Hessian: \n");
1339 debug_out_hessian(LIF(sys),sys);
1340 #endif
1341 /**
1342 *** forward substitute
1343 **/
1344 for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1345 int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1346 for( nz.col = 0; nz.col < nz.row; nz.col++ ) {
1347 sys->nullspace.vec[nz.row+offset] -=
1348 sys->nullspace.vec[nz.col+offset]*
1349 sys->ZBZ.mtx[nz.row][nz.col];
1350 }
1351 sys->nullspace.vec[nz.row+offset] /=
1352 sys->ZBZ.mtx[nz.row][nz.row];
1353 }
1354
1355 /**
1356 *** backward substitute
1357 **/
1358 for( nz.row = sys->ZBZ.order-1; nz.row >= 0; nz.row-- ) {
1359 int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1360 for( nz.col = nz.row+1; nz.col < sys->ZBZ.order; nz.col++ ) {
1361 sys->nullspace.vec[nz.row+offset] -=
1362 sys->nullspace.vec[nz.col+offset]*
1363 sys->ZBZ.mtx[nz.row][nz.col];
1364 }
1365 sys->nullspace.vec[nz.row+offset] /=
1366 sys->ZBZ.mtx[nz.row][nz.row];
1367 }
1368 square_norm( &(sys->nullspace) );
1369 }
1370 sys->nullspace.accurate = TRUE;
1371 #if DEBUG
1372 FPRINTF(LIF(sys),"Nullspace: ");
1373 debug_out_vector(LIF(sys),sys,&(sys->nullspace));
1374 #endif
1375 }
1376
1377 static void calc_varstep1(sys)
1378 slv0_system_t sys;
1379 /**
1380 *** Calculate the 1st order descent direction for phi
1381 *** in the variables.
1382 **/
1383 {
1384 if( sys->varstep1.accurate )
1385 return;
1386
1387 if( !OPTIMIZING(sys) ) {
1388 copy_vector(&(sys->gamma),&(sys->varstep1));
1389 sys->varstep1.norm2 = sys->gamma.norm2;
1390 } else {
1391 int32 col;
1392 col = sys->varstep1.rng->low;
1393 for( ; col <= sys->varstep1.rng->high; col++ )
1394 sys->varstep1.vec[col] = sys->p.rho*sys->gamma.vec[col] -
1395 sys->stationary.vec[col];
1396 square_norm( &(sys->varstep1) );
1397 }
1398 sys->varstep1.accurate = TRUE;
1399 #if DEBUG
1400 FPRINTF(LIF(sys),"Varstep1: ");
1401 debug_out_vector(LIF(sys),sys,&(sys->varstep1));
1402 #endif
1403 }
1404
1405
1406 static void calc_Bvarstep1(sys)
1407 slv0_system_t sys;
1408 /**
1409 *** Computes an update to the product B and varstep1.
1410 **/
1411 {
1412 if( sys->Bvarstep1.accurate )
1413 return;
1414
1415 if ( !OPTIMIZING(sys) ) {
1416 zero_vector(&(sys->Bvarstep1));
1417 sys->Bvarstep1.norm2 = 0.0;
1418 } else {
1419 struct hessian_data *update;
1420 copy_vector(&(sys->varstep1),&(sys->Bvarstep1));
1421 for( update=sys->B; update != NULL; update = update->next ) {
1422 int32 col;
1423 real64 yv = inner_product( &(update->y),&(sys->varstep1) );
1424 real64 sBv = inner_product( &(update->Bs),&(sys->varstep1) );
1425 col = sys->Bvarstep1.rng->low;
1426 for( ; col <= sys->Bvarstep1.rng->high; col++ ) {
1427 sys->Bvarstep1.vec[col] += update->ys > 0.0 ?
1428 (update->y.vec[col])*yv/update->ys : 0.0;
1429 sys->Bvarstep1.vec[col] -= update->sBs > 0.0 ?
1430 (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1431 }
1432 }
1433 square_norm( &(sys->Bvarstep1) );
1434 }
1435 sys->Bvarstep1.accurate = TRUE;
1436 #if DEBUG
1437 FPRINTF(LIF(sys),"Bvarstep1: ");
1438 debug_out_vector(LIF(sys),sys,&(sys->Bvarstep1));
1439 #endif
1440 }
1441
1442
1443 static void calc_varstep2(sys)
1444 slv0_system_t sys;
1445 /**
1446 *** Calculate the 2nd order descent direction for phi
1447 *** in the variables.
1448 **/
1449 {
1450 if( sys->varstep2.accurate )
1451 return;
1452
1453 if( !OPTIMIZING(sys) ) {
1454 copy_vector(&(sys->newton),&(sys->varstep2));
1455 sys->varstep2.norm2 = sys->newton.norm2;
1456 } else {
1457 int32 col;
1458 col = sys->varstep2.rng->low;
1459 for( ; col <= sys->varstep2.rng->high - sys->ZBZ.order ; ++col ) {
1460 int32 dep;
1461 int32 ind = mtx_col_to_org(sys->J.mtx,col);
1462 sys->varstep2.vec[col] = sys->newton.vec[col];
1463 if( set_is_member(sys->J.varpivots,ind) ) {
1464 dep = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1465 for( ; dep <= sys->varstep2.rng->high; dep++ )
1466 sys->varstep2.vec[col] += sys->nullspace.vec[dep]*
1467 (-linsol_org_col_dependency(sys->J.sys,dep,ind));
1468 }
1469 }
1470 col = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1471 for( ; col <= sys->varstep2.rng->high; ++col )
1472 sys->varstep2.vec[col] = sys->nullspace.vec[col] +
1473 sys->newton.vec[col];
1474 square_norm( &(sys->varstep2) );
1475 }
1476 sys->varstep2.accurate = TRUE;
1477 #if DEBUG
1478 FPRINTF(LIF(sys),"Varstep2: ");
1479 debug_out_vector(LIF(sys),sys,&(sys->varstep2));
1480 #endif
1481 }
1482
1483
1484 static void calc_Bvarstep2(sys)
1485 slv0_system_t sys;
1486 /**
1487 *** Computes an update to the product B and varstep2.
1488 **/
1489 {
1490 if( sys->Bvarstep2.accurate )
1491 return;
1492
1493 if ( !OPTIMIZING(sys) ) {
1494 zero_vector(&(sys->Bvarstep2));
1495 sys->Bvarstep2.norm2 = 0.0;
1496 } else {
1497 struct hessian_data *update;
1498 copy_vector(&(sys->varstep2),&(sys->Bvarstep2));
1499 for( update=sys->B; update != NULL; update = update->next ) {
1500 int32 col;
1501 real64 yv = inner_product( &(update->y),&(sys->varstep2) );
1502 real64 sBv = inner_product( &(update->Bs),&(sys->varstep2) );
1503 col = sys->Bvarstep2.rng->low;
1504 for( ; col <= sys->Bvarstep2.rng->high; col++ ) {
1505 sys->Bvarstep2.vec[col] += update->ys > 0.0 ?
1506 (update->y.vec[col])*yv/update->ys : 0.0;
1507 sys->Bvarstep2.vec[col] -= update->sBs > 0.0 ?
1508 (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1509 }
1510 }
1511 square_norm( &(sys->Bvarstep2) );
1512 }
1513 sys->Bvarstep2.accurate = TRUE;
1514 #if DEBUG
1515 FPRINTF(LIF(sys),"Bvarstep2: ");
1516 debug_out_vector(LIF(sys),sys,&(sys->Bvarstep2));
1517 #endif
1518 }
1519
1520
1521 static void calc_mulstep1(sys)
1522 slv0_system_t sys;
1523 /**
1524 *** Calculate the negative gradient direction of phi in the
1525 *** multipliers.
1526 **/
1527 {
1528 if( sys->mulstep1.accurate )
1529 return;
1530
1531 if( !OPTIMIZING(sys) ) {
1532 zero_vector(&(sys->mulstep1));
1533 sys->mulstep1.norm2 = 0.0;
1534 } else {
1535 int32 row;
1536 row = sys->mulstep1.rng->low;
1537 for( ; row <= sys->mulstep1.rng->high; row++ )
1538 sys->mulstep1.vec[row] = -sys->residuals.vec[row];
1539 square_norm( &(sys->mulstep1) );
1540 }
1541 sys->mulstep1.accurate = TRUE;
1542 #if DEBUG
1543 FPRINTF(LIF(sys),"Mulstep1: ");
1544 debug_out_vector(LIF(sys),sys,&(sys->mulstep1));
1545 #endif
1546 }
1547
1548
1549 static void calc_mulstep2(sys)
1550 slv0_system_t sys;
1551 /**
1552 *** Calculate the mulstep2 direction of phi in the
1553 *** multipliers.
1554 **/
1555 {
1556 if( sys->mulstep2.accurate )
1557 return;
1558
1559 if( !OPTIMIZING(sys) ) {
1560 zero_vector(&(sys->mulstep2));
1561 sys->mulstep2.norm2 = 0.0;
1562 } else {
1563 linsol_system_t lsys = sys->J.sys;
1564 int32 row;
1565 sys->J.rhs = linsol_get_rhs(lsys,2);
1566 mem_zero_byte_cast(sys->J.rhs,0.0,sys->cap*sizeof(real64));
1567 calc_rhs(sys, &(sys->Bvarstep2), -1.0, TRUE);
1568 calc_rhs(sys, &(sys->stationary), -1.0, TRUE);
1569 linsol_solve(lsys,sys->J.rhs);
1570 row = sys->mulstep2.rng->low;
1571 for( ; row <= sys->mulstep2.rng->high; row++ )
1572 sys->mulstep2.vec[row] = linsol_var_value
1573 (lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row));
1574 if (sys->iarray[SP0_SAVLIN]) {
1575 FILE *ldat;
1576 int32 ov;
1577 sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1578 ldat=fopen(savlinfilename,"w");
1579 FPRINTF(ldat,
1580 "================= mulstep2rhs (orgcoled) itn %d =======\n",
1581 sys->s.iteration);
1582 debug_write_array(ldat,sys->J.rhs,sys->cap);
1583 FPRINTF(ldat,
1584 "================= mulstep2vars (orgrowed) ============\n");
1585 for(ov=0 ; ov < sys->cap; ov++ )
1586 FPRINTF(ldat,"%.20g\n",linsol_var_value(lsys,sys->J.rhs,ov));
1587 fclose(ldat);
1588 }
1589 square_norm( &(sys->mulstep2) );
1590 }
1591 sys->mulstep2.accurate = TRUE;
1592 #if DEBUG
1593 FPRINTF(LIF(sys),"Mulstep2: ");
1594 debug_out_vector(LIF(sys),sys,&(sys->mulstep2));
1595 #endif
1596 }
1597
1598
1599 static void calc_phi(sys)
1600 slv0_system_t sys;
1601 /**
1602 *** Computes the global minimizing function Phi.
1603 **/
1604 {
1605 if( !OPTIMIZING(sys) )
1606 sys->phi = 0.5*sys->residuals.norm2;
1607 else {
1608 sys->phi = sys->objective;
1609 sys->phi += inner_product( &(sys->multipliers),&(sys->residuals) );
1610 sys->phi += 0.5*sys->p.rho*sys->residuals.norm2;
1611 }
1612 }
1613
1614
1615 /**
1616 *** OK. Here's where we compute the actual step to be taken. It will
1617 *** be some linear combination of the 1st order and 2nd order steps.
1618 **/
1619
1620 typedef real64 sym_2x2_t[3]; /* Stores symmetric 2x2 matrices */
1621
1622 struct parms_t {
1623 real64 low,high,guess; /* Used to search for parameter */
1624 };
1625
1626 struct calc_step_vars {
1627 sym_2x2_t coef1, coef2;
1628 real64 rhs[2]; /* RHS for 2x2 system */
1629 struct parms_t parms;
1630 real64 alpha1, alpha2;
1631 real64 error; /* Error between step norm and sys->maxstep */
1632 };
1633
1634 static void calc_2x2_system(sys,vars)
1635 slv0_system_t sys;
1636 struct calc_step_vars *vars;
1637 /**
1638 *** Calculates 2x2 system (coef1,coef2,rhs).
1639 **/
1640 {
1641 vars->coef1[0] = (2.0*sys->phi/sys->newton.norm2)*
1642 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1643 vars->coef1[1] = 1.0;
1644 vars->coef1[2] = (sys->Jgamma.norm2/sys->gamma.norm2)*
1645 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1646
1647 vars->coef2[0] = 1.0;
1648 vars->coef2[1] = 2.0*sys->phi/
1649 calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
1650 vars->coef2[2] = 1.0;
1651
1652 vars->rhs[0] = 2.0*sys->phi/
1653 sys->maxstep/calc_sqrt_D0(sys->gamma.norm2);
1654 vars->rhs[1] = calc_sqrt_D0(sys->newton.norm2)/sys->maxstep;
1655 }
1656
1657 static void coefs_from_parm(sys,vars)
1658 slv0_system_t sys;
1659 struct calc_step_vars *vars;
1660 /**
1661 *** Determines alpha1 and alpha2 from the parameter (guess).
1662 **/
1663 {
1664 #define DETZERO (1e-8)
1665
1666 sym_2x2_t coef; /* Actual coefficient matrix */
1667 real64 det; /* Determinant of coefficient matrix */
1668 int i;
1669
1670 for( i=0 ; i<3 ; ++i ) coef[i] =
1671 vars->coef1[i] + vars->parms.guess * vars->coef2[i];
1672 det = coef[0]*coef[2] - coef[1]*coef[1];
1673 if( det < 0.0 )
1674 FPRINTF(MIF(sys),"%-40s ---> %g\n",
1675 " Unexpected negative determinant!",det);
1676 if( det <= DETZERO ) {
1677 /**
1678 *** varstep2 and varstep1 are essentially parallel:
1679 *** adjust length of either
1680 **/
1681 vars->alpha2 = 0.0;
1682 vars->alpha1 = 1.0;
1683 } else {
1684 vars->alpha2 = (vars->rhs[0]*coef[2] - vars->rhs[1]*coef[1])/det;
1685 vars->alpha1 = (vars->rhs[1]*coef[0] - vars->rhs[0]*coef[1])/det;
1686 }
1687 #undef DETZERO
1688 }
1689
1690 static real64 step_norm2(sys,vars)
1691 slv0_system_t sys;
1692 struct calc_step_vars *vars;
1693 /**
1694 *** Computes step vector length based on 1st order and 2nd order
1695 *** vectors and their coefficients.
1696 **/
1697 {
1698 return sys->maxstep*sys->maxstep*
1699 (vars->alpha2 * vars->alpha2 +
1700 vars->alpha2 * vars->alpha1 * sys->phi/
1701 calc_sqrt_D0(sys->varstep2.norm2 + sys->mulstep2.norm2)/
1702 calc_sqrt_D0(sys->varstep1.norm2 + sys->mulstep1.norm2) +
1703 vars->alpha1 * vars->alpha1);
1704 }
1705
1706 static void adjust_parms(sys,vars)
1707 slv0_system_t sys;
1708 struct calc_step_vars *vars;
1709 /**
1710 *** Re-guesses the parameters based on
1711 *** step size vs. target value.
1712 **/
1713 {
1714 vars->error = (calc_sqrt_D0(step_norm2(sys,vars))/sys->maxstep) - 1.0;
1715 if( vars->error > 0.0 ) {
1716 /* Increase parameter (to decrease step length) */
1717 vars->parms.low = vars->parms.guess;
1718 vars->parms.guess = (vars->parms.high>3.0*vars->parms.guess)
1719 ? 2.0*vars->parms.guess
1720 : 0.5*(vars->parms.low + vars->parms.high);
1721 } else {
1722 /* Decrease parameter (to increase step norm) */
1723 vars->parms.high = vars->parms.guess;
1724 vars->parms.guess = 0.5*(vars->parms.low + vars->parms.high);
1725 }
1726 }
1727
1728 static void compute_step(sys,vars)
1729 slv0_system_t sys;
1730 struct calc_step_vars *vars;
1731 /**
1732 *** Computes the step based on the coefficients in vars.
1733 **/
1734 {
1735 int32 row,col;
1736 real64 tot1_norm2, tot2_norm2;
1737
1738 tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
1739 tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
1740 if( !sys->varstep.accurate ) {
1741 for( col=sys->varstep.rng->low ; col<=sys->varstep.rng->high ; ++col )
1742 if( (vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
1743 sys->varstep.vec[col] = sys->maxstep *
1744 sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2);
1745 } else if( (vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
1746 sys->varstep.vec[col] = sys->maxstep *
1747 sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2);
1748 } else if( (vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
1749 sys->varstep.vec[col] = sys->maxstep*
1750 (
1751 vars->alpha2*sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2) +
1752 vars->alpha1*sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2)
1753 );
1754 }
1755 sys->varstep.accurate = TRUE;
1756 }
1757 if( !sys->mulstep.accurate ) {
1758 for( row=sys->mulstep.rng->low ; row<=sys->mulstep.rng->high ; ++row )
1759 if( (vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
1760 sys->mulstep.vec[row] = sys->maxstep *
1761 sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2);
1762 } else if( (vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
1763 sys->mulstep.vec[row] = sys->maxstep *
1764 sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2);
1765 } else if( (vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
1766 sys->mulstep.vec[row] = sys->maxstep*
1767 (
1768 vars->alpha2*sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2) +
1769 vars->alpha1*sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2)
1770 );
1771 }
1772 sys->mulstep.accurate = TRUE;
1773 }
1774 #if DEBUG
1775 FPRINTF(LIF(sys),"Varstep: ");
1776 debug_out_vector(LIF(sys),sys,&(sys->varstep));
1777 FPRINTF(LIF(sys),"Mulstep: ");
1778 debug_out_vector(LIF(sys),sys,&(sys->mulstep));
1779 #endif
1780 }
1781
1782
1783 static void calc_step(sys, minor)
1784 slv0_system_t sys;
1785 int minor;
1786 /**
1787 *** Calculates step vector, based on sys->maxstep, and the varstep2/
1788 *** varstep1 and mulstep2/mulstep1 vectors. Nothing is assumed to be
1789 *** calculated, except the weights and the jacobian (scaled). Also,
1790 *** the step is not checked for legitimacy.
1791 *** NOTE: the step is scaled.
1792 **/
1793 {
1794 #define STEPSIZEERR_MAX 1e-4
1795 /* Step size must be determined this precisely */
1796 #define PARMRNG_MIN 1e-12
1797 /* OR parameter range must be this narrow to exit */
1798
1799 struct calc_step_vars vars;
1800 real64 tot1_norm2, tot2_norm2;
1801
1802 if( sys->varstep.accurate && sys->mulstep.accurate )
1803 return;
1804 if (sys->p.output.less_important) {
1805 FPRINTF(LIF(sys),"\n%-40s ---> %d\n", " Step trial",minor);
1806 }
1807
1808 tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
1809 tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
1810 if( (tot1_norm2 == 0.0) && (tot2_norm2 == 0.0) ) {
1811 /* Take no step at all */
1812 vars.alpha1 = 0.0;
1813 vars.alpha2 = 0.0;
1814 sys->maxstep = 0.0;
1815 sys->varstep.norm2 = 0.0;
1816 sys->mulstep.norm2 = 0.0;
1817
1818 } else if( (tot2_norm2 > 0.0) && OPTIMIZING(sys) ) {
1819 /* Stay in varstep2 direction */
1820 vars.alpha1 = 0.0;
1821 vars.alpha2 = 1.0;
1822 sys->maxstep = MIN(sys->maxstep,calc_sqrt_D0(tot2_norm2));
1823 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
1824 sys->varstep2.norm2/tot2_norm2;
1825 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
1826 sys->mulstep2.norm2/tot2_norm2;
1827
1828 } else if( (tot2_norm2>0.0)&&(calc_sqrt_D0(tot2_norm2)<=sys->maxstep) ) {
1829 /* Attempt step in varstep2 direction */
1830 vars.alpha1 = 0.0;
1831 vars.alpha2 = 1.0;
1832 sys->maxstep = calc_sqrt_D0(tot2_norm2);
1833 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
1834 sys->varstep2.norm2/tot2_norm2;
1835 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
1836 sys->mulstep2.norm2/tot2_norm2;
1837
1838 } else if( (tot2_norm2==0.0 || sys->s.block.current_size==1) &&
1839 (tot1_norm2 > 0.0) ) {
1840 /* Attempt step in varstep1 direction */
1841 vars.alpha1 = 1.0;
1842 vars.alpha2 = 0.0;
1843 if ( (sys->gamma.norm2/sys->Jgamma.norm2)*
1844 calc_sqrt_D0(sys->gamma.norm2) <= sys->maxstep )
1845 sys->maxstep = (sys->gamma.norm2/sys->Jgamma.norm2)*
1846 calc_sqrt_D0(sys->gamma.norm2);
1847 sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
1848 sys->varstep1.norm2/tot1_norm2;
1849 sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
1850 sys->mulstep1.norm2/tot1_norm2;
1851
1852 } else {
1853 /* Attempt step in varstep1-varstep2 direction */
1854 vars.parms.low = 0.0;
1855 vars.parms.high = MAXDOUBLE;
1856 vars.parms.guess = 1.0;
1857 calc_2x2_system(sys,&vars);
1858 do {
1859 coefs_from_parm(sys, &vars);
1860 adjust_parms(sys, &vars);
1861 } while( fabs(vars.error) > STEPSIZEERR_MAX &&
1862 vars.parms.high - vars.parms.low > PARMRNG_MIN );
1863 if (sys->p.output.less_important) {
1864 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1865 " parameter high", vars.parms.high);
1866 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1867 " parameter low", vars.parms.low);
1868 FPRINTF(LIF(sys),"%-40s ---> %g\n",
1869 " Error in step length", vars.error);
1870 }
1871 sys->varstep.norm2 = step_norm2(sys, &vars);
1872 sys->mulstep.norm2 = 0.0;
1873 }
1874
1875 if (sys->p.output.less_important) {
1876 FPRINTF(LIF(sys),"%-40s ---> %g\n", " Alpha1 coefficient (normalized)",
1877 vars.alpha1);
1878 FPRINTF(LIF(sys),"%-40s ---> %g\n", " Alpha2 coefficient (normalized)",
1879 vars.alpha2);
1880 }
1881 compute_step(sys,&vars);
1882 return;
1883
1884 #undef STEPSIZEERR_MAX
1885 #undef PARMRNG_MIN
1886 }
1887
1888
1889
1890 /**
1891 *** Variable values maintenance
1892 *** ---------------------------
1893 *** restore_variables(sys)
1894 *** coef = required_coef_to_stay_inbounds(sys)
1895 *** apply_step(sys,coef)
1896 *** step_accepted(sys)
1897 *** change_maxstep(sys,maxstep)
1898 **/
1899
1900 static void restore_variables(sys)
1901 slv0_system_t sys;
1902 /**
1903 *** Restores the values of the variables before applying
1904 *** a step.
1905 **/
1906 {
1907 int32 col;
1908
1909 for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
1910 struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1911 var_set_value(var,sys->variables.vec[col]*sys->nominals.vec[col]);
1912 }
1913 }
1914
1915
1916 static real64 required_coef_to_stay_inbounds(sys)
1917 slv0_system_t sys;
1918 /**
1919 *** Calculates the maximum fraction of the step which can be
1920 *** taken without going out of bounds. If the entire step can be
1921 *** taken, 1.0 is returned. Otherwise a value less than 1 is
1922 *** returned. It is assumed that the current variable values
1923 *** are within their bounds. The step must be calculated.
1924 **/
1925 {
1926 real64 mincoef;
1927 int32 col;
1928
1929 if( sys->p.ignore_bounds )
1930 return(1.0);
1931
1932 mincoef = 1.0;
1933 for( col=sys->varstep.rng->low; col <= sys->varstep.rng->high; col++ ) {
1934 struct var_variable *var;
1935 real64 coef,dx,val,bnd;
1936 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1937 coef = 1.0;
1938 dx = sys->varstep.vec[col] * sys->nominals.vec[col];
1939 bnd = var_upper_bound(var);
1940 if( (val=var_value(var)) + dx > bnd )
1941 coef = MIN((bnd-val)/dx, 1.0);
1942 bnd = var_lower_bound(var);
1943 if( val + dx < bnd )
1944 coef = MIN((bnd-val)/dx, 1.0);
1945 if( coef < mincoef )
1946 mincoef = coef;
1947 }
1948 return(mincoef);
1949 }
1950
1951
1952 #define TRUNCATE FALSE /* override projection */
1953 static void apply_step(sys)
1954 slv0_system_t sys;
1955 /**
1956 *** Adds sys->varstep to the variable values in block: projecting
1957 *** near bounds.
1958 **/
1959 {
1960 FILE *lif = LIF(sys);
1961 int nproj = 0;
1962 real64 bounds_coef = 1.0;
1963 int32 col;
1964
1965 if (TRUNCATE && (!sys->p.ignore_bounds))
1966 bounds_coef = required_coef_to_stay_inbounds(sys);
1967
1968 for( col=sys->varstep.rng->low; col <= sys->varstep.rng->high; col++ ) {
1969 struct var_variable *var;
1970 real64 dx,val,bnd;
1971 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1972 dx = sys->nominals.vec[col]*sys->varstep.vec[col];
1973 val = var_value(var);
1974 if (bounds_coef < 1.0) {
1975 dx = dx*TOWARD_BOUNDS*bounds_coef;
1976 sys->varstep.vec[col] = dx/sys->nominals.vec[col];
1977 } else {
1978 if( !sys->p.ignore_bounds ) {
1979 if( val + dx > (bnd=var_upper_bound(var)) ) {
1980 dx = TOWARD_BOUNDS*(bnd-val);
1981 sys->varstep.vec[col] = dx/sys->nominals.vec[col];
1982 if (sys->p.output.less_important) {
1983 FPRINTF(lif,"%-40s ---> ",
1984 " Variable projected to upper bound");
1985 print_var_name(lif,sys,var); PUTC('\n',lif);
1986 }
1987 ++nproj;
1988 } else if( val + dx < (bnd=var_lower_bound(var)) ) {
1989 dx = TOWARD_BOUNDS*(bnd-val);
1990 sys->varstep.vec[col] = dx/sys->nominals.vec[col];
1991 if (sys->p.output.less_important) {
1992 FPRINTF(lif,"%-40s ---> ",
1993 " Variable projected to lower bound");
1994 print_var_name(lif,sys,var); PUTC('\n',lif);
1995 }
1996 ++nproj;
1997 }
1998 }
1999 }
2000 var_set_value(var,val+dx);
2001 }
2002
2003 if( !sys->p.ignore_bounds ) {
2004 if (nproj > 0) {
2005 square_norm(&(sys->varstep));
2006 sys->progress = calc_sqrt_D0
2007 (calc_sqrt_D0((sys->varstep.norm2 + sys->mulstep.norm2)*
2008 (sys->varstep1.norm2 + sys->mulstep1.norm2)));
2009 if (sys->p.output.less_important) {
2010 FPRINTF(lif,"%-40s ---> %g\n", " Projected step length (scaled)",
2011 calc_sqrt_D0(sys->varstep.norm2 + sys->mulstep.norm2));
2012 FPRINTF(lif,"%-40s ---> %g\n",
2013 " Projected progress", sys->progress);
2014 }
2015 }
2016 if (bounds_coef < 1.0) {
2017 square_norm(&(sys->varstep));
2018 if (sys->p.output.less_important) {
2019 FPRINTF(lif,"%-40s ---> %g\n",
2020 " Truncated step length (scaled)",
2021 calc_sqrt_D0(sys->varstep.norm2 + sys->mulstep.norm2));
2022 }
2023 sys->progress = calc_sqrt_D0
2024 (calc_sqrt_D0((sys->varstep.norm2 + sys->mulstep.norm2)*
2025 (sys->varstep1.norm2 + sys->mulstep1.norm2)));
2026 (calc_sqrt_D0(sys->varstep.norm2*sys->varstep1.norm2));
2027 if (sys->p.output.less_important) {
2028 FPRINTF(lif,"%-40s ---> %g\n",
2029 " Truncated progress", sys->progress);
2030 }
2031 }
2032 }
2033
2034 /* Allow weighted residuals to be recalculated at new point */
2035 sys->residuals.accurate = FALSE;
2036
2037 return;
2038 }
2039 #undef TRUNCATE
2040
2041 static void step_accepted(sys)
2042 slv0_system_t sys;
2043 /**
2044 *** This function should be called when the step is accepted.
2045 **/
2046 {
2047 /* Maintain update status on jacobian and weights */
2048 if (--(sys->update.jacobian) <= 0)
2049 sys->J.accurate = FALSE;
2050
2051 sys->ZBZ.accurate = FALSE;
2052 sys->variables.accurate = FALSE;
2053 sys->gradient.accurate = FALSE;
2054 sys->multipliers.accurate = FALSE;
2055 sys->stationary.accurate = FALSE;
2056 sys->newton.accurate = FALSE;
2057 sys->Bnewton.accurate = FALSE;
2058 sys->nullspace.accurate = FALSE;
2059 sys->gamma.accurate = FALSE;
2060 sys->Jgamma.accurate = FALSE;
2061 sys->varstep1.accurate = FALSE;
2062 sys->Bvarstep1.accurate = FALSE;
2063 sys->varstep2.accurate = FALSE;
2064 sys->Bvarstep2.accurate = FALSE;
2065 sys->mulstep1.accurate = FALSE;
2066 sys->mulstep2.accurate = FALSE;
2067 sys->varstep.accurate = FALSE;
2068 sys->mulstep.accurate = FALSE;
2069
2070 if( !OPTIMIZING(sys) ) {
2071 sys->ZBZ.accurate = TRUE;
2072 sys->gradient.accurate = TRUE;
2073 sys->multipliers.accurate = TRUE;
2074 sys->stationary.accurate = TRUE;
2075 sys->Bnewton.accurate = TRUE;
2076 sys->nullspace.accurate = TRUE;
2077 sys->Bvarstep1.accurate = TRUE;
2078 sys->Bvarstep2.accurate = TRUE;
2079 }
2080 }
2081
2082 static void change_maxstep(sys,maxstep)
2083 slv0_system_t sys;
2084 real64 maxstep;
2085 /**
2086 *** This function changes sys->maxstep to the given number and should be
2087 *** called whenever sys->maxstep is to be changed.
2088 **/
2089 {
2090 sys->maxstep = maxstep;
2091 sys->varstep.accurate = FALSE;
2092 if( OPTIMIZING(sys) ) sys->mulstep.accurate = FALSE;
2093 }
2094
2095
2096
2097 /**
2098 *** Block routines
2099 *** --------------
2100 *** feas = block_feasible(sys)
2101 *** insure_bounds(sys,var)
2102 *** move_to_next_block(sys)
2103 *** find_next_unconverged_block(sys)
2104 **/
2105
2106 static boolean block_feasible(sys)
2107 slv0_system_t sys;
2108 /**
2109 *** Returns TRUE if the current block is feasible, FALSE otherwise.
2110 *** It is assumed that the residuals have been computed.
2111 **/
2112 {
2113 int32 row;
2114
2115 if( !sys->s.calc_ok )
2116 return(FALSE);
2117
2118 for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2119 struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2120 if( !rel_satisfied(rel) ) return FALSE;
2121 }
2122 return TRUE;
2123 }
2124
2125 static void insure_bounds(sys,var)
2126 slv0_system_t sys;
2127 struct var_variable *var;
2128 /**
2129 *** Insures that the variable value is within its bounds.
2130 **/
2131 {
2132 real64 val,low,high;
2133 FILE *mif = MIF(sys);
2134
2135 if( sys->p.ignore_bounds )
2136 return;
2137
2138 low = var_lower_bound(var);
2139 high = var_upper_bound(var);
2140 val = var_value(var);
2141 if( low > high ) {
2142 FPRINTF(mif,"Bounds for variable ");
2143 print_var_name(mif,sys,var);
2144 FPRINTF(mif," are inconsistent [%g,%g].\n",low,high);
2145 FPRINTF(mif,"Bounds will be swapped.\n");
2146 var_set_upper_bound(var, low);
2147 var_set_lower_bound(var, high);
2148 low = var_lower_bound(var);
2149 high = var_upper_bound(var);
2150 }
2151
2152 if( low > val ) {
2153 FPRINTF(mif,"Variable ");
2154 print_var_name(mif,sys,var);
2155 FPRINTF(mif," was initialized below its lower bound.\n");
2156 FPRINTF(mif,"It will be moved to its lower bound.\n");
2157 var_set_value(var, low);
2158 } else if( val > high ) {
2159 FPRINTF(mif,"Variable ");
2160 print_var_name(mif,sys,var);
2161 FPRINTF(mif," was initialized above its upper bound.\n");
2162 FPRINTF(mif,"It will be moved to its upper bound.\n");
2163 var_set_value(var, high);
2164 }
2165 }
2166
2167 static void move_to_next_block(sys)
2168 slv0_system_t sys;
2169 /**
2170 *** Moves on to the next block, updating all of the solver information.
2171 *** To move to the first block, set sys->s.block.current_block to -1 before
2172 *** calling. If already at the last block, then sys->s.block.current_block
2173 *** will equal the number of blocks and the system will be declared
2174 *** converged. Otherwise, the residuals for the new block will be computed
2175 *** and sys->s.calc_ok set according.
2176 **/
2177 {
2178 #ifdef KAA_DEBUG
2179 int32 klowcol, khighcol;
2180 int32 klowrow, khighrow;
2181 struct rel_relation *krel;
2182 struct var_variable *kvar;
2183 double ktime =0;
2184 int k, ktmp;
2185 #endif
2186
2187 /* sys->clock = tm_cpu_time(); */
2188
2189 if( sys->s.block.current_block >= 0 ) {
2190
2191 int32 row;
2192 int32 col;
2193 int32 ci;
2194 struct var_variable *var;
2195 struct rel_relation *rel;
2196
2197 /* Record cost accounting info here. */
2198 ci=sys->s.block.current_block;
2199 sys->s.cost[ci].size=sys->s.block.current_size;
2200 sys->s.cost[ci].iterations=sys->s.block.iteration;
2201 sys->s.cost[ci].funcs=sys->s.block.funcs;
2202 sys->s.cost[ci].jacs=sys->s.block.jacs;
2203 sys->s.cost[ci].functime=sys->s.block.functime;
2204 sys->s.cost[ci].jactime=sys->s.block.jactime;
2205 sys->s.cost[ci].time=sys->s.block.cpu_elapsed;
2206 sys->s.cost[ci].resid=sys->s.block.residual;
2207
2208 #undef KAA_DEBUG
2209 #ifdef KAA_DEBUG
2210 if (sys->s.block.current_size>1) {
2211 FPRINTF(stderr,"Block %d took %d iterations\n",
2212 sys->s.block.current_block,sys->s.block.iteration);
2213 }
2214 #endif /* KAA_DEBUG */
2215
2216 /* De-initialize previous block */
2217 if (sys->p.output.less_important && (sys->s.block.current_size >1 ||
2218 sys->iarray[SP0_LIFDS])) {
2219 FPRINTF(LIF(sys),"Block %d converged.\n",
2220 sys->s.block.current_block);
2221 }
2222 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
2223 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2224 var_set_in_block(var,FALSE);
2225 }
2226 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2227 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2228 rel_set_in_block(rel,FALSE);
2229 }
2230 sys->s.block.previous_total_size += sys->s.block.current_size;
2231 }
2232
2233 sys->s.block.current_block += 1;
2234 if( sys->s.block.current_block < sys->s.block.number_of ) {
2235 int32 row;
2236 int32 col;
2237 boolean ok;
2238
2239 /* Initialize next block */
2240 if( OPTIMIZING(sys) )
2241 mtx_region(&(sys->J.reg), 0, sys->rank-1, 0, sys->vused-1 );
2242 else
2243 mtx_block(sys->J.mtx,sys->s.block.current_block,&(sys->J.reg));
2244
2245 row = sys->J.reg.row.high - sys->J.reg.row.low + 1;
2246 col = sys->J.reg.col.high - sys->J.reg.col.low + 1;
2247 sys->s.block.current_size = MAX(row,col);
2248
2249 sys->s.block.iteration = 0;
2250 sys->s.block.funcs = 0;
2251 sys->s.block.jacs = 0;
2252 sys->s.block.cpu_elapsed = 0.0;
2253 sys->s.block.functime = 0.0;
2254 sys->s.block.jactime = 0.0;
2255
2256 /*************************************************************************/
2257 for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
2258 struct var_variable *var;
2259 var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
2260 var_set_in_block(var,TRUE);
2261 insure_bounds(sys,var);
2262 }
2263 for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
2264 struct rel_relation *rel;
2265 rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
2266 rel_set_in_block(rel,TRUE);
2267 }
2268 /* REORDER */
2269 if( sys->s.block.current_size > 1L )
2270 linsol_reorder(sys->J.sys,&(sys->J.reg));
2271
2272 if(sys->p.output.less_important && (sys->iarray[SP0_LIFDS] ||
2273 sys->s.block.current_size > 1)) {
2274 debug_delimiter(LIF(sys));
2275 debug_delimiter(LIF(sys));
2276 }
2277 if(sys->p.output.less_important && sys->iarray[SP0_LIFDS]) {
2278 FPRINTF(LIF(sys),"\n%-40s ---> %d in [%d..%d]\n",
2279 "Current block number", sys->s.block.current_block,
2280 0, sys->s.block.number_of-1);
2281 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Current block size",
2282 sys->s.block.current_size);
2283 }
2284 sys->s.calc_ok = TRUE;
2285
2286 if( !(ok = calc_objective(sys)) )
2287 FPRINTF(MIF(sys),"Objective calculation errors detected.\n");
2288 if(sys->p.output.less_important && sys->obj) {
2289 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Objective", sys->objective);
2290 }
2291 sys->s.calc_ok = sys->s.calc_ok && ok;
2292
2293 sys->residuals.accurate = FALSE;
2294 if( !(ok = calc_residuals(sys)) )
2295 FPRINTF(MIF(sys),"Residual calculation errors detected.\n");
2296 if(sys->p.output.less_important && (sys->s.block.current_size >1 ||
2297 sys->iarray[SP0_LIFDS])) {
2298 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
2299 sys->s.block.residual);
2300 }
2301 sys->s.calc_ok = sys->s.calc_ok && ok;
2302
2303 /* Must be updated as soon as required */
2304 sys->J.accurate = FALSE;
2305 sys->update.weights = 0;
2306 sys->update.nominals = 0;
2307 sys->ZBZ.accurate = FALSE;
2308 sys->variables.accurate = FALSE;
2309 sys->gradient.accurate = FALSE;
2310 sys->multipliers.accurate = FALSE;
2311 sys->stationary.accurate = FALSE;
2312 sys->newton.accurate = FALSE;
2313 sys->Bnewton.accurate = FALSE;
2314 sys->nullspace.accurate = FALSE;
2315 sys->gamma.accurate = FALSE;
2316 sys->Jgamma.accurate = FALSE;
2317 sys->varstep1.accurate = FALSE;
2318 sys->Bvarstep1.accurate = FALSE;
2319 sys->varstep2.accurate = FALSE;
2320 sys->Bvarstep2.accurate = FALSE;
2321 sys->mulstep1.accurate = FALSE;
2322 sys->mulstep2.accurate = FALSE;
2323 sys->varstep.accurate = FALSE;
2324 sys->mulstep.accurate = FALSE;
2325
2326 if( !OPTIMIZING(sys) ) {
2327 sys->ZBZ.accurate = TRUE;
2328 sys->gradient.accurate = TRUE;
2329 sys->multipliers.accurate = TRUE;
2330 sys->stationary.accurate = TRUE;
2331 sys->Bnewton.accurate = TRUE;
2332 sys->nullspace.accurate = TRUE;
2333 sys->Bvarstep1.accurate = TRUE;
2334 sys->Bvarstep2.accurate = TRUE;
2335 }
2336
2337 } else {
2338 boolean ok;
2339 /**
2340 *** Before we claim convergence, we must check if we left behind
2341 *** some unassigned relations. If and only if they happen to be
2342 *** satisfied at the current point, convergence has been obtained.
2343 ***
2344 *** Also insures that all included relations have valid residuals.
2345 *** Included inequalities will have correct residuals.
2346 *** Unsatisfied included inequalities cause inconsistency.
2347 ***
2348 *** This of course ignores that fact an objective function might
2349 *** be present. Then, feasibility isn't enough, is it now.
2350 **/
2351 if( sys->s.struct_singular ) {
2352 /* black box w/singletons provoking bug here, maybe */
2353 sys->s.block.current_size = sys->rused - sys->rank;
2354 if(sys->p.output.less_important) {
2355 debug_delimiter(LIF(sys));
2356 FPRINTF(LIF(sys),"%-40s ---> %d\n", "Unassigned Relations",
2357 sys->s.block.current_size);
2358 }
2359 sys->J.reg.row.low = sys->J.reg.col.low = sys->rank;
2360 sys->J.reg.row.high = sys->J.reg.col.high = sys->rused - 1;
2361 sys->residuals.accurate = FALSE;
2362 if( !(ok=calc_residuals(sys)) )
2363 FPRINTF(MIF(sys),"Residual calculation errors detected.\n");
2364 if(sys->p.output.less_important) {
2365 FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
2366 sys->s.block.residual);
2367 }
2368 if( block_feasible(sys) ) {
2369 if(sys->p.output.less_important) {
2370 FPRINTF(LIF(sys),"\nUnassigned relations ok. You lucked out.\n");
2371 }
2372 sys->s.converged = TRUE;
2373 } else {
2374 if(sys->p.output.less_important) {
2375 FPRINTF(LIF(sys),"\nProblem inconsistent: %s.\n",
2376 "Unassigned relations not satisfied");
2377 }
2378 sys->s.inconsistent = TRUE;
2379 }
2380 if(sys->p.output.less_important) {
2381 debug_delimiter(LIF(sys));
2382 }
2383 } else {
2384 sys->s.converged = TRUE;
2385 }
2386 /* nearly done checking. Must verify included inequalities if
2387 we think equalities are ok. */
2388 if (sys->s.converged) {
2389 sys->s.inconsistent=(!calc_inequalities(sys));
2390 }
2391 }
2392 }
2393
2394 static void find_next_unconverged_block(sys)
2395 slv0_system_t sys;
2396 /**
2397 *** Moves to next unconverged block, assuming that the current block has
2398 *** converged (or is -1, to start).
2399 **/
2400 {
2401 do {
2402 move_to_next_block(sys);
2403 } while( !sys->s.converged && block_feasible(sys) && !OPTIMIZING(sys) );
2404 }
2405
2406
2407 /**
2408 *** Iteration begin/end routines
2409 *** ----------------------------
2410 *** iteration_begins(sys)
2411 *** iteration_ends(sys)
2412 **/
2413
2414 static void iteration_begins(sys)
2415 slv0_system_t sys;
2416 /**
2417 *** Prepares sys for entering an iteration, increasing the iteration counts
2418 *** and starting the clock.
2419 **/
2420 {
2421 sys->clock = tm_cpu_time();
2422 ++(sys->s.block.iteration);
2423 ++(sys->s.iteration);
2424 if(sys->p.output.less_important&& (sys->s.block.current_size >1 ||
2425 sys->iarray[SP0_LIFDS])) {
2426 FPRINTF(LIF(sys),"\n%-40s ---> %d\n",
2427 "Iteration", sys->s.block.iteration);
2428 FPRINTF(LIF(sys),"%-40s ---> %d\n",
2429 "Total iteration", sys->s.iteration);
2430 }
2431 }
2432
2433 static void iteration_ends(sys)
2434 slv0_system_t sys;
2435 /**
2436 *** Prepares sys for exiting an iteration, stopping the clock and recording
2437 *** the cpu time.
2438 **/
2439 {
2440 double cpu_elapsed; /* elapsed this iteration */
2441
2442 cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
2443 sys->s.block.cpu_elapsed += cpu_elapsed;
2444 sys->s.cpu_elapsed += cpu_elapsed;
2445 if(sys->p.output.less_important && (sys->s.block.current_size >1 ||
2446 sys->iarray[SP0_LIFDS])) {
2447 FPRINTF(LIF(sys),"%-40s ---> %g\n",
2448 "Elapsed time", sys->s.block.cpu_elapsed);
2449 FPRINTF(LIF(sys),"%-40s ---> %g\n",
2450 "Total elapsed time", sys->s.cpu_elapsed);
2451 }
2452 }
2453
2454 static void update_status(sys)
2455 slv0_system_t sys;
2456 /**
2457 *** Updates the solver status.
2458 **/
2459 {
2460 boolean unsuccessful;
2461
2462 if( !sys->s.converged ) {
2463 sys->s.time_limit_exceeded =
2464 abs(sys->s.block.cpu_elapsed >= sys->p.time_limit);
2465 sys->s.iteration_limit_exceeded =
2466 abs(sys->s.block.iteration >= sys->p.iteration_limit);
2467 }
2468
2469 unsuccessful = sys->s.diverged || sys->s.inconsistent ||
2470 sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
2471
2472 sys->s.ready_to_solve =abs( !unsuccessful && !sys->s.converged);
2473 sys->s.ok = abs(!unsuccessful &&
2474 sys->s.calc_ok && !sys->s.struct_singular);
2475 }
2476
2477 /**
2478 *** External routines
2479 *** -----------------
2480 *** See slv.h
2481 **/
2482
2483 static SlvClientToken slv0_create(slv_system_t server, int *statusindex)
2484 {
2485 slv0_system_t sys;
2486 FILE *fplif, *fpmif;
2487
2488 fpmif = stderr; /* will set up a proper file handle eventually */
2489 fplif = fpmif;
2490
2491 sys = (slv0_system_t)ascmalloc( sizeof(struct slv0_system_structure) );
2492 if (sys==NULL) {
2493 *statusindex = 1;
2494 return sys;
2495 }
2496 mem_zero_byte_cast(sys,0,sizeof(struct slv0_system_structure));
2497 sys->integrity = OK;
2498 sys->p.output.more_important = stdout;
2499 sys->p.output.less_important = NULL;
2500 sys->p.tolerance.pivot = 0.1; /* scaled */
2501 sys->p.tolerance.singular = 1e-12; /* scaled */
2502 sys->p.tolerance.feasible = 1e-8; /* unscaled */
2503 sys->p.tolerance.stationary = 1e-8; /* scaled */
2504 sys->p.tolerance.termination = 1e-12; /* scaled */
2505 sys->p.time_limit = 1500.0;
2506 sys->p.iteration_limit = 100;
2507 sys->p.partition = TRUE;
2508 sys->p.ignore_bounds = FALSE;
2509 sys->p.whose = slv0_solver_number;
2510 sys->p.rho = 1.0;
2511 sys->p.sp.iap=&(sys->iarray[0]); /* all defaults in iarray are 0 */
2512 sys->s.ok = TRUE;
2513 sys->s.calc_ok = TRUE;
2514 sys->s.costsize = 0;
2515 sys->s.cost = NULL; /*redundant, but sanity preserving */
2516
2517 return(sys);
2518 }
2519
2520 static void destroy_matrices(sys)
2521 slv0_system_t sys;
2522 {
2523 if( sys->J.sys ) {
2524 int count = linsol_number_of_rhs(sys->J.sys)-1;
2525 for( ; count >= 0; count-- )
2526 destroy_array(linsol_get_rhs(sys->J.sys,count));
2527 mtx_destroy(linsol_get_matrix(sys->J.sys));
2528 linsol_destroy(sys->J.sys);
2529 if( sys->J.relpivots ) set_destroy( sys->J.relpivots );
2530 if( sys->J.varpivots ) set_destroy( sys->J.varpivots );
2531 sys->J.sys = NULL;
2532 }
2533 if( sys->B ) {
2534 struct hessian_data *update;
2535 for( update=sys->B; update != NULL; ) {
2536 struct hessian_data *handle;
2537 handle = update;
2538 update = update->next;
2539 destroy_array(handle->y.vec);
2540 destroy_array(handle->Bs.vec);
2541 ascfree(handle);
2542 }
2543 sys->B = NULL;
2544 }
2545 if( sys->ZBZ.order > 0 ) {
2546 int i;
2547 for( i = 0; i < sys->ZBZ.order; i++ )
2548 destroy_array(sys->ZBZ.mtx[i]);
2549 destroy_array(sys->ZBZ.mtx);
2550 destroy_array(sys->ZBZ.ZBs);
2551 destroy_array(sys->ZBZ.Zy);
2552 sys->ZBZ.order = 0;
2553 }
2554 }
2555
2556 static void destroy_vectors(sys)
2557 slv0_system_t sys;
2558 {
2559 destroy_array(sys->nominals.vec);
2560 destroy_array(sys->weights.vec);
2561 destroy_array(sys->variables.vec);
2562 destroy_array(sys->residuals.vec);
2563 destroy_array(sys->gradient.vec);
2564 destroy_array(sys->multipliers.vec);
2565 destroy_array(sys->stationary.vec);
2566 destroy_array(sys->newton.vec);
2567 destroy_array(sys->Bnewton.vec);
2568 destroy_array(sys->nullspace.vec);
2569 destroy_array(sys->gamma.vec);
2570 destroy_array(sys->Jgamma.vec);
2571 destroy_array(sys->varstep1.vec);
2572 destroy_array(sys->Bvarstep1.vec);
2573 destroy_array(sys->varstep2.vec);
2574 destroy_array(sys->Bvarstep2.vec);
2575 destroy_array(sys->mulstep1.vec);
2576 destroy_array(sys->mulstep2.vec);
2577 destroy_array(sys->varstep.vec);
2578 destroy_array(sys->mulstep.vec);
2579 }
2580
2581 static void slv0_closefiles(sys)
2582 slv0_system_t sys;
2583 {
2584 return;
2585 }
2586
2587 static void slv0_set_var_list(slv0_system_t sys,struct var_variable **vlist)
2588 {
2589 static struct var_variable *empty_list[] = {NULL};
2590 check_system(sys);
2591 if( sys->vlist_user == NULL )
2592 if( sys->vlist != NULL && pl_length(sys->vlist) > 0 )
2593 ascfree( (POINTER)(sys->vlist) );
2594 sys->vlist_user = vlist;
2595 sys->vlist = (vlist==NULL ? empty_list : vlist);
2596 sys->s.ready_to_solve = FALSE;
2597 }
2598
2599 static struct var_variable **slv0_get_var_list(sys)
2600 slv0_system_t sys;
2601 {
2602 check_system(sys);
2603 return( sys->vlist_user );
2604 }
2605
2606 static void slv0_set_rel_list(slv0_system_t sys, struct rel_relation **rlist)
2607 {
2608 static struct rel_relation *empty_list[] = {NULL};
2609 check_system(sys);
2610 sys->rlist_user = rlist;
2611 sys->rlist = (rlist==NULL ? empty_list : rlist);
2612 sys->s.ready_to_solve = FALSE;
2613 }
2614
2615 static struct rel_relation **slv0_get_rel_list(sys)
2616 slv0_system_t sys;
2617 {
2618 check_system(sys);
2619 return( sys->rlist_user );
2620 }
2621
2622 static int slv0_count_vars(sys,vfilter)
2623 slv0_system_t sys;
2624 var_filter_t *vfilter;
2625 {
2626 struct var_variable **vp;
2627 int32 count = 0;
2628 check_system(sys);
2629 for( vp=sys->vlist; *vp != NULL; vp++ )
2630 if( var_apply_filter(*vp,vfilter) ) ++count;
2631 return( count );
2632 }
2633
2634 static int slv0_count_rels(slv0_system_t sys,rel_filter_t *rfilter)
2635 {
2636 struct rel_relation **rp;
2637 int32 count = 0;
2638 check_system(sys);
2639 for( rp=sys->rlist; *rp != NULL; rp++ )
2640 if( rel_apply_filter(*rp,rfilter) ) ++count;
2641 return( count );
2642 }
2643
2644 static int slv0_eligible_solver(slv_system_t server)
2645 {
2646 struct rel_relation **rp;
2647 for( rp=slv_get_solvers_rel_list(server); *rp != NULL ; ++rp ) {
2648 if( rel_less(*rp) || rel_greater(*rp) ) return(FALSE);
2649 }
2650 return(TRUE);
2651 }
2652
2653 static void slv0_get_parameters(slv_system_t server, SlvClientToken asys,
2654 slv_parameters_t *parameters)
2655 {
2656 slv0_system_t sys;
2657 sys = SLV0(asys);
2658 if (check_system(sys)) return;
2659 mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
2660 }
2661
2662 static void slv0_set_parameters(slv_system_t server, SlvClientToken asys,
2663 slv_parameters_t *parameters)
2664 {
2665 slv0_system_t sys;
2666 sys = SLV0(asys);
2667 if (check_system(sys)) return;
2668 mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
2669 }
2670
2671 static void slv0_get_status(slv_system_t server, SlvClientToken asys,
2672 slv_status_t *status)
2673 {
2674 slv0_system_t sys;
2675 sys = SLV0(asys);
2676 if (check_system(sys)) return;
2677 mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
2678 }
2679
2680 static linsolqr_system_t slv0_get_linsolqr_sys(slv_system_t server,
2681 SlvClientToken asys)
2682 {
2683 return NULL;
2684 }
2685
2686 static linsol_system_t slv0_get_linsol_sys(slv_system_t server,
2687 SlvClientToken asys)
2688 {
2689 slv0_system_t sys;
2690 sys = SLV0(asys);
2691 if (check_system(sys)) return NULL;
2692 return(sys->J.sys);
2693 }
2694
2695 static void sort_unassigned_rows(slv0_system_t sys) {
2696 /**
2697 *** Rearranges unassigned rows of the matrix as follows
2698 *** (note that unassigned rows may be from unassigned relations
2699 *** included or unincluded or empty rows that are unused):
2700 *** rows included && !assigned | rows of unincluded eqns | empty rows.
2701 *** Empty in this context means row unassociated with eqn, not whether
2702 *** there is incidence in that row.
2703 *** (mtx_output_assign no longer puts empties up into the 0 .. rtot-1
2704 *** region.)
2705 *** Joe correctly contends that this is a solver job and not an output
2706 *** assign postprocessing job mtx could do.
2707 *** (mtx_output assign does in fact leave unassigned rows _with_incidence
2708 *** up against the assigned region, but rows with no incidence
2709 *** are mixed up with unincluded rows.
2710 **/
2711 int32 i,j,k,rank,rinc,rtot,cap;
2712 mtx_matrix_t mtx;
2713 mtx=sys->J.mtx;
2714 rank=sys->rank;
2715 rinc=sys->rused;
2716 rtot=sys->rtot;
2717 cap=sys->cap;
2718 k=cap-1;
2719 j=rtot-1;
2720 /* now move all included but unassigned rows up against assigned rows */
2721 for (i=rank; i<rinc; i++) {
2722 if (!rel_included(sys->rlist[mtx_row_to_org(mtx,i)]) ||
2723 !rel_active(sys->rlist[mtx_row_to_org(mtx,i)]) ) {
2724 /* if unincluded found in range we will check for convergence ...*/
2725 while (!rel_included(sys->rlist[mtx_row_to_org(mtx,j)]) ||
2726 !rel_active(sys->rlist[mtx_row_to_org(mtx,j)]) ) j--;
2727 /* find last included row */
2728 if (j > i) mtx_swap_rows(mtx,i,j--);
2729 else FPRINTF(stderr,"Bug in slv0.c:sort_unassigned_rows\n");
2730 /* if this message is seen we would have swapped included into unincls
2731 an accounting error has taken place elsewhere in presolve*/
2732 }
2733 } /* all included but unassigned rows now against bottom of assigned region*/
2734 }
2735
2736 static void structural_analysis(slv_system_t server, slv0_system_t sys)
2737 /**
2738 *** Performs structural analysis on the system, setting the flags in
2739 *** status. The problem must be set up, the relation/variable list
2740 *** must be non-NULL, and their sizes must be in sys->s.v/r.total. The
2741 *** jacobian (linear) system must be created and have the correct order
2742 *** (stored in sys->cap). Everything else will be determined here.
2743 **/
2744 {
2745 struct var_variable **vp;
2746 var_filter_t vfilter;
2747 struct rel_relation **rp;
2748 rel_filter_t rfilter;
2749
2750 /**
2751 *** The server has marked incidenct flags already.
2752 **/
2753 /* count included equalities */
2754 rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
2755 rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
2756 sys->rused = slv_count_solvers_rels(server,&rfilter);
2757
2758 /* count free and incident vars */
2759 vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
2760 vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
2761 sys->vused = slv_count_solvers_vars(server,&vfilter);
2762
2763 /* Generate incidence matrix */
2764 mtx_clear(sys->J.mtx);
2765 for( rp = sys->rlist ; *rp != NULL ; ++rp )
2766 if( rel_apply_filter(*rp,&rfilter) ) {
2767 relman_get_incidence(*rp,&vfilter,sys->J.mtx);
2768 }
2769
2770 /* Symbolic analysis */
2771 sys->rtot = slv_get_num_solvers_rels(server);
2772 sys->vtot = slv_get_num_solvers_vars(server);
2773 mtx_output_assign(sys->J.mtx,sys->rtot,sys->vtot);
2774 /* symbolic rank is set */
2775 sys->rank = mtx_symbolic_rank(sys->J.mtx);
2776 sys->ZBZ.order = sys->obj ? (sys->vused - sys->rank) : 0;
2777 /* eventually we should let another client do this */
2778 if( sys->p.partition && !OPTIMIZING(sys) ) {
2779 mtx_partition(sys->J.mtx);
2780 }
2781 sort_unassigned_rows(sys);
2782
2783 /* Initialize Status */
2784 sys->s.over_defined = (sys->rused > sys->vused);
2785 sys->s.under_defined = (sys->rused < sys->vused);
2786 sys->s.struct_singular = (sys->rank < sys->rused);
2787 sys->s.block.number_of = mtx_number_of_blocks(sys->J.mtx);
2788 }
2789
2790 static void create_matrices(slv_system_t server, slv0_system_t sys)
2791 {
2792 sys->J.sys = linsol_create();
2793 sys->J.mtx = mtx_create();
2794 mtx_set_order(sys->J.mtx,sys->cap);
2795 linsol_set_matrix(sys->J.sys,sys->J.mtx);
2796 linsol_set_pivot_zero(sys->J.sys, sys->p.tolerance.singular);
2797 linsol_set_pivot_tolerance(sys->J.sys, sys->p.tolerance.pivot);
2798 /* rhs 0 for sys->multipliers */
2799 sys->J.rhs = create_array(sys->cap,real64);
2800 linsol_add_rhs(sys->J.sys,sys->J.rhs,TRUE);
2801 /* rhs 1 for sys->newton */
2802 sys->J.rhs = create_array(sys->cap,real64);
2803 linsol_add_rhs(sys->J.sys,sys->J.rhs,FALSE);
2804 /* rhs 2 for sys->mulstep2 */
2805 sys->J.rhs = create_array(sys->cap,real64);
2806 linsol_add_rhs(sys->J.sys,sys->J.rhs,TRUE);
2807 sys->J.relpivots = set_create(sys->cap);
2808 sys->J.varpivots = set_create(sys->cap);
2809
2810 structural_analysis(server,sys);
2811 sys->ZBZ.mtx = create_array(sys->ZBZ.order,real64 *);
2812 if( sys->ZBZ.mtx ) {
2813 int i;
2814 for( i=0; i<sys->ZBZ.order; i++ )
2815 sys->ZBZ.mtx[i] = create_array(sys->ZBZ.order,real64);
2816 }
2817 sys->ZBZ.ZBs = create_array(sys->ZBZ.order,real64);
2818 sys->ZBZ.Zy = create_array(sys->ZBZ.order,real64);
2819 }
2820
2821 static void create_vectors(sys)
2822 slv0_system_t sys;
2823 {
2824
2825 sys->nominals.vec = create_array(sys->cap,real64);
2826 sys->nominals.rng = &(sys->J.reg.col);
2827 sys->weights.vec = create_array(sys->cap,real64);
2828 sys->weights.rng = &(sys->J.reg.row);
2829 sys->variables.vec = create_array(sys->cap,real64);
2830 sys->variables.rng = &(sys->J.reg.col);
2831 sys->residuals.vec = create_array(sys->cap,real64);
2832 sys->residuals.rng = &(sys->J.reg.row);
2833 if( OPTIMIZING(sys) ) {
2834 sys->gradient.vec = create_array(sys->cap,real64);
2835 sys->gradient.rng = &(sys->J.reg.col);
2836 sys->multipliers.vec = create_array(sys->cap,real64);
2837 sys->multipliers.rng = &(sys->J.reg.row);
2838 sys->stationary.vec = create_array(sys->cap,real64);
2839