/[ascend]/branches/jacob/disused/solver/slv0.c
ViewVC logotype

Contents of /branches/jacob/disused/solver/slv0.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2933 - (show annotations) (download) (as text)
Mon Jun 1 20:38:42 2015 UTC (9 years, 1 month ago) by jacob
File MIME type: text/x-csrc
File size: 105821 byte(s)
Create a new branch from trunk r2932, with name 'jacob'
To be used by Jacob Shealy during GSOC2015 for development of ideal 
solution capabilities in FPROPS

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