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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (19 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 106022 byte(s)
Setting up web subdirectory in repository
1 aw0a 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