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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 136104 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.87 $
6     * Version control file: $RCSfile: slv3.c,v $
7     * Date last modified: $Date: 2000/01/25 02:27:32 $
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 <math.h>
34     #include <stdarg.h>
35     #include "utilities/ascConfig.h"
36     #include "utilities/ascSignal.h"
37     #include "utilities/ascMalloc.h"
38     #include "utilities/set.h"
39     #include "general/tm_time.h"
40     #include "utilities/mem.h"
41     /* #include "compiler/compiler.h" */
42     #include "utilities/ascPanic.h"
43     #include "general/list.h"
44     #include "compiler/fractions.h"
45     #include "compiler/dimen.h"
46     #include "compiler/functype.h"
47     #include "compiler/func.h"
48     #include "solver/mtx.h"
49     #include "solver/linsol.h"
50     #include "solver/linsolqr.h"
51     #include "solver/slv_types.h"
52     #include "solver/var.h"
53     #include "solver/rel.h"
54     #include "solver/discrete.h"
55     #include "solver/conditional.h"
56     #include "solver/logrel.h"
57     #include "solver/bnd.h"
58     #include "solver/calc.h"
59     #include "solver/relman.h"
60     #include "solver/slv_common.h"
61     #include "solver/slv_client.h"
62     #include "solver/slv3.h"
63     #include "solver/slv_stdcalls.h"
64    
65     #if !defined(STATIC_QRSLV) && !defined(DYNAMIC_QRSLV)
66     int slv3_register(SlvFunctionsT *f)
67     {
68     (void)f; /* stop gcc whine about unused parameter */
69    
70     FPRINTF(stderr,"QRSlv not compiled in this ASCEND IV.\n");
71     return 1;
72     }
73     #else /* either STATIC_QRSLV or DYNAMIC_QRSLV is defined */
74     #ifdef DYNAMIC_QRSLV
75     /* do dynamic loading stuff. yeah, right */
76     #else /* following is used if STATIC_QRSLV is defined */
77    
78    
79     #define KILL 0
80     /* code that needs to be deleted compiles only with kill = 1 */
81     #define CANOPTIMIZE FALSE
82     /* TRUE iff optimization code completed, meaning relman_diff fixed. */
83     #define DEBUG FALSE
84     /* makes lots of extra spew */
85    
86     #define SLV3(s) ((slv3_system_t)(s))
87     #define SERVER (sys->slv)
88     #define slv3_PA_SIZE 43 /* MUST INCREMENT WHEN ADDING PARAMETERS */
89     #define slv3_RA_SIZE 11
90    
91     /* do not delete (or extend) this array definition.
92     */
93     #define IEX(n) slv3_iaexpln[(n)]
94     #define slv3_IA_SIZE 16
95     static char *slv3_iaexpln[slv3_IA_SIZE] = {
96     "If lifds != 0 and showlessimportant is TRUE, show direct solve details",
97     "If savlin != 0, write out matrix data file at each iteration to SlvLinsol.dat",
98     "scale residuals by relation nominals for evaluating progress",
99     "Cutoff is the block size cutoff for MODEL-based reordering of partitions",
100     "Update jacobian every this many major iterations",
101     "Update row scalings every this many major iterations",
102     "Update column scalings every this many major iterations",
103     "Require misunderstood reduction somewhere in the stepping algorithm",
104     "Require residual >= some other number in the stepping algorithm",
105     "Check jacobian for poorly scaled columns and whine if found",
106     "Truncate whole step vector rather than componentwise at variable bound",
107     "Reorder option. 0 = MODEL based, 1 = MODEL based2, 2 = simple spk1",
108     "Use safe calculation routines",
109     "Update relation nominal scalings every this many major iterations",
110     "Max iterations for iterative scaling",
111     "scaleopt = 0: 2norm,= 1: relnom,= 2 2norm + iterative,= 3: relnom + iterative,= 4: iterative"
112     };
113    
114     /* change slv3_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
115     #define UPDATE_JACOBIAN_PTR (sys->parm_array[0])
116     #define UPDATE_JACOBIAN ((*(int *)UPDATE_JACOBIAN_PTR))
117     #define UPDATE_WEIGHTS_PTR (sys->parm_array[1])
118     #define UPDATE_WEIGHTS ((*(int *)UPDATE_WEIGHTS_PTR))
119     #define UPDATE_NOMINALS_PTR (sys->parm_array[2])
120     #define UPDATE_NOMINALS ((*(int *)UPDATE_NOMINALS_PTR))
121     #define REDUCE_PTR (sys->parm_array[3])
122     #define REDUCE ((*(int *)REDUCE_PTR))
123     #define EXACT_LINE_SEARCH_PTR (sys->parm_array[4])
124     #define EXACT_LINE_SEARCH ((*(int *)EXACT_LINE_SEARCH_PTR))
125     #define DUMPCNORM_PTR (sys->parm_array[5])
126     #define DUMPCNORM ((*(int *)DUMPCNORM_PTR))
127     #define TRUNCATE_PTR (sys->parm_array[6])
128     #define TRUNCATE ((*(int *)TRUNCATE_PTR))
129     #define SAFE_CALC_PTR (sys->parm_array[7])
130     #define SAFE_CALC ((*(int *)SAFE_CALC_PTR))
131     #define SCALEOPT_PTR (sys->parm_array[8])
132     #define SCALEOPT ((*(char **)SCALEOPT_PTR))
133     #define UPDATE_RELNOMS_PTR (sys->parm_array[9])
134     #define UPDATE_RELNOMS ((*(int *)UPDATE_RELNOMS_PTR))
135     #define ITSCALELIM_PTR (sys->parm_array[10])
136     #define ITSCALELIM ((*(int *)ITSCALELIM_PTR))
137     #define RELNOMSCALE_PTR (sys->parm_array[11])
138     #define RELNOMSCALE ((*(int *)RELNOMSCALE_PTR))
139     #define TOO_SMALL_PTR (sys->parm_array[12])
140     #define TOO_SMALL ((*(real64 *)TOO_SMALL_PTR))
141     #define CNLOW_PTR (sys->parm_array[13])
142     #define CNLOW ((*(real64 *)CNLOW_PTR))
143     #define CNHIGH_PTR (sys->parm_array[14])
144     #define CNHIGH ((*(real64 *)CNHIGH_PTR))
145     #define TOWARD_BOUNDS_PTR (sys->parm_array[15])
146     #define TOWARD_BOUNDS ((*(real64 *)TOWARD_BOUNDS_PTR))
147     #define POSITIVE_DEFINITE_PTR (sys->parm_array[16])
148     #define POSITIVE_DEFINITE ((*(real64 *)POSITIVE_DEFINITE_PTR))
149     #define DETZERO_PTR (sys->parm_array[17])
150     #define DETZERO ((*(real64 *)DETZERO_PTR))
151     #define STEPSIZEERR_MAX_PTR (sys->parm_array[18])
152     #define STEPSIZEERR_MAX ((*(real64 *)STEPSIZEERR_MAX_PTR))
153     #define PARMRNG_MIN_PTR (sys->parm_array[19])
154     #define PARMRNG_MIN ((*(real64 *)PARMRNG_MIN_PTR))
155     #define MIN_COEF_PTR (sys->parm_array[20])
156     #define MIN_COEF ((*(real64 *)MIN_COEF_PTR))
157     #define MAX_COEF_PTR (sys->parm_array[21])
158     #define MAX_COEF ((*(real64 *)MAX_COEF_PTR))
159     #define ITSCALETOL_PTR (sys->parm_array[22])
160     #define ITSCALETOL ((*(real64 *)ITSCALETOL_PTR))
161     #define IGNORE_BOUNDS_PTR (sys->parm_array[23])
162     #define IGNORE_BOUNDS ((*(int32 *)IGNORE_BOUNDS_PTR))
163     #define SHOW_MORE_IMPT_PTR (sys->parm_array[24])
164     #define SHOW_MORE_IMPT ((*(int32 *)SHOW_MORE_IMPT_PTR))
165     #define RHO_PTR (sys->parm_array[25])
166     #define RHO ((*(real64 *)RHO_PTR))
167     #define PARTITION_PTR (sys->parm_array[26])
168     #define PARTITION ((*(int32 *)PARTITION_PTR))
169     #define SHOW_LESS_IMPT_PTR (sys->parm_array[27])
170     #define SHOW_LESS_IMPT ((*(int32 *)SHOW_LESS_IMPT_PTR))
171     #define AUTO_RESOLVE_PTR (sys->parm_array[28])
172     #define AUTO_RESOLVE ((*(int32 *)AUTO_RESOLVE_PTR))
173     #define TIME_LIMIT_PTR (sys->parm_array[29])
174     #define TIME_LIMIT ((*(int32 *)TIME_LIMIT_PTR))
175     #define ITER_LIMIT_PTR (sys->parm_array[30])
176     #define ITER_LIMIT ((*(int32 *)ITER_LIMIT_PTR))
177     #define STAT_TOL_PTR (sys->parm_array[31])
178     #define STAT_TOL ((*(real64 *)STAT_TOL_PTR))
179     #define TERM_TOL_PTR (sys->parm_array[32])
180     #define TERM_TOL ((*(real64 *)TERM_TOL_PTR))
181     #define SING_TOL_PTR (sys->parm_array[33])
182     #define SING_TOL ((*(real64 *)SING_TOL_PTR))
183     #define PIVOT_TOL_PTR (sys->parm_array[34])
184     #define PIVOT_TOL ((*(real64 *)PIVOT_TOL_PTR))
185     #define FEAS_TOL_PTR (sys->parm_array[35])
186     #define FEAS_TOL ((*(real64 *)FEAS_TOL_PTR))
187     #define LIFDS_PTR (sys->parm_array[36])
188     #define LIFDS ((*(int32 *)LIFDS_PTR))
189     #define SAVLIN_PTR (sys->parm_array[37])
190     #define SAVLIN ((*(int32 *)SAVLIN_PTR))
191     #define REORDER_OPTION_PTR (sys->parm_array[38])
192     #define REORDER_OPTION ((*(char **)REORDER_OPTION_PTR))
193     #define CUTOFF_PTR (sys->parm_array[39])
194     #define CUTOFF ((*(int32 *)CUTOFF_PTR))
195     #define FACTOR_OPTION_PTR (sys->parm_array[40])
196     #define FACTOR_OPTION ((*(char **)FACTOR_OPTION_PTR))
197     #define CONVOPT_PTR (sys->parm_array[41])
198     #define CONVOPT ((*(char **)CONVOPT_PTR))
199     #define LINTIME_PTR (sys->parm_array[42])
200     #define LINTIME ((*(int *)LINTIME_PTR))
201     /* change slv3_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
202    
203    
204     #define REX(n) slv3_raexpln[(n)]
205     static
206     char *slv3_raexpln[slv3_RA_SIZE] = {
207     "Var nominal to use if user specifies 0.0",
208     "Smallest column norm we won't complain about if checking",
209     "Largest column norm we won't complain about if checking",
210     "If bound is in the way, we go this fraction toward it",
211     "Hessian fudge number when optimizing",
212     "Minimum 2x2 determinant of newton/gradient we consider non-parallel",
213     "Step size must be determined this precisely, or prngmin happy",
214     "Parameter range must be this narrow to exit inner loop if step size unhappy",
215     "'Largest' drop in maxstep allowed",
216     "'Smallest' drop in maxstep allowed",
217     "scale termination ratio for iterative method"};
218    
219     /*********************************************************************\
220     Subparameters implemented: (value/meaning)
221     LIFDS 0=>do not show full detail info for singletons
222     1=>do (this value ignored if detailed solve info not on.
223     SAVLIN 0=>do not append linearizations arising in the newton
224     scheme to the file SlvLinsol.dat.
225     1=>do.
226     SCALEOPT
227     0=>Use variable nominals and row two-norms for scaling
228     the Jacobian and rhs.
229     Use variable nominals and relation nominals for
230     scaling the Jacobian and rhs.
231     2=>Prescale by option 0 and then apply Fourer's
232     iterative scaling routine.
233     3=>Prescale by option 1 and then apply Fourer's
234     iterative scaling routine.
235     4=>Scale using only Fourer's iterative routine.
236     RELNOMSCALE
237     0=>use Jacobian row scaling for scaling residuals
238     for purpose of detecting descent.
239     1=>use most recently recorded relation nominals
240     for scaling residuals for purpose of
241     detecting descent.
242     The residuals will also be scaled by the
243     relation nominals AT THE CURRENT POINT
244     for determining constraint satisfaction.
245     UPRELNOM
246     0-INF=> Set number of iterations to wait
247     before updating vector of relation nominals.
248     CUTOFF] MODEL tearing/reordering cutoff number.
249    
250     [*] Generally cryptic parameters left by Joe. Someone
251     should play with and document them. See the defaults.
252    
253     \*********************************************************************/
254     struct update_data {
255     int jacobian; /* Countdown on jacobian updating */
256     int weights; /* Countdown on weights updating */
257     int nominals; /* Countdown on nominals updating */
258     int relnoms; /* Countdown on relnom updating */
259     int iterative; /* Countdown on iterative scale update */
260     };
261    
262     /* varpivots, relpivots used only in optimizing, if we rewrite calc_pivots
263     without them. */
264     struct jacobian_data {
265     linsolqr_system_t sys; /* Linear system */
266     mtx_matrix_t mtx; /* Transpose gradient of residuals */
267     real64 *rhs; /* RHS from linear system */
268     unsigned *varpivots; /* Pivoted variables */
269     unsigned *relpivots; /* Pivoted relations */
270     unsigned *subregions; /* Set of subregion indeces */
271     dof_t *dofdata; /* dof data pointer from server */
272     mtx_region_t reg; /* Current block region */
273     int32 rank; /* Numerical rank of the jacobian */
274     enum factor_method fm; /* Linear factorization method */
275     boolean accurate; /* ? Recalculate matrix */
276     boolean singular; /* ? Can matrix be inverted */
277     boolean old_partition; /* old value of partition flag */
278     };
279    
280     struct hessian_data {
281     struct vector_data Bs; /* Product of B and s */
282     struct vector_data y; /* Difference in stationaries */
283     real64 ys; /* inner product of y and s */
284     real64 sBs; /* inner product of s and Bs */
285     struct hessian_data *next; /* previous iteration data */
286     };
287    
288     struct reduced_data {
289     real64 **mtx; /* Dense matrix */
290     real64 *ZBs; /* Reduced Bs */
291     real64 *Zy; /* Reduced y */
292     int32 order; /* Degrees of freedom */
293     boolean accurate; /* Ready to re-compute ? */
294     };
295    
296     struct slv3_system_structure {
297    
298     /**
299     *** Problem definition
300     **/
301     slv_system_t slv; /* slv_system_t back-link */
302     struct rel_relation *obj; /* Objective function: NULL = none */
303     struct var_variable **vlist; /* Variable list (NULL terminated) */
304     struct rel_relation **rlist; /* Relation list (NULL terminated) */
305    
306     /**
307     *** Solver information
308     **/
309     int integrity; /* ? Has the system been created */
310     int32 presolved; /* ? Has the system been presolved */
311     slv_parameters_t p; /* Parameters */
312     slv_status_t s; /* Status (as of iteration end) */
313     struct update_data update; /* Jacobian frequency counters */
314     int32 cap; /* Order of matrix/vectors */
315     int32 rank; /* Symbolic rank of problem */
316     int32 vused; /* Free and incident variables */
317     int32 vtot; /* length of varlist */
318     int32 rused; /* Included relations */
319     int32 rtot; /* length of rellist */
320     double clock; /* CPU time */
321     void *parm_array[slv3_PA_SIZE]; /* array of pointers to param values */
322     struct slv_parameter pa[slv3_PA_SIZE];/* &pa[0] => sys->p.parms */
323    
324     /**
325     *** Calculated data (scaled)
326     **/
327     struct jacobian_data J; /* linearized system */
328     struct hessian_data *B; /* Curvature information */
329     struct reduced_data ZBZ; /* Reduced hessian */
330    
331     struct vector_data nominals; /* Variable nominals */
332     struct vector_data weights; /* Relation weights */
333     struct vector_data relnoms; /* Relation nominals */
334     struct vector_data variables; /* Variable values */
335     struct vector_data residuals; /* Relation residuals */
336     struct vector_data gradient; /* Objective gradient */
337     struct vector_data multipliers; /* Relation multipliers */
338     struct vector_data stationary; /* Lagrange gradient */
339     struct vector_data gamma; /* Feasibility steepest descent */
340     struct vector_data Jgamma; /* Product of J and gamma */
341     struct vector_data newton; /* Dependent variables */
342     struct vector_data Bnewton; /* Product of B and newton */
343     struct vector_data nullspace; /* Independent variables */
344     struct vector_data varstep1; /* 1st order in variables */
345     struct vector_data Bvarstep1; /* Product of B and varstep1 */
346     struct vector_data varstep2; /* 2nd order in variables */
347     struct vector_data Bvarstep2; /* Product of B and varstep2 */
348     struct vector_data mulstep1; /* 1st order in multipliers */
349     struct vector_data mulstep2; /* 2nd order in multipliers */
350     struct vector_data varstep; /* Step in variables */
351     struct vector_data mulstep; /* Step in multipliers */
352    
353     real64 objective; /* Objective function evaluation */
354     real64 phi; /* Unconstrained minimizer */
355     real64 maxstep; /* Maximum step size allowed */
356     real64 progress; /* Steepest directional derivative */
357     };
358    
359    
360     /**
361     *** Integrity checks
362     *** ----------------
363     *** check_system(sys)
364     **/
365    
366     #define OK ((int)813029392)
367     #define DESTROYED ((int)103289182)
368     static int check_system(slv3_system_t sys)
369     /**
370     *** Checks sys for NULL and for integrity.
371     **/
372     {
373     if( sys == NULL ) {
374     FPRINTF(stderr,"ERROR: (slv3) check_system\n");
375     FPRINTF(stderr," NULL system handle.\n");
376     return 1;
377     }
378    
379     switch( sys->integrity ) {
380     case OK:
381     return 0;
382     case DESTROYED:
383     FPRINTF(stderr,"ERROR: (slv3) check_system\n");
384     FPRINTF(stderr," System was recently destroyed.\n");
385     return 1;
386     default:
387     FPRINTF(stderr,"ERROR: (slv3) check_system\n");
388     FPRINTF(stderr," System reused or never allocated.\n");
389     return 1;
390     }
391     }
392    
393     /**
394     *** General input/output routines
395     *** -----------------------------
396     *** print_var_name(out,sys,var)
397     *** print_rel_name(out,sys,rel)
398     **/
399    
400     #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
401     #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
402    
403     /**
404     *** Debug output routines
405     *** ---------------------
406     *** debug_delimiter(fp)
407     *** debug_out_vector(fp,vec)
408     *** debug_out_var_values(fp,sys)
409     *** debug_out_rel_residuals(fp,sys)
410     *** debug_out_jacobian(fp,sys)
411     *** debug_out_hessian(fp,sys)
412     *** debug_write_array(fp,real64 *,length)
413     *** debug_out_parameters(fp)
414     **/
415    
416     static void debug_delimiter( FILE *fp)
417     /**
418     *** Outputs a hyphenated line.
419     **/
420     {
421     int i;
422     for( i=0; i<60; i++ ) PUTC('-',fp);
423     PUTC('\n',fp);
424     }
425    
426     #if DEBUG
427     static void debug_out_vector( FILE *fp, slv3_system_t sys,
428     struct vector_data *vec)
429     /**
430     *** Outputs a vector.
431     **/
432     {
433     int32 ndx;
434     FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n",
435     calc_sqrt_D0(vec->norm2), vec->accurate?"TRUE":"FALSE",
436     vec->rng->low,vec->rng->high);
437     FPRINTF(fp,"Vector --> ");
438     for( ndx=vec->rng->low ; ndx<=vec->rng->high ; ++ndx )
439     FPRINTF(fp, "%g ", vec->vec[ndx]);
440     PUTC('\n',fp);
441     }
442    
443     static void debug_out_var_values( FILE *fp, slv3_system_t sys)
444     /**
445     *** Outputs all variable values in current block.
446     **/
447     {
448     int32 col;
449     struct var_variable *var;
450    
451     FPRINTF(fp,"Var values --> \n");
452     for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high ; col++ ) {
453     var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
454     print_var_name(fp,sys,var);
455     FPRINTF(fp, "\nI Lb Value Ub Scale Col INom\n");
456     FPRINTF(fp,"%d\t%.4g\t%.4g\t%.4g\t%.4g\t%d\t%.4g\n",
457     var_sindex(var),var_lower_bound(var),var_value(var),
458     var_upper_bound(var),var_nominal(var),
459     col,sys->nominals.vec[col]);
460     }
461     }
462    
463     static void debug_out_rel_residuals( FILE *fp, slv3_system_t sys)
464     /**
465     *** Outputs all relation residuals in current block.
466     **/
467     {
468     int32 row;
469    
470     FPRINTF(fp,"Rel residuals --> \n");
471     for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high ; row++ ) {
472     struct rel_relation *rel;
473     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
474     FPRINTF(fp," %g : ",rel_residual(rel));
475     print_rel_name(fp,sys,rel);
476     PUTC('\n',fp);
477     }
478     PUTC('\n',fp);
479     }
480    
481     static void debug_out_jacobian( FILE *fp, slv3_system_t sys)
482     /**
483     *** Outputs permutation and values of the nonzero elements in the
484     *** the jacobian matrix.
485     **/
486     {
487     mtx_coord_t nz;
488     real64 value;
489    
490     nz.row = sys->J.reg.row.low;
491     for( ; nz.row <= sys->J.reg.row.high; ++(nz.row) ) {
492     FPRINTF(fp," Row %d (rel %d)\n", nz.row,
493     mtx_row_to_org(sys->J.mtx,nz.row));
494     nz.col = mtx_FIRST;
495     while( value = mtx_next_in_row(sys->J.mtx,&nz,&(sys->J.reg.col)),
496     nz.col != mtx_LAST ) {
497     FPRINTF(fp," Col %d (var %d) has value %g\n", nz.col,
498     mtx_col_to_org(sys->J.mtx,nz.col), value);
499     }
500     }
501     }
502    
503     static void debug_out_hessian( FILE *fp, slv3_system_t sys)
504     /**
505     *** Outputs permutation and values of the nonzero elements in the
506     *** reduced hessian matrix.
507     **/
508     {
509     mtx_coord_t nz;
510    
511     for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
512     nz.col = nz.row + sys->J.reg.col.high + 1 - sys->ZBZ.order;
513     FPRINTF(fp," ZBZ[%d (var %d)] = ",
514     nz.row, mtx_col_to_org(sys->J.mtx,nz.col));
515     for( nz.col = 0; nz.col < sys->ZBZ.order; nz.col++ ) {
516     FPRINTF(fp,"%10g ",sys->ZBZ.mtx[nz.row][nz.col]);
517     }
518     PUTC('\n',fp);
519     }
520     }
521    
522     #endif
523    
524     static void debug_write_array(FILE *fp,real64 *vec, int32 length)
525     {
526     int32 i;
527     for (i=0; i< length;i++)
528     FPRINTF(fp,"%.20g\n",vec[i]);
529     }
530    
531     static char savlinfilename[]="SlvLinsol.dat. \0";
532     static char savlinfilebase[]="SlvLinsol.dat.\0";
533     static int savlinnum=0;
534     /** The number to postfix to savlinfilebase. increases with file accesses. **/
535    
536     /**
537     *** Array/vector operations
538     *** ----------------------------
539     *** destroy_array(p)
540     *** create_array(len,type)
541     ***
542     *** zero_vector(vec)
543     *** copy_vector(vec1,vec2)
544     *** prod = inner_product(vec1,vec2)
545     *** norm2 = square_norm(vec)
546     *** matrix_product(mtx,vec,prod,scale,transpose)
547     **/
548    
549     #define destroy_array(p) \
550     if( (p) != NULL ) ascfree((p))
551     #define create_array(len,type) \
552     ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
553     #define create_zero_array(len,type) \
554     ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
555    
556     #define zero_vector(v) slv_zero_vector(v)
557     #define copy_vector(v,t) slv_copy_vector((v),(t))
558     #define inner_product(v,u) slv_inner_product((v),(u))
559     #define square_norm(v) slv_square_norm(v)
560     #define matrix_product(m,v,p,s,t) slv_matrix_product((m),(v),(p),(s),(t))
561    
562     /**
563     *** Calculation routines
564     *** --------------------
565     *** ok = calc_objective(sys)
566     *** ok = calc_boundaries(sys)
567     *** ok = calc_residuals(sys)
568     *** ok = calc_J(sys)
569     *** calc_nominals(sys)
570     *** calc_weights(sys)
571     *** scale_J(sys)
572     *** scale_variables(sys)
573     *** scale_residuals(sys)
574     *** ok = calc_gradient(sys)
575     *** calc_B(sys)
576     *** calc_ZBZ(sys)
577     *** calc_pivots(sys)
578     *** calc_rhs(sys)
579     *** calc_multipliers(sys)
580     *** calc_stationary(sys)
581     *** calc_newton(sys)
582     *** calc_Bnewton(sys)
583     *** calc_nullspace(sys)
584     *** calc_gamma(sys)
585     *** calc_Jgamma(sys)
586     *** calc_varstep1(sys)
587     *** calc_Bvarstep1(sys)
588     *** calc_varstep2(sys)
589     *** calc_Bvarstep2(sys)
590     *** calc_mulstep1(sys)
591     *** calc_mulstep2(sys)
592     *** calc_varstep(sys)
593     *** calc_mulstep(sys)
594     *** calc_phi(sys)
595     **/
596    
597    
598     #define OPTIMIZING(sys) ((sys)->ZBZ.order > 0)
599    
600     static boolean calc_objective( slv3_system_t sys)
601     /**
602     *** Evaluate the objective function.
603     **/
604     {
605     calc_ok = TRUE;
606     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
607     sys->objective = (sys->obj ? relman_eval(sys->obj,&calc_ok,SAFE_CALC) : 0.0);
608     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
609     return calc_ok;
610     }
611    
612     static boolean calc_objectives( slv3_system_t sys)
613     /**
614     *** Evaluate all objectives.
615     **/
616     {
617     int32 len,i;
618     static rel_filter_t rfilter;
619     struct rel_relation **rlist=NULL;
620     rfilter.matchbits = (REL_INCLUDED);
621     rfilter.matchvalue =(REL_INCLUDED);
622     rlist = slv_get_solvers_obj_list(SERVER);
623     len = slv_get_num_solvers_objs(SERVER);
624     calc_ok = TRUE;
625     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
626     for (i = 0; i < len; i++) {
627     if (rel_apply_filter(rlist[i],&rfilter)) {
628     relman_eval(rlist[i],&calc_ok,SAFE_CALC);
629     #if DEBUG
630     if (calc_ok == FALSE) {
631     FPRINTF(stderr,"error in calc_objectives\n");
632     calc_ok = TRUE;
633     }
634     #endif
635     }
636     }
637     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
638     return calc_ok;
639     }
640    
641    
642     static boolean calc_inequalities( slv3_system_t sys)
643     /**
644     *** Calculates all of the residuals of included inequalities.
645     *** Returns true iff (calculations preceded without error and
646     *** all inequalities are satisfied.)
647     **/
648     {
649     struct rel_relation **rp;
650     boolean satisfied=TRUE;
651     static rel_filter_t rfilter;
652     rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
653     rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
654    
655     calc_ok = TRUE;
656     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
657     for (rp=sys->rlist;*rp != NULL; rp++) {
658     if (rel_apply_filter(*rp,&rfilter)) {
659     relman_eval(*rp,&calc_ok,SAFE_CALC);
660     satisfied= satisfied &&
661     relman_calc_satisfied(*rp,FEAS_TOL);
662     }
663     }
664     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
665     return (calc_ok && satisfied);
666     }
667    
668     static boolean calc_residuals( slv3_system_t sys)
669     /**
670     *** Calculates all of the residuals in the current block and computes
671     *** the residual norm for block status. Returns true iff calculations
672     *** preceded without error.
673     **/
674     {
675     int32 row;
676     struct rel_relation *rel;
677     double time0;
678    
679     if( sys->residuals.accurate ) return TRUE;
680    
681     calc_ok = TRUE;
682     row = sys->residuals.rng->low;
683     time0=tm_cpu_time();
684     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
685     for( ; row <= sys->residuals.rng->high; row++ ) {
686     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
687     #if DEBUG
688     if (!rel) {
689     int r;
690     r=mtx_row_to_org(sys->J.mtx,row);
691     FPRINTF(stderr,"NULL relation found !!\n");
692     FPRINTF(stderr,"at row %d rel %d in calc_residuals\n",(int)row,r);
693     FFLUSH(stderr);
694     }
695     #endif
696     sys->residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
697    
698     if (strcmp(CONVOPT,"ABSOLUTE") == 0) {
699     relman_calc_satisfied(rel,FEAS_TOL);
700     } else if (strcmp(CONVOPT,"RELNOM_SCALE") == 0) {
701     relman_calc_satisfied_scaled(rel,FEAS_TOL);
702     }
703     }
704     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
705     sys->s.block.functime += (tm_cpu_time() -time0);
706     sys->s.block.funcs++;
707     square_norm( &(sys->residuals) );
708     sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
709     return(calc_ok);
710     }
711    
712    
713     static boolean calc_J( slv3_system_t sys)
714     /**
715     *** Calculates the current block of the jacobian.
716     *** It is initially unscaled.
717     **/
718     {
719     int32 row;
720     var_filter_t vfilter;
721     double time0;
722     real64 resid;
723    
724     if( sys->J.accurate )
725     return TRUE;
726    
727     calc_ok = TRUE;
728     vfilter.matchbits = (VAR_INBLOCK | VAR_ACTIVE);
729     vfilter.matchvalue = (VAR_INBLOCK | VAR_ACTIVE);
730     time0=tm_cpu_time();
731     mtx_clear_region(sys->J.mtx,&(sys->J.reg));
732     for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
733     struct rel_relation *rel;
734     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
735     relman_diffs(rel,&vfilter,sys->J.mtx,&resid,SAFE_CALC);
736     }
737     sys->s.block.jactime += (tm_cpu_time() - time0);
738     sys->s.block.jacs++;
739    
740     if( --(sys->update.nominals) <= 0 ) sys->nominals.accurate = FALSE;
741     if( --(sys->update.weights) <= 0 ) sys->weights.accurate = FALSE;
742    
743     linsolqr_matrix_was_changed(sys->J.sys);
744     return(calc_ok);
745     }
746    
747    
748     static void calc_nominals( slv3_system_t sys)
749     /**
750     *** Retrieves the nominal values of all of the block variables,
751     *** insuring that they are all strictly positive.
752     **/
753     {
754     int32 col;
755     FILE *fp = MIF(sys);
756     if( sys->nominals.accurate ) return;
757     fp = MIF(sys);
758     col = sys->nominals.rng->low;
759     if(strcmp(SCALEOPT,"NONE") == 0 ||
760     strcmp(SCALEOPT,"ITERATIVE") == 0){
761     for( ; col <= sys->nominals.rng->high; col++ ) {
762     sys->nominals.vec[col] = 1;
763     }
764     } else {
765     for( ; col <= sys->nominals.rng->high; col++ ) {
766     struct var_variable *var;
767     real64 n;
768     var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
769     n = var_nominal(var);
770     if( n <= 0.0 ) {
771     if( n == 0.0 ) {
772     n = TOO_SMALL;
773     FPRINTF(fp,"ERROR: (slv3) calc_nominals\n");
774     FPRINTF(fp," Variable ");
775     print_var_name(fp,sys,var);
776     FPRINTF(fp," \nhas nominal value of zero.\n");
777     FPRINTF(fp," Resetting to %g.\n",n);
778     var_set_nominal(var,n);
779     } else {
780     n = -n;
781     FPRINTF(fp,"ERROR: (slv3) calc_nominals\n");
782     FPRINTF(fp," Variable ");
783     print_var_name(fp,sys,var);
784     FPRINTF(fp," \nhas negative nominal value.\n");
785     FPRINTF(fp," Resetting to %g.\n",n);
786     var_set_nominal(var,n);
787     }
788     }
789     #if DEBUG
790     FPRINTF(fp,"Column %d is");
791     print_var_name(fp,sys,var);
792     FPRINTF(fp,"\nScaling of column %d is %g\n",col,n);
793     #endif
794     sys->nominals.vec[col] = n;
795     }
796     }
797     square_norm( &(sys->nominals) );
798     sys->update.nominals = UPDATE_NOMINALS;
799     sys->nominals.accurate = TRUE;
800     }
801    
802     static void calc_weights( slv3_system_t sys)
803     /**
804     *** Calculates the weights of all of the block relations
805     *** to scale the rows of the Jacobian.
806     **/
807     {
808     mtx_coord_t nz;
809     real64 sum;
810    
811     if( sys->weights.accurate )
812     return;
813    
814     nz.row = sys->weights.rng->low;
815     if(strcmp(SCALEOPT,"NONE") == 0 ||
816     strcmp(SCALEOPT,"ITERATIVE") == 0) {
817     for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
818     sys->weights.vec[nz.row] = 1;
819     }
820     } else if (strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
821     strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0) {
822     for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
823     sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
824     sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
825     }
826     } else if (strcmp(SCALEOPT,"RELNOM") == 0 ||
827     strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0) {
828     for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
829     sys->weights.vec[nz.row] =
830     1.0/rel_nominal(sys->rlist[mtx_row_to_org(sys->J.mtx,nz.row)]);
831     }
832     }
833     square_norm( &(sys->weights) );
834     sys->update.weights = UPDATE_WEIGHTS;
835     sys->residuals.accurate = FALSE;
836     sys->weights.accurate = TRUE;
837     }
838    
839     static void scale_J( slv3_system_t sys)
840     /**
841     *** Scales the jacobian.
842     **/
843     {
844     int32 row;
845     int32 col;
846    
847     if( sys->J.accurate ) return;
848    
849     calc_nominals(sys);
850     for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ )
851     mtx_mult_col(sys->J.mtx,col,sys->nominals.vec[col],&(sys->J.reg.row));
852    
853     calc_weights(sys);
854     for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ )
855     mtx_mult_row(sys->J.mtx,row,sys->weights.vec[row],&(sys->J.reg.col));
856     }
857    
858     static void jacobian_scaled(slv3_system_t sys)
859     {
860     int32 col;
861     if (DUMPCNORM) {
862     for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
863     real64 cnorm;
864     cnorm =
865     calc_sqrt_D0(mtx_sum_sqrs_in_col(sys->J.mtx,col,&(sys->J.reg.row)));
866     if (cnorm >CNHIGH || cnorm <CNLOW) {
867     FPRINTF(stderr,"[col %d org %d] %g\n", col,
868     mtx_col_to_org(sys->J.mtx,col), cnorm);
869     }
870     }
871     }
872    
873     sys->update.jacobian = UPDATE_JACOBIAN;
874     sys->J.accurate = TRUE;
875     sys->J.singular = FALSE; /* yet to be determined */
876     #if DEBUG
877     FPRINTF(LIF(sys),"\nJacobian: \n");
878     debug_out_jacobian(LIF(sys),sys);
879     #endif
880     }
881    
882     static void scale_variables( slv3_system_t sys)
883     {
884     int32 col;
885    
886     if( sys->variables.accurate ) return;
887    
888     col = sys->variables.rng->low;
889     for( ; col <= sys->variables.rng->high; col++ ) {
890     struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
891     sys->variables.vec[col] = var_value(var)/sys->nominals.vec[col];
892     }
893     square_norm( &(sys->variables) );
894     sys->variables.accurate = TRUE;
895     #if DEBUG
896     FPRINTF(LIF(sys),"Variables: ");
897     debug_out_vector(LIF(sys),sys,&(sys->variables));
898     #endif
899     }
900    
901     static void scale_residuals( slv3_system_t sys)
902     /**
903     *** Scales the previously calculated residuals.
904     **/
905     {
906     int32 row;
907    
908     if( sys->residuals.accurate ) return;
909    
910     row = sys->residuals.rng->low;
911     for( ; row <= sys->residuals.rng->high; row++ ) {
912     struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
913     sys->residuals.vec[row] = rel_residual(rel)*sys->weights.vec[row];
914     }
915     square_norm( &(sys->residuals) );
916     sys->residuals.accurate = TRUE;
917     #if DEBUG
918     FPRINTF(LIF(sys),"Residuals: ");
919     debug_out_vector(LIF(sys),sys,&(sys->residuals));
920     #endif
921     }
922    
923     static void calc_relnoms(slv3_system_t sys)
924     /**
925     *** Calculates relnoms for all relations in sys
926     *** using variable nominals.
927     **/
928     {
929     int32 row, col;
930     struct var_variable *var;
931     struct rel_relation *rel;
932     real64 *var_list;
933     var_list = create_array(sys->cap,real64);
934     col = 0;
935     var = sys->vlist[col];
936     /* store current variable values and
937     set variable value to nominal value */
938     while(var != NULL){
939     var_list[col] = var_value(var);
940     var_set_value(var, var_nominal(var));
941     col++;
942     var = sys->vlist[col];
943     }
944     row = 0;
945     rel = sys->rlist[row];
946     /* calculate relation nominal */
947     while(rel != NULL){
948     relman_scale(rel);
949     row++;
950     rel = sys->rlist[row];
951     }
952     col = 0;
953     var = sys->vlist[col];
954     /* restore variable values */
955     while(var != NULL){
956     var_set_value(var, var_list[col]);
957     col++;
958     var = sys->vlist[col];
959     }
960     destroy_array(var_list);
961     }
962    
963    
964     static real64 col_max_ratio(mtx_matrix_t *mtx,
965     mtx_region_t *reg)
966     /**
967     *** Returns the maximum ratio of magnitudes of any two nonzero
968     *** elements in the same column of mtx. Only considers elements
969     *** in region reg.
970     **/
971     {
972     real64 ratio;
973     real64 max_ratio;
974     real64 num, denom, dummy;
975     mtx_coord_t coord;
976     max_ratio = 0;
977     for(coord.col = reg->col.low;
978     coord.col <= reg->col.high; coord.col++) {
979     ratio = 0;
980     num = mtx_col_max(*mtx,&(coord),&(reg->row),&(dummy));
981     denom = mtx_col_min(*mtx,&(coord),&(reg->row),&(dummy),1e-7);
982     if(denom >0){
983     ratio = num/denom;
984     }
985     if(ratio > 10000000){
986     /* FPRINTF(stderr,"HELPME\n");*/
987     }
988     if(ratio > max_ratio){
989     max_ratio = ratio;
990     }
991     }
992     if(max_ratio == 0){
993     max_ratio = 1;
994     }
995     return max_ratio;
996     }
997     static real64 row_max_ratio(mtx_matrix_t *mtx,
998     mtx_region_t *reg)
999     /**
1000     *** Returns the maximum ratio of magnitudes of any two nonzero
1001     *** elements in the same row of mtx. Only considers elements
1002     *** in region reg.
1003     **/
1004     {
1005     real64 ratio;
1006     real64 max_ratio;
1007     real64 num, denom, dummy;
1008     mtx_coord_t coord;
1009     max_ratio = 0;
1010     for(coord.row = reg->row.low;
1011     coord.row <= reg->row.high; coord.row++) {
1012     ratio = 0;
1013     num = mtx_row_max(*mtx,&(coord),&(reg->col),&(dummy));
1014     denom = mtx_row_min(*mtx,&(coord),&(reg->col),&(dummy),1e-7);
1015     if(denom >0){
1016     ratio = num/denom;
1017     }
1018     if(ratio > 10000000){
1019     /* FPRINTF(stderr,"HELPME\n");*/
1020     }
1021     if(ratio > max_ratio){
1022     max_ratio = ratio;
1023     }
1024     }
1025     if(max_ratio == 0){
1026     max_ratio = 1;
1027     }
1028     return max_ratio;
1029     }
1030    
1031     static real64 calc_fourer_scale(mtx_matrix_t mtx,
1032     mtx_region_t reg,
1033     int32 loc,
1034     int32 option)
1035     /**
1036     *** Calculates scaling factor suggested by Fourer.
1037     *** For option = 0, returns scaling factor for
1038     *** row number loc.
1039     *** For option = 1, returns scaling factor for
1040     *** col number loc.
1041     **/
1042     {
1043     mtx_coord_t coord;
1044     real64 min, max, dummy;
1045     real64 scale;
1046     if(option == 0){
1047     if((loc < reg.row.low) || (loc > reg.row.high)){
1048     return 1;
1049     }
1050     coord.row = loc;
1051     min = mtx_row_min(mtx,&(coord),&(reg.col),&(dummy),1e-7);
1052     max = mtx_row_max(mtx,&(coord),&(reg.col),&(dummy));
1053     scale = min*max;
1054     if(scale > 0){
1055     scale = sqrt(scale);
1056     } else {
1057     scale = 1;
1058     }
1059     return scale;
1060     } else {
1061     if(loc < reg.col.low || loc > reg.col.high){
1062     return 1;
1063     }
1064     coord.col = loc;
1065     min = mtx_col_min(mtx,&(coord),&(reg.row),&(dummy),1e-7);
1066     max = mtx_col_max(mtx,&(coord),&(reg.row),&(dummy));
1067     scale = min*max;
1068     if(scale > 0){
1069     scale = sqrt(scale);
1070     } else {
1071     scale = 1;
1072     }
1073     return scale;
1074     }
1075     }
1076    
1077     static void scale_J_iterative(slv3_system_t sys)
1078     /**
1079     *** This funcion is an implementation of the scaling
1080     *** routine by Fourer on p304 of Mathematical Programing
1081     *** vol 23, (1982).
1082     *** This function will scale the Jacobian and store the scaling
1083     *** factors in sys->nominals and sys->weights.
1084     *** If the Jacobian has been previously scaled
1085     *** by another method (during this iteration) then these vectors
1086     *** should contain the scale factors used in that scaling.
1087     **/
1088     {
1089     real64 rho_col_old, rho_col_new;
1090     real64 rho_row_old, rho_row_new;
1091     int32 k;
1092     int32 done;
1093     int32 row, col;
1094     real64 *colvec = sys->nominals.vec;
1095     real64 *rowvec = sys->weights.vec;
1096     real64 rowscale, colscale;
1097     rho_col_old = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1098     rho_row_old = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1099     k = 0;
1100     done = 0;
1101     while(done == 0){
1102     k++;
1103     for(row = sys->J.reg.row.low;
1104     row <= sys->J.reg.row.high; row++){
1105     rowscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,row,0);
1106     mtx_mult_row(sys->J.mtx,row,rowscale,&(sys->J.reg.col));
1107     rowvec[row] *= rowscale;
1108     }
1109     for(col = sys->J.reg.col.low;
1110     col <= sys->J.reg.col.high; col++){
1111     colscale = 1/calc_fourer_scale(sys->J.mtx,sys->J.reg,col,1);
1112     mtx_mult_col(sys->J.mtx,col,colscale,&(sys->J.reg.row));
1113     colvec[col] *= colscale;
1114     }
1115     rho_col_new = col_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1116     rho_row_new = row_max_ratio(&(sys->J.mtx),&(sys->J.reg));
1117     if((rho_col_new >= ITSCALETOL*rho_col_old &&
1118     rho_row_new >= ITSCALETOL*rho_row_old)
1119     || k >= ITSCALELIM){
1120     done = 1;
1121     /* FPRINTF(stderr,"%d ITERATIVE SCALING ITERATIONS\n",k);*/
1122     } else {
1123     rho_row_old = rho_row_new;
1124     rho_col_old = rho_col_new;
1125     }
1126     }
1127     square_norm( &(sys->nominals) );
1128     sys->update.nominals = UPDATE_NOMINALS;
1129     sys->nominals.accurate = TRUE;
1130    
1131     square_norm( &(sys->weights) );
1132     sys->update.weights = UPDATE_WEIGHTS;
1133     sys->residuals.accurate = FALSE;
1134     sys->weights.accurate = TRUE;
1135     }
1136    
1137     static void scale_system( slv3_system_t sys )
1138     /**
1139     *** Scale system dependent on interface parameters
1140     **/
1141     {
1142     if(strcmp(SCALEOPT,"NONE") == 0){
1143     if(sys->J.accurate == FALSE){
1144     calc_nominals(sys);
1145     calc_weights(sys);
1146     jacobian_scaled(sys);
1147     }
1148     scale_variables(sys);
1149     scale_residuals(sys);
1150     return;
1151     }
1152     if(strcmp(SCALEOPT,"ROW_2NORM") == 0 ||
1153     strcmp(SCALEOPT,"RELNOM") == 0){
1154     if(sys->J.accurate == FALSE){
1155     scale_J(sys);
1156     jacobian_scaled(sys);
1157     }
1158     scale_variables(sys);
1159     scale_residuals(sys);
1160     return;
1161     }
1162     if(strcmp(SCALEOPT,"2NORM+ITERATIVE") == 0 ||
1163     strcmp(SCALEOPT,"RELNOM+ITERATIVE") == 0){
1164     if(sys->J.accurate == FALSE){
1165     --sys->update.iterative;
1166     if(sys->update.iterative <= 0) {
1167     scale_J(sys);
1168     scale_J_iterative(sys);
1169     sys->update.iterative =
1170     UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1171     } else {
1172     sys->weights.accurate = TRUE;
1173     sys->nominals.accurate = TRUE;
1174     scale_J(sys); /* will use current scaling vectors */
1175     }
1176     jacobian_scaled(sys);
1177     }
1178     scale_variables(sys);
1179     scale_residuals(sys);
1180     return;
1181     }
1182     if(strcmp(SCALEOPT,"ITERATIVE") == 0){
1183     if(sys->J.accurate == FALSE){
1184     --sys->update.iterative;
1185     if(sys->update.iterative <= 0) {
1186     calc_nominals(sys);
1187     calc_weights(sys);
1188     scale_J_iterative(sys);
1189     sys->update.iterative =
1190     UPDATE_WEIGHTS < UPDATE_NOMINALS ? UPDATE_WEIGHTS : UPDATE_NOMINALS;
1191     } else {
1192     sys->weights.accurate = TRUE;
1193     sys->nominals.accurate = TRUE;
1194     scale_J(sys); /* will use current scaling vectors */
1195     }
1196     jacobian_scaled(sys);
1197     }
1198     scale_variables(sys);
1199     scale_residuals(sys);
1200     }
1201     return;
1202     }
1203    
1204     static boolean calc_gradient(slv3_system_t sys)
1205     /**
1206     *** Calculate scaled gradient of the objective function.
1207     **/
1208     {
1209     /*
1210     *
1211     * This entire function needs to be reimplemented with relman_diffs.
1212     *
1213     */
1214     if( sys->gradient.accurate ) return TRUE;
1215    
1216     calc_ok = TRUE;
1217     if ( !OPTIMIZING(sys) ) {
1218     zero_vector(&(sys->gradient));
1219     sys->gradient.norm2 = 0.0;
1220     } else {
1221     Asc_Panic(2, "calc_gradient", "Not implemented");
1222     #if CANOPTIMIZE
1223     real64 pd;
1224     const struct var_variable **vp;
1225     var_filter_t vfilter;
1226     vfilter.matchbits = (VAR_INBLOCK | VAR_SVAR | VAR_ACTIVE);
1227     vfilter.matchvalue = (VAR_INBLOCK | VAR_SVAR | VAR_ACTIVE);
1228     zero_vector(&(sys->gradient));
1229     /* the next line will core dump anyway since vp not null-terminated*/
1230     for( vp = rel_incidence_list(sys->obj) ; *vp != NULL ; ++vp ) {
1231     int32 col;
1232     col = mtx_org_to_col(sys->J.mtx,var_sindex(*vp));
1233     if( var_apply_filter(*vp,&vfilter) ) {
1234     /* the next line will core dump anyway since _diff not implemented */
1235     relman_diff(sys->obj,*vp,&pd,SAFE_CALC); /* barf */
1236     sys->gradient.vec[col] = sys->nominals.vec[col]*pd;
1237     }
1238     }
1239     #endif
1240     square_norm( &(sys->gradient) );
1241     }
1242     sys->gradient.accurate = TRUE;
1243     #if DEBUG
1244     FPRINTF(LIF(sys),"Gradient: ");
1245     debug_out_vector(LIF(sys),sys,&(sys->gradient));
1246     #endif
1247     return calc_ok;
1248     }
1249    
1250     static void create_update( slv3_system_t sys)
1251     /**
1252     *** Create a new hessian_data structure for storing
1253     *** latest update information.
1254     **/
1255     {
1256     struct hessian_data *update;
1257    
1258     if( !OPTIMIZING(sys) )
1259     return;
1260    
1261     update = (struct hessian_data *)ascmalloc(sizeof(struct hessian_data));
1262     update->y.vec = create_array(sys->cap,real64);
1263     update->y.rng = &(sys->J.reg.col);
1264     update->y.accurate = FALSE;
1265     update->Bs.vec = create_array(sys->cap,real64);
1266     update->Bs.rng = &(sys->J.reg.col);
1267     update->Bs.accurate = FALSE;
1268     update->next = sys->B;
1269     sys->B = update;
1270     }
1271    
1272    
1273     static void calc_B( slv3_system_t sys)
1274     /**
1275     *** Computes a rank 2 BFGS update to the hessian matrix
1276     *** B which accumulates curvature.
1277     **/
1278     {
1279     if( sys->s.block.iteration > 1 ) {
1280     create_update(sys);
1281     } else {
1282     if( sys->B ) {
1283     struct hessian_data *update;
1284     for( update=sys->B; update != NULL; ) {
1285     struct hessian_data *handle;
1286     handle = update;
1287     update = update->next;
1288     destroy_array(handle->y.vec);
1289     destroy_array(handle->Bs.vec);
1290     ascfree(handle);
1291     }
1292     sys->B = NULL;
1293     }
1294     }
1295     if( sys->B ) {
1296     real64 theta;
1297     /**
1298     *** The y vector
1299     **/
1300     if( !sys->B->y.accurate ) {
1301     int32 col;
1302     matrix_product(sys->J.mtx, &(sys->multipliers),
1303     &(sys->B->y), 1.0, TRUE);
1304     col = sys->B->y.rng->low;
1305     for( ; col <= sys->B->y.rng->high; col++ ) {
1306     sys->B->y.vec[col] += sys->gradient.vec[col] -
1307     sys->stationary.vec[col];
1308     }
1309     square_norm( &(sys->B->y) );
1310     sys->B->y.accurate = TRUE;
1311     }
1312    
1313     /**
1314     *** The Bs vector
1315     **/
1316     if( !sys->B->Bs.accurate ) {
1317     struct hessian_data *update;
1318     copy_vector(&(sys->varstep),&(sys->B->Bs));
1319     for( update=sys->B->next; update != NULL; update = update->next ) {
1320     int32 col;
1321     real64 ys = inner_product( &(update->y),&(sys->varstep) );
1322     real64 sBs = inner_product( &(update->Bs),&(sys->varstep) );
1323     col = sys->B->Bs.rng->low;
1324     for( ; col<=sys->B->Bs.rng->high; col++) {
1325     sys->B->Bs.vec[col] += update->ys > 0.0 ?
1326     (update->y.vec[col])*ys/update->ys : 0.0;
1327     sys->B->Bs.vec[col] -= update->sBs > 0.0 ?
1328     (update->Bs.vec[col])*sBs/update->sBs : 0.0;
1329     }
1330     }
1331     square_norm( &(sys->B->Bs) );
1332     sys->B->Bs.accurate = TRUE;
1333     }
1334    
1335     sys->B->ys = inner_product( &(sys->B->y),&(sys->varstep) );
1336     sys->B->sBs = inner_product( &(sys->B->Bs),&(sys->varstep) );
1337    
1338     if( sys->B->ys == 0.0 && sys->B->sBs == 0.0 ) {
1339     theta = 0.0;
1340     } else {
1341     theta = sys->B->ys < POSITIVE_DEFINITE*sys->B->sBs ?
1342     (1.0-POSITIVE_DEFINITE)*sys->B->sBs/(sys->B->sBs - sys->B->ys):1.0;
1343     }
1344     #if DEBUG
1345     FPRINTF(LIF(sys),"ys, sBs, PD, theta = %g, %g, %g, %g\n",
1346     sys->B->ys,
1347     sys->B->sBs,
1348     POSITIVE_DEFINITE,
1349     theta);
1350     #endif
1351     if( theta < 1.0 ) {
1352     int32 col;
1353     col = sys->B->y.rng->low;
1354     for( ; col <= sys->B->y.rng->high; col++ )
1355     sys->B->y.vec[col] = theta*sys->B->y.vec[col] +
1356     (1.0-theta)*sys->B->Bs.vec[col];
1357     square_norm( &(sys->B->y) );
1358     sys->B->ys = theta*sys->B->ys + (1.0-theta)*sys->B->sBs;
1359     }
1360     }
1361     }
1362    
1363    
1364     static int calc_pivots(slv3_system_t sys)
1365     /**
1366     *** Obtain the equations and variables which
1367     *** are able to be pivoted.
1368     *** return value is the row rank deficiency, which we hope is 0.
1369     **/
1370     {
1371     int row_rank_defect=0, oldtiming;
1372     linsolqr_system_t lsys = sys->J.sys;
1373     FILE *fp = LIF(sys);
1374    
1375     oldtiming = g_linsolqr_timing;
1376     g_linsolqr_timing =LINTIME;
1377     linsolqr_factor(lsys,sys->J.fm); /* factor */
1378     g_linsolqr_timing = oldtiming;
1379    
1380     if (OPTIMIZING(sys)) {
1381     /* need things for nullspace move. don't care about
1382     * dependency coefficiency in any circumstances at present.
1383     */
1384     linsolqr_calc_col_dependencies(lsys);
1385     set_null(sys->J.relpivots,sys->cap);
1386     set_null(sys->J.varpivots,sys->cap);
1387     linsolqr_get_pivot_sets(lsys,sys->J.relpivots,sys->J.varpivots);
1388     }
1389    
1390     sys->J.rank = linsolqr_rank(lsys);
1391     sys->J.singular = FALSE;
1392     row_rank_defect = sys->J.reg.row.high -
1393     sys->J.reg.row.low+1 - sys->J.rank;
1394     if( row_rank_defect > 0 ) {
1395     int32 row,krow;
1396     mtx_sparse_t *uprows=NULL;
1397     sys->J.singular = TRUE;
1398     uprows = linsolqr_unpivoted_rows(lsys);
1399     if (uprows !=NULL) {
1400     for( krow=0; krow < uprows->len ; krow++ ) {
1401     int32 org_row;
1402     struct rel_relation *rel;
1403    
1404     org_row = uprows->idata[krow];
1405     row = mtx_org_to_row(sys->J.mtx,org_row);
1406     rel = sys->rlist[org_row];
1407     FPRINTF(fp,"%-40s ---> ","Relation not pivoted");
1408     print_rel_name(fp,sys,rel);
1409     PUTC('\n',fp);
1410    
1411     /**
1412     *** assign zeros to the corresponding weights
1413     *** so that subsequent calls to "scale_residuals"
1414     *** will only measure the pivoted equations.
1415     **/
1416     sys->weights.vec[row] = 0.0;
1417     sys->residuals.vec[row] = 0.0;
1418     sys->residuals.accurate = FALSE;
1419     mtx_mult_row(sys->J.mtx,row,0.0,&(sys->J.reg.col));
1420     }
1421     mtx_destroy_sparse(uprows);
1422     }
1423     if( !sys->residuals.accurate ) {
1424     square_norm( &(sys->residuals) );
1425     sys->residuals.accurate = TRUE;
1426     sys->update.weights = 0; /* re-compute weights next iteration. */
1427     }
1428     }
1429     if( sys->J.rank < sys->J.reg.col.high-sys->J.reg.col.low+1 ) {
1430     int32 col,kcol;
1431     mtx_sparse_t *upcols=NULL;
1432     if (NOTNULL(upcols)) {
1433     for( kcol=0; upcols != NULL && kcol < upcols->len ; kcol++ ) {
1434     int32 org_col;
1435     struct var_variable *var;
1436    
1437     org_col = upcols->idata[kcol];
1438     col = mtx_org_to_col(sys->J.mtx,org_col);
1439     var = sys->vlist[org_col];
1440     FPRINTF(fp,"%-40s ---> ","Variable not pivoted");
1441     print_var_name(fp,sys,var);
1442     PUTC('\n',fp);
1443     /**
1444     *** If we're not optimizing (everything should be
1445     *** pivotable) or this was one of the dependent variables,
1446     *** consider this variable as if it were fixed.
1447     **/
1448     if( col <= sys->J.reg.col.high - sys->ZBZ.order ) {
1449     mtx_mult_col(sys->J.mtx,col,0.0,&(sys->J.reg.row));
1450     }
1451     }
1452     mtx_destroy_sparse(upcols);
1453     }
1454     }
1455     if (SHOW_LESS_IMPT) {
1456     FPRINTF(LIF(sys),"%-40s ---> %d (%s)\n","Jacobian rank", sys->J.rank,
1457     sys->J.singular ? "deficient":"full");
1458     FPRINTF(LIF(sys),"%-40s ---> %g\n","Smallest pivot",
1459     linsolqr_smallest_pivot(sys->J.sys));
1460     }
1461     return row_rank_defect;
1462     }
1463    
1464     static void calc_ZBZ(slv3_system_t sys)
1465     /**
1466     *** Updates the reduced hessian matrix.
1467     *** if !OPTIMIZING just sets zbz.accurate true and returns.
1468     **/
1469     {
1470     mtx_coord_t nz;
1471    
1472     if( sys->ZBZ.accurate ) return;
1473    
1474     for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1475     for( nz.col = 0; nz.col <= nz.row; nz.col++ ) {
1476     int32 col, depr, depc;
1477     col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
1478     depr = mtx_col_to_org(sys->J.mtx,col);
1479     col = nz.col+sys->J.reg.col.high+1-sys->ZBZ.order;
1480     depc = mtx_col_to_org(sys->J.mtx,col);
1481     sys->ZBZ.mtx[nz.row][nz.col] = (nz.row==nz.col ? 1.0 : 0.0);
1482     col = sys->J.reg.col.low;
1483     for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
1484     int32 ind = mtx_col_to_org(sys->J.mtx,col);
1485     if( set_is_member(sys->J.varpivots,ind) )
1486     sys->ZBZ.mtx[nz.row][nz.col] +=
1487     (-linsolqr_org_col_dependency(sys->J.sys,depr,ind))*
1488     (-linsolqr_org_col_dependency(sys->J.sys,depc,ind));
1489     }
1490     if( nz.row != nz.col ) {
1491     sys->ZBZ.mtx[nz.col][nz.row] =
1492     sys->ZBZ.mtx[nz.row][nz.col];
1493     }
1494     }
1495     }
1496     if( OPTIMIZING(sys) ) {
1497     struct hessian_data *update;
1498     for( update=sys->B; update != NULL; update = update->next ) {
1499     for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1500     int32 col, dep;
1501     col = nz.row + sys->J.reg.col.high + 1 - sys->ZBZ.order;
1502     dep = mtx_col_to_org(sys->J.mtx,col);
1503     sys->ZBZ.Zy[nz.row] = update->y.vec[col];
1504     sys->ZBZ.ZBs[nz.row] = update->Bs.vec[col];
1505     col = sys->J.reg.col.low;
1506     for( ; col <= sys->J.reg.col.high - sys->ZBZ.order; col++ ) {
1507     int32 ind = mtx_col_to_org(sys->J.mtx,col);
1508     if( set_is_member(sys->J.varpivots,ind) ) {
1509     sys->ZBZ.Zy[nz.row] += update->y.vec[col]*
1510     (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1511     sys->ZBZ.ZBs[nz.row] += update->Bs.vec[col]*
1512     (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1513     }
1514     }
1515     for( nz.col=0; nz.col <= nz.row; nz.col++ ) {
1516     sys->ZBZ.mtx[nz.row][nz.col] += update->ys > 0.0 ?
1517     sys->ZBZ.Zy[nz.row]*sys->ZBZ.Zy[nz.col]/update->ys : 0.0;
1518     sys->ZBZ.mtx[nz.row][nz.col] -= update->sBs > 0.0 ?
1519     sys->ZBZ.ZBs[nz.row]*sys->ZBZ.ZBs[nz.col]/update->sBs : 0.0;
1520     if( nz.row != nz.col ) {
1521     sys->ZBZ.mtx[nz.col][nz.row] =
1522     sys->ZBZ.mtx[nz.row][nz.col];
1523     }
1524     }
1525     }
1526     }
1527     }
1528     sys->ZBZ.accurate = TRUE;
1529     #if DEBUG
1530     FPRINTF(LIF(sys),"\nReduced Hessian: \n");
1531     debug_out_hessian(LIF(sys),sys);
1532     #endif
1533     }
1534    
1535    
1536     static void calc_rhs(slv3_system_t sys, struct vector_data *vec,
1537     real64 scalar, boolean transpose)
1538     /**
1539     *** Calculates just the jacobian RHS. This function should be used to
1540     *** supplement calculation of the jacobian. The vector vec must
1541     *** already be calculated and scaled so as to simply be added to the
1542     *** rhs. Caller is responsible for initially zeroing the rhs vector.
1543     **/
1544     {
1545     if( transpose ) { /* vec is indexed by col */
1546     int32 col;
1547     for( col=vec->rng->low; col<=vec->rng->high; col++ ) {
1548     sys->J.rhs[mtx_col_to_org(sys->J.mtx,col)] += scalar*vec->vec[col];
1549     }
1550     } else { /* vec is indexed by row */
1551     int32 row;
1552     for( row=vec->rng->low; row<=vec->rng->high; row++ ) {
1553     sys->J.rhs[mtx_row_to_org(sys->J.mtx,row)] += scalar*vec->vec[row];
1554     }
1555     }
1556     linsolqr_rhs_was_changed(sys->J.sys,sys->J.rhs);
1557     }
1558    
1559    
1560     static void calc_multipliers(slv3_system_t sys)
1561     /**
1562     *** Computes the lagrange multipliers for the equality constraints.
1563     **/
1564     {
1565     if( sys->multipliers.accurate )
1566     return;
1567    
1568     if ( !OPTIMIZING(sys) ) {
1569     zero_vector(&(sys->multipliers));
1570     sys->multipliers.norm2 = 0.0;
1571     } else {
1572     linsolqr_system_t lsys = sys->J.sys;
1573     int32 row;
1574     sys->J.rhs = linsolqr_get_rhs(lsys,0);
1575     mtx_zero_real64(sys->J.rhs,sys->cap);
1576     calc_rhs(sys, &(sys->gradient), -1.0, TRUE );
1577     linsolqr_solve(lsys,sys->J.rhs);
1578     row = sys->multipliers.rng->low;
1579     for( ; row <= sys->multipliers.rng->high; row++ ) {
1580     struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1581     sys->multipliers.vec[row] = linsolqr_var_value
1582     (lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row));
1583     rel_set_multiplier(rel,sys->multipliers.vec[row]*
1584     sys->weights.vec[row]);
1585    
1586     }
1587     if (SAVLIN) {
1588     FILE *ldat;
1589     int32 ov;
1590     sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1591     ldat=fopen(savlinfilename,"w");
1592     FPRINTF(ldat,
1593     "================= multipliersrhs (orgcoled) itn %d =====\n",
1594     sys->s.iteration);
1595     debug_write_array(ldat,sys->J.rhs,sys->cap);
1596     FPRINTF(ldat,
1597     "================= multipliers (orgrowed) ============\n");
1598     for(ov=0 ; ov < sys->cap; ov++ )
1599     FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1600     fclose(ldat);
1601     }
1602     square_norm( &(sys->multipliers) );
1603     }
1604     sys->multipliers.accurate = TRUE;
1605     #if DEBUG
1606     FPRINTF(LIF(sys),"Multipliers: ");
1607     debug_out_vector(LIF(sys),sys,&(sys->multipliers));
1608     #endif
1609     }
1610    
1611    
1612     static void calc_stationary( slv3_system_t sys)
1613     /**
1614     *** Computes the gradient of the lagrangian which
1615     *** should be zero at the optimum solution.
1616     **/
1617     {
1618     if( sys->stationary.accurate )
1619     return;
1620    
1621     if ( !OPTIMIZING(sys) ) {
1622     zero_vector(&(sys->stationary));
1623     sys->stationary.norm2 = 0.0;
1624     } else {
1625     int32 col;
1626     matrix_product(sys->J.mtx, &(sys->multipliers),
1627     &(sys->stationary), 1.0, TRUE);
1628     col = sys->stationary.rng->low;
1629     for( ; col <= sys->stationary.rng->high; col++ )
1630     sys->stationary.vec[col] += sys->gradient.vec[col];
1631     square_norm( &(sys->stationary) );
1632     }
1633     sys->stationary.accurate = TRUE;
1634     #if DEBUG
1635     FPRINTF(LIF(sys),"Stationary: ");
1636     debug_out_vector(LIF(sys),sys,&(sys->stationary));
1637     #endif
1638     }
1639    
1640    
1641     static void calc_gamma( slv3_system_t sys)
1642     /**
1643     *** Calculate the gamma vector.
1644     **/
1645     {
1646     if( sys->gamma.accurate )
1647     return;
1648    
1649     matrix_product(sys->J.mtx, &(sys->residuals),
1650     &(sys->gamma), -1.0, TRUE);
1651     square_norm( &(sys->gamma) );
1652     sys->gamma.accurate = TRUE;
1653     #if DEBUG
1654     FPRINTF(LIF(sys),"Gamma: ");
1655     debug_out_vector(LIF(sys),sys,&(sys->gamma));
1656     #endif
1657     }
1658    
1659     static void calc_Jgamma( slv3_system_t sys)
1660     /**
1661     *** Calculate the Jgamma vector.
1662     **/
1663     {
1664     if( sys->Jgamma.accurate )
1665     return;
1666    
1667     matrix_product(sys->J.mtx, &(sys->gamma),
1668     &(sys->Jgamma), 1.0, FALSE);
1669     square_norm( &(sys->Jgamma) );
1670     sys->Jgamma.accurate = TRUE;
1671     #if DEBUG
1672     FPRINTF(LIF(sys),"Jgamma: ");
1673     debug_out_vector(LIF(sys),sys,&(sys->Jgamma));
1674     #endif
1675     }
1676    
1677    
1678     static void calc_newton( slv3_system_t sys)
1679     /**
1680     *** Computes a step to solve the linearized equations.
1681     **/
1682     {
1683     linsolqr_system_t lsys = sys->J.sys;
1684     int32 col;
1685    
1686     if( sys->newton.accurate )
1687     return;
1688    
1689     sys->J.rhs = linsolqr_get_rhs(lsys,1);
1690     mtx_zero_real64(sys->J.rhs,sys->cap);
1691     calc_rhs(sys, &(sys->residuals), -1.0, FALSE);
1692     linsolqr_solve(lsys,sys->J.rhs);
1693     col = sys->newton.rng->low;
1694     for( ; col <= sys->newton.rng->high; col++ ) {
1695     sys->newton.vec[col] =
1696     linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1697     }
1698     if (SAVLIN) {
1699     FILE *ldat;
1700     int32 ov;
1701     sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1702     ldat=fopen(savlinfilename,"w");
1703     FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1704     sys->s.iteration);
1705     debug_write_array(ldat,sys->J.rhs,sys->cap);
1706     FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1707     for(ov=0 ; ov < sys->cap; ov++ )
1708     FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1709     fclose(ldat);
1710     }
1711     square_norm( &(sys->newton) );
1712     sys->newton.accurate = TRUE;
1713     #if DEBUG
1714     FPRINTF(LIF(sys),"Newton: ");
1715     debug_out_vector(LIF(sys),sys,&(sys->newton));
1716     #endif
1717     }
1718    
1719    
1720     static void calc_Bnewton( slv3_system_t sys)
1721     /**
1722     *** Computes an update to the product B and newton.
1723     **/
1724     {
1725     if( sys->Bnewton.accurate )
1726     return;
1727    
1728     if ( !OPTIMIZING(sys) ) {
1729     zero_vector(&(sys->Bnewton));
1730     sys->Bnewton.norm2 = 0.0;
1731     } else {
1732     struct hessian_data *update;
1733     copy_vector(&(sys->newton),&(sys->Bnewton));
1734     for( update=sys->B; update != NULL; update = update->next ) {
1735     int32 col;
1736     real64 Yn = inner_product( &(update->y),&(sys->newton) );
1737     real64 sBn = inner_product( &(update->Bs),&(sys->newton) );
1738     col = sys->Bnewton.rng->low;
1739     for( ; col <= sys->Bnewton.rng->high; col++ ) {
1740     sys->Bnewton.vec[col] += update->ys > 0.0 ?
1741     (update->y.vec[col])*Yn/update->ys : 0.0;
1742     sys->Bnewton.vec[col] -= update->sBs > 0.0 ?
1743     (update->Bs.vec[col])*sBn/update->sBs : 0.0;
1744     }
1745     }
1746     square_norm( &(sys->Bnewton) );
1747     }
1748     sys->Bnewton.accurate = TRUE;
1749     #if DEBUG
1750     FPRINTF(LIF(sys),"Bnewton: ");
1751     debug_out_vector(LIF(sys),sys,&(sys->Bnewton));
1752     #endif
1753     }
1754    
1755    
1756     static void calc_nullspace( slv3_system_t sys)
1757     /**
1758     *** Calculate the nullspace move if OPTIMIZING.
1759     **/
1760     {
1761     if( sys->nullspace.accurate )
1762     return;
1763    
1764     if( !OPTIMIZING(sys) ) {
1765     zero_vector(&(sys->nullspace));
1766     sys->nullspace.norm2 = 0.0;
1767     } else {
1768     mtx_coord_t nz;
1769     zero_vector(&(sys->nullspace));
1770     for( nz.row=0; nz.row < sys->ZBZ.order; nz.row++ ) {
1771     int32 col, dep, ndx;
1772     col = nz.row+sys->J.reg.col.high+1-sys->ZBZ.order;
1773     dep = mtx_col_to_org(sys->J.mtx,col);
1774     sys->nullspace.vec[col] = -sys->stationary.vec[col] -
1775     sys->Bnewton.vec[col];
1776     ndx = sys->J.reg.col.low;
1777     for( ; ndx <= sys->J.reg.col.high - sys->ZBZ.order; ndx++ ) {
1778     int32 ind = mtx_col_to_org(sys->J.mtx,ndx);
1779     if( set_is_member(sys->J.varpivots,ind) )
1780     sys->nullspace.vec[col] -=
1781     (sys->stationary.vec[ndx] + sys->Bnewton.vec[ndx])*
1782     (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1783     }
1784     }
1785     /**
1786     *** Must invert ZBZ first. It's symmetric so
1787     *** can find Cholesky factors. Essentially, find
1788     *** the "square root" of the matrix such that
1789     ***
1790     *** T
1791     *** L U = U U = ZBZ, where U is an upper triangular
1792     *** matrix.
1793     **/
1794     for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1795     for( nz.col = nz.row; nz.col < sys->ZBZ.order; nz.col++ ) {
1796     int32 col;
1797     for( col = 0; col < nz.row; col++ )
1798     sys->ZBZ.mtx[nz.row][nz.col] -=
1799     sys->ZBZ.mtx[nz.row][col]*
1800     sys->ZBZ.mtx[col][nz.col];
1801     if( nz.row == nz.col )
1802     sys->ZBZ.mtx[nz.row][nz.col] =
1803     calc_sqrt_D0(sys->ZBZ.mtx[nz.row][nz.col]);
1804     else {
1805     sys->ZBZ.mtx[nz.row][nz.col] /=
1806     sys->ZBZ.mtx[nz.row][nz.row];
1807     sys->ZBZ.mtx[nz.col][nz.row] =
1808     sys->ZBZ.mtx[nz.row][nz.col];
1809     }
1810     }
1811     }
1812     #if DEBUG
1813     FPRINTF(LIF(sys),"\nInverse Reduced Hessian: \n");
1814     debug_out_hessian(LIF(sys),sys);
1815     #endif
1816     /**
1817     *** forward substitute
1818     **/
1819     for( nz.row = 0; nz.row < sys->ZBZ.order; nz.row++ ) {
1820     int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1821     for( nz.col = 0; nz.col < nz.row; nz.col++ ) {
1822     sys->nullspace.vec[nz.row+offset] -=
1823     sys->nullspace.vec[nz.col+offset]*
1824     sys->ZBZ.mtx[nz.row][nz.col];
1825     }
1826     sys->nullspace.vec[nz.row+offset] /=
1827     sys->ZBZ.mtx[nz.row][nz.row];
1828     }
1829    
1830     /**
1831     *** backward substitute
1832     **/
1833     for( nz.row = sys->ZBZ.order-1; nz.row >= 0; nz.row-- ) {
1834     int32 offset = sys->J.reg.col.high+1-sys->ZBZ.order;
1835     for( nz.col = nz.row+1; nz.col < sys->ZBZ.order; nz.col++ ) {
1836     sys->nullspace.vec[nz.row+offset] -=
1837     sys->nullspace.vec[nz.col+offset]*
1838     sys->ZBZ.mtx[nz.row][nz.col];
1839     }
1840     sys->nullspace.vec[nz.row+offset] /=
1841     sys->ZBZ.mtx[nz.row][nz.row];
1842     }
1843     square_norm( &(sys->nullspace) );
1844     }
1845     sys->nullspace.accurate = TRUE;
1846     #if DEBUG
1847     FPRINTF(LIF(sys),"Nullspace: ");
1848     debug_out_vector(LIF(sys),sys,&(sys->nullspace));
1849     #endif
1850     }
1851    
1852     static void calc_varstep1( slv3_system_t sys)
1853     /**
1854     *** Calculate the 1st order descent direction for phi
1855     *** in the variables.
1856     **/
1857     {
1858     if( sys->varstep1.accurate )
1859     return;
1860    
1861     if( !OPTIMIZING(sys) ) {
1862     copy_vector(&(sys->gamma),&(sys->varstep1));
1863     sys->varstep1.norm2 = sys->gamma.norm2;
1864     } else {
1865     int32 col;
1866     col = sys->varstep1.rng->low;
1867     for( ; col <= sys->varstep1.rng->high; col++ )
1868     sys->varstep1.vec[col] = RHO*sys->gamma.vec[col] -
1869     sys->stationary.vec[col];
1870     square_norm( &(sys->varstep1) );
1871     }
1872     sys->varstep1.accurate = TRUE;
1873     #if DEBUG
1874     FPRINTF(LIF(sys),"Varstep1: ");
1875     debug_out_vector(LIF(sys),sys,&(sys->varstep1));
1876     #endif
1877     }
1878    
1879    
1880     static void calc_Bvarstep1( slv3_system_t sys)
1881     /**
1882     *** Computes an update to the product B and varstep1.
1883     **/
1884     {
1885     if( sys->Bvarstep1.accurate )
1886     return;
1887    
1888     if ( !OPTIMIZING(sys) ) {
1889     zero_vector(&(sys->Bvarstep1));
1890     sys->Bvarstep1.norm2 = 0.0;
1891     } else {
1892     struct hessian_data *update;
1893     copy_vector(&(sys->varstep1),&(sys->Bvarstep1));
1894     for( update=sys->B; update != NULL; update = update->next ) {
1895     int32 col;
1896     real64 yv = inner_product( &(update->y),&(sys->varstep1) );
1897     real64 sBv = inner_product( &(update->Bs),&(sys->varstep1) );
1898     col = sys->Bvarstep1.rng->low;
1899     for( ; col <= sys->Bvarstep1.rng->high; col++ ) {
1900     sys->Bvarstep1.vec[col] += update->ys > 0.0 ?
1901     (update->y.vec[col])*yv/update->ys : 0.0;
1902     sys->Bvarstep1.vec[col] -= update->sBs > 0.0 ?
1903     (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1904     }
1905     }
1906     square_norm( &(sys->Bvarstep1) );
1907     }
1908     sys->Bvarstep1.accurate = TRUE;
1909     #if DEBUG
1910     FPRINTF(LIF(sys),"Bvarstep1: ");
1911     debug_out_vector(LIF(sys),sys,&(sys->Bvarstep1));
1912     #endif
1913     }
1914    
1915    
1916     static void calc_varstep2( slv3_system_t sys)
1917     /**
1918     *** Calculate the 2nd order descent direction for phi
1919     *** in the variables.
1920     **/
1921     {
1922     if( sys->varstep2.accurate )
1923     return;
1924    
1925     if( !OPTIMIZING(sys) ) {
1926     copy_vector(&(sys->newton),&(sys->varstep2));
1927     sys->varstep2.norm2 = sys->newton.norm2;
1928     } else {
1929     int32 col;
1930     col = sys->varstep2.rng->low;
1931     for( ; col <= sys->varstep2.rng->high - sys->ZBZ.order ; ++col ) {
1932     int32 dep;
1933     int32 ind = mtx_col_to_org(sys->J.mtx,col);
1934     sys->varstep2.vec[col] = sys->newton.vec[col];
1935     if( set_is_member(sys->J.varpivots,ind) ) {
1936     dep = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1937     for( ; dep <= sys->varstep2.rng->high; dep++ )
1938     sys->varstep2.vec[col] += sys->nullspace.vec[dep]*
1939     (-linsolqr_org_col_dependency(sys->J.sys,dep,ind));
1940     }
1941     }
1942     col = sys->varstep2.rng->high + 1 - sys->ZBZ.order;
1943     for( ; col <= sys->varstep2.rng->high; ++col )
1944     sys->varstep2.vec[col] = sys->nullspace.vec[col] +
1945     sys->newton.vec[col];
1946     square_norm( &(sys->varstep2) );
1947     }
1948     sys->varstep2.accurate = TRUE;
1949     #if DEBUG
1950     FPRINTF(LIF(sys),"Varstep2: ");
1951     debug_out_vector(LIF(sys),sys,&(sys->varstep2));
1952     #endif
1953     }
1954    
1955    
1956     static void calc_Bvarstep2( slv3_system_t sys)
1957     /**
1958     *** Computes an update to the product B and varstep2.
1959     **/
1960     {
1961     if( sys->Bvarstep2.accurate )
1962     return;
1963    
1964     if ( !OPTIMIZING(sys) ) {
1965     zero_vector(&(sys->Bvarstep2));
1966     sys->Bvarstep2.norm2 = 0.0;
1967     } else {
1968     struct hessian_data *update;
1969     copy_vector(&(sys->varstep2),&(sys->Bvarstep2));
1970     for( update=sys->B; update != NULL; update = update->next ) {
1971     int32 col;
1972     real64 yv = inner_product( &(update->y),&(sys->varstep2) );
1973     real64 sBv = inner_product( &(update->Bs),&(sys->varstep2) );
1974     col = sys->Bvarstep2.rng->low;
1975     for( ; col <= sys->Bvarstep2.rng->high; col++ ) {
1976     sys->Bvarstep2.vec[col] += update->ys > 0.0 ?
1977     (update->y.vec[col])*yv/update->ys : 0.0;
1978     sys->Bvarstep2.vec[col] -= update->sBs > 0.0 ?
1979     (update->Bs.vec[col])*sBv/update->sBs : 0.0;
1980     }
1981     }
1982     square_norm( &(sys->Bvarstep2) );
1983     }
1984     sys->Bvarstep2.accurate = TRUE;
1985     #if DEBUG
1986     FPRINTF(LIF(sys),"Bvarstep2: ");
1987     debug_out_vector(LIF(sys),sys,&(sys->Bvarstep2));
1988     #endif
1989     }
1990    
1991    
1992     static void calc_mulstep1( slv3_system_t sys)
1993     /**
1994     *** Calculate the negative gradient direction of phi in the
1995     *** multipliers.
1996     **/
1997     {
1998     if( sys->mulstep1.accurate )
1999     return;
2000    
2001     if( !OPTIMIZING(sys) ) {
2002     zero_vector(&(sys->mulstep1));
2003     sys->mulstep1.norm2 = 0.0;
2004     } else {
2005     int32 row;
2006     row = sys->mulstep1.rng->low;
2007     for( ; row <= sys->mulstep1.rng->high; row++ )
2008     sys->mulstep1.vec[row] = -sys->residuals.vec[row];
2009     square_norm( &(sys->mulstep1) );
2010     }
2011     sys->mulstep1.accurate = TRUE;
2012     #if DEBUG
2013     FPRINTF(LIF(sys),"Mulstep1: ");
2014     debug_out_vector(LIF(sys),sys,&(sys->mulstep1));
2015     #endif
2016     }
2017    
2018    
2019     static void calc_mulstep2( slv3_system_t sys)
2020     /**
2021     *** Calculate the mulstep2 direction of phi in the
2022     *** multipliers.
2023     **/
2024     {
2025     if( sys->mulstep2.accurate )
2026     return;
2027    
2028     if( !OPTIMIZING(sys) ) {
2029     zero_vector(&(sys->mulstep2));
2030     sys->mulstep2.norm2 = 0.0;
2031     } else {
2032     linsolqr_system_t lsys = sys->J.sys;
2033     int32 row;
2034     sys->J.rhs = linsolqr_get_rhs(lsys,2);
2035     mtx_zero_real64(sys->J.rhs,sys->cap);
2036     calc_rhs(sys, &(sys->Bvarstep2), -1.0, TRUE);
2037     calc_rhs(sys, &(sys->stationary), -1.0, TRUE);
2038     linsolqr_solve(lsys,sys->J.rhs);
2039     row = sys->mulstep2.rng->low;
2040     for( ; row <= sys->mulstep2.rng->high; row++ )
2041     sys->mulstep2.vec[row] = linsolqr_var_value
2042     (lsys,sys->J.rhs,mtx_row_to_org(sys->J.mtx,row));
2043     if (SAVLIN) {
2044     FILE *ldat;
2045     int32 ov;
2046     sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
2047     ldat=fopen(savlinfilename,"w");
2048     FPRINTF(ldat,
2049     "================= mulstep2rhs (orgcoled) itn %d =======\n",
2050     sys->s.iteration);
2051     debug_write_array(ldat,sys->J.rhs,sys->cap);
2052     FPRINTF(ldat,
2053     "================= mulstep2vars (orgrowed) ============\n");
2054     for(ov=0 ; ov < sys->cap; ov++ )
2055     FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
2056     fclose(ldat);
2057     }
2058     square_norm( &(sys->mulstep2) );
2059     }
2060     sys->mulstep2.accurate = TRUE;
2061     #if DEBUG
2062     FPRINTF(LIF(sys),"Mulstep2: ");
2063     debug_out_vector(LIF(sys),sys,&(sys->mulstep2));
2064     #endif
2065     }
2066    
2067    
2068     static void calc_phi( slv3_system_t sys)
2069     /**
2070     *** Computes the global minimizing function Phi.
2071     **/
2072     {
2073     if( !OPTIMIZING(sys) )
2074     sys->phi = 0.5*sys->residuals.norm2;
2075     else {
2076     sys->phi = sys->objective;
2077     sys->phi += inner_product( &(sys->multipliers),&(sys->residuals) );
2078     sys->phi += 0.5*RHO*sys->residuals.norm2;
2079     }
2080     }
2081    
2082     /**
2083     *** OK. Here's where we compute the actual step to be taken. It will
2084     *** be some linear combination of the 1st order and 2nd order steps.
2085     **/
2086    
2087     typedef real64 sym_2x2_t[3]; /* Stores symmetric 2x2 matrices */
2088    
2089     struct parms_t {
2090     real64 low,high,guess; /* Used to search for parameter */
2091     };
2092    
2093     struct calc_step_vars {
2094     sym_2x2_t coef1, coef2;
2095     real64 rhs[2]; /* RHS for 2x2 system */
2096     struct parms_t parms;
2097     real64 alpha1, alpha2;
2098     real64 error; /* Error between step norm and sys->maxstep */
2099     };
2100    
2101     static void calc_2x2_system(slv3_system_t sys, struct calc_step_vars *vars)
2102     /**
2103     *** Calculates 2x2 system (coef1,coef2,rhs).
2104     **/
2105     {
2106     vars->coef1[0] = (2.0*sys->phi/sys->newton.norm2)*
2107     calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
2108     vars->coef1[1] = 1.0;
2109     vars->coef1[2] = (sys->Jgamma.norm2/sys->gamma.norm2)*
2110     calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
2111    
2112     vars->coef2[0] = 1.0;
2113     vars->coef2[1] = 2.0*sys->phi/
2114     calc_sqrt_D0(sys->newton.norm2)/calc_sqrt_D0(sys->gamma.norm2);
2115     vars->coef2[2] = 1.0;
2116    
2117     vars->rhs[0] = 2.0*sys->phi/
2118     sys->maxstep/calc_sqrt_D0(sys->gamma.norm2);
2119     vars->rhs[1] = calc_sqrt_D0(sys->newton.norm2)/sys->maxstep;
2120     }
2121    
2122     static void coefs_from_parm( slv3_system_t sys, struct calc_step_vars *vars)
2123     /**
2124     *** Determines alpha1 and alpha2 from the parameter (guess).
2125     **/
2126     {
2127    
2128     sym_2x2_t coef; /* Actual coefficient matrix */
2129     real64 det; /* Determinant of coefficient matrix */
2130     int i;
2131    
2132     for( i=0 ; i<3 ; ++i ) coef[i] =
2133     vars->coef1[i] + vars->parms.guess * vars->coef2[i];
2134     det = coef[0]*coef[2] - coef[1]*coef[1];
2135     if( det < 0.0 )
2136     FPRINTF(MIF(sys),"%-40s ---> %g\n",
2137     " Unexpected negative determinant!",det);
2138     if( det <= DETZERO ) {
2139     /**
2140     *** varstep2 and varstep1 are essentially parallel:
2141     *** adjust length of either
2142     **/
2143     vars->alpha2 = 0.0;
2144     vars->alpha1 = 1.0;
2145     } else {
2146     vars->alpha2 = (vars->rhs[0]*coef[2] - vars->rhs[1]*coef[1])/det;
2147     vars->alpha1 = (vars->rhs[1]*coef[0] - vars->rhs[0]*coef[1])/det;
2148     }
2149     }
2150    
2151     static real64 step_norm2( slv3_system_t sys, struct calc_step_vars *vars)
2152     /**
2153     *** Computes step vector length based on 1st order and 2nd order
2154     *** vectors and their coefficients.
2155     **/
2156     {
2157     return sys->maxstep*sys->maxstep*
2158     (vars->alpha2 * vars->alpha2 +
2159     vars->alpha2 * vars->alpha1 * sys->phi/
2160     calc_sqrt_D0(sys->varstep2.norm2 + sys->mulstep2.norm2)/
2161     calc_sqrt_D0(sys->varstep1.norm2 + sys->mulstep1.norm2) +
2162     vars->alpha1 * vars->alpha1);
2163     }
2164    
2165     static void adjust_parms( slv3_system_t sys, struct calc_step_vars *vars)
2166     /**
2167     *** Re-guesses the parameters based on
2168     *** step size vs. target value.
2169     **/
2170     {
2171     vars->error = (calc_sqrt_D0(step_norm2(sys,vars))/sys->maxstep) - 1.0;
2172     if( vars->error > 0.0 ) {
2173     /* Increase parameter (to decrease step length) */
2174     vars->parms.low = vars->parms.guess;
2175     vars->parms.guess = (vars->parms.high>3.0*vars->parms.guess)
2176     ? 2.0*vars->parms.guess
2177     : 0.5*(vars->parms.low + vars->parms.high);
2178     } else {
2179     /* Decrease parameter (to increase step norm) */
2180     vars->parms.high = vars->parms.guess;
2181     vars->parms.guess = 0.5*(vars->parms.low + vars->parms.high);
2182     }
2183     }
2184    
2185     static void compute_step( slv3_system_t sys, struct calc_step_vars *vars)
2186     /**
2187     *** Computes the step based on the coefficients in vars.
2188     **/
2189     {
2190     int32 row,col;
2191     real64 tot1_norm2, tot2_norm2;
2192    
2193     tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
2194     tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
2195     if( !sys->varstep.accurate ) {
2196     for( col=sys->varstep.rng->low ; col<=sys->varstep.rng->high ; ++col )
2197     if( (vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
2198     sys->varstep.vec[col] = sys->maxstep *
2199     sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2);
2200     } else if( (vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
2201     sys->varstep.vec[col] = sys->maxstep *
2202     sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2);
2203     } else if( (vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
2204     sys->varstep.vec[col] = sys->maxstep*
2205     (
2206     vars->alpha2*sys->varstep2.vec[col]/calc_sqrt_D0(tot2_norm2) +
2207     vars->alpha1*sys->varstep1.vec[col]/calc_sqrt_D0(tot1_norm2)
2208     );
2209     }
2210     sys->varstep.accurate = TRUE;
2211     }
2212     if( !sys->mulstep.accurate ) {
2213     for( row=sys->mulstep.rng->low ; row<=sys->mulstep.rng->high ; ++row )
2214     if( (vars->alpha2 == 1.0) && (vars->alpha1 == 0.0) ) {
2215     sys->mulstep.vec[row] = sys->maxstep *
2216     sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2);
2217     } else if( (vars->alpha2 == 0.0) && (vars->alpha1 == 1.0) ) {
2218     sys->mulstep.vec[row] = sys->maxstep *
2219     sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2);
2220     } else if( (vars->alpha2 != 0.0) && (vars->alpha1 != 0.0) ) {
2221     sys->mulstep.vec[row] = sys->maxstep*
2222     (
2223     vars->alpha2*sys->mulstep2.vec[row]/calc_sqrt_D0(tot2_norm2) +
2224     vars->alpha1*sys->mulstep1.vec[row]/calc_sqrt_D0(tot1_norm2)
2225     );
2226     }
2227     sys->mulstep.accurate = TRUE;
2228     }
2229     #if DEBUG
2230     FPRINTF(LIF(sys),"Varstep: ");
2231     debug_out_vector(LIF(sys),sys,&(sys->varstep));
2232     FPRINTF(LIF(sys),"Mulstep: ");
2233     debug_out_vector(LIF(sys),sys,&(sys->mulstep));
2234     #endif
2235     }
2236    
2237    
2238     static void calc_step( slv3_system_t sys, int minor)
2239     /**
2240     *** Calculates step vector, based on sys->maxstep, and the varstep2/
2241     *** varstep1 and mulstep2/mulstep1 vectors. Nothing is assumed to be
2242     *** calculated, except the weights and the jacobian (scaled). Also,
2243     *** the step is not checked for legitimacy.
2244     *** NOTE: the step is scaled.
2245     **/
2246     {
2247    
2248     struct calc_step_vars vars;
2249     real64 tot1_norm2, tot2_norm2;
2250    
2251     if( sys->varstep.accurate && sys->mulstep.accurate )
2252     return;
2253     if (SHOW_LESS_IMPT) {
2254     FPRINTF(LIF(sys),"\n%-40s ---> %d\n", " Step trial",minor);
2255     }
2256    
2257     tot1_norm2 = sys->varstep1.norm2 + sys->mulstep1.norm2;
2258     tot2_norm2 = sys->varstep2.norm2 + sys->mulstep2.norm2;
2259     if( (tot1_norm2 == 0.0) && (tot2_norm2 == 0.0) ) {
2260     /* Take no step at all */
2261     vars.alpha1 = 0.0;
2262     vars.alpha2 = 0.0;
2263     sys->maxstep = 0.0;
2264     sys->varstep.norm2 = 0.0;
2265     sys->mulstep.norm2 = 0.0;
2266    
2267     } else if( (tot2_norm2 > 0.0) && OPTIMIZING(sys) ) {
2268     /* Stay in varstep2 direction */
2269     vars.alpha1 = 0.0;
2270     vars.alpha2 = 1.0;
2271     sys->maxstep = MIN(sys->maxstep,calc_sqrt_D0(tot2_norm2));
2272     sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2273     sys->varstep2.norm2/tot2_norm2;
2274     sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2275     sys->mulstep2.norm2/tot2_norm2;
2276    
2277     } else if( (tot2_norm2>0.0)&&(calc_sqrt_D0(tot2_norm2)<=sys->maxstep) ) {
2278     /* Attempt step in varstep2 direction */
2279     vars.alpha1 = 0.0;
2280     vars.alpha2 = 1.0;
2281     sys->maxstep = calc_sqrt_D0(tot2_norm2);
2282     sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2283     sys->varstep2.norm2/tot2_norm2;
2284     sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2285     sys->mulstep2.norm2/tot2_norm2;
2286    
2287     } else if( (tot2_norm2==0.0 || sys->s.block.current_size==1) &&
2288     (tot1_norm2 > 0.0) ) {
2289     /* Attempt step in varstep1 direction */
2290     vars.alpha1 = 1.0;
2291     vars.alpha2 = 0.0;
2292     if ( (sys->gamma.norm2/sys->Jgamma.norm2)*
2293     calc_sqrt_D0(sys->gamma.norm2) <= sys->maxstep )
2294     sys->maxstep = (sys->gamma.norm2/sys->Jgamma.norm2)*
2295     calc_sqrt_D0(sys->gamma.norm2);
2296     sys->varstep.norm2 = calc_sqr_D0(sys->maxstep)*
2297     sys->varstep1.norm2/tot1_norm2;
2298     sys->mulstep.norm2 = calc_sqr_D0(sys->maxstep)*
2299     sys->mulstep1.norm2/tot1_norm2;
2300    
2301     } else {
2302     /* Attempt step in varstep1-varstep2 direction */
2303     vars.parms.low = 0.0;
2304     vars.parms.high = MAXDOUBLE;
2305     vars.parms.guess = 1.0;
2306     calc_2x2_system(sys,&vars);
2307     do {
2308     coefs_from_parm(sys, &vars);
2309     adjust_parms(sys, &vars);
2310     } while( fabs(vars.error) > STEPSIZEERR_MAX &&
2311     vars.parms.high - vars.parms.low > PARMRNG_MIN );
2312     if (SHOW_LESS_IMPT) {
2313     FPRINTF(LIF(sys),"%-40s ---> %g\n",
2314     " parameter high", vars.parms.high);
2315     FPRINTF(LIF(sys),"%-40s ---> %g\n",
2316     " parameter low", vars.parms.low);
2317     FPRINTF(LIF(sys),"%-40s ---> %g\n",
2318     " Error in step length", vars.error);
2319     }
2320     sys->varstep.norm2 = step_norm2(sys, &vars);
2321     sys->mulstep.norm2 = 0.0;
2322     }
2323    
2324