/[ascend]/trunk/base/generic/solver/slv3.c
ViewVC logotype

Annotation of /trunk/base/generic/solver/slv3.c

Parent Directory Parent Directory | Revision Log Revision Log


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