/[ascend]/trunk/ipslv-temp/slv5.montero.c
ViewVC logotype

Annotation of /trunk/ipslv-temp/slv5.montero.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 916 - (hide annotations) (download) (as text)
Tue Oct 31 02:39:16 2006 UTC (17 years, 1 month ago) by johnpye
File MIME type: text/x-csrc
File size: 86427 byte(s)
Adding Vicente's IPSlv code (needs revising and testing against the new code)\
1 johnpye 916 /*
2     * IPSlv ASCEND Interior Point Method Solver
3     * by Vicente Rico-Ramirez based on QRSlv
4     * Created:
5     * Version: $ $
6     * Version control file: $ $
7     * Date last modified: $ $
8     * Last modified by: $Author: $
9     *
10     * This file is part of the SLV solver.
11     *
12     * The SLV solver is free software; you can redistribute
13     * it and/or modify it under the terms of the GNU General Public License as
14     * published by the Free Software Foundation; either version 2 of the
15     * License, or (at your option) any later version.
16     *
17     * The SLV solver is distributed in hope that it will be
18     * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
19     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20     * General Public License for more details.
21     *
22     * You should have received a copy of the GNU General Public License
23     * along with the program; if not, write to the Free Software Foundation,
24     * Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
25     * COPYING. COPYING is found in ../compiler.
26     *
27     */
28    
29     #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/time.h"
36     #include "utilities/mem.h"
37     #include "utilities/ascPanic.h"
38     #include "general/list.h"
39     #include "compiler/fractions.h"
40     #include "compiler/dimen.h"
41     #include "compiler/functype.h"
42     #include "compiler/func.h"
43     #include "solver/mtx.h"
44     #include "solver/linsol.h"
45     #include "solver/linsolqr.h"
46     #include "solver/slv_types.h"
47     #include "solver/var.h"
48     #include "solver/rel.h"
49     #include "solver/discrete.h"
50     #include "solver/conditional.h"
51     #include "solver/logrel.h"
52     #include "solver/bnd.h"
53     #include "solver/calc.h"
54     #include "solver/relman.h"
55     #include "solver/slv_common.h"
56     #include "solver/slv_client.h"
57     #include "solver/slv5.h"
58     #include "solver/slv_stdcalls.h"
59    
60     #if !defined(STATIC_IPSLV) && !defined(DYNAMIC_IPSLV)
61     int slv5_register(SlvFunctionsT *f)
62     {
63     (void)f; /* stop gcc whine about unused parameter */
64    
65     FPRINTF(ASCERR,"IPSlv not compiled in this ASCEND IV.\n");
66     return 1;
67     }
68     #else /* either STATIC_IPSLV or DYNAMIC_IPSLV is defined */
69     #ifdef DYNAMIC_IPSLV
70     /* do dynamic loading stuff. yeah, right */
71     #else /* following is used if STATIC_IPSLV is defined */
72    
73     #define DEBUG TRUE
74     #define DEBUG_COMPLEMENTARY_VAR TRUE
75     #define DEBUG_OBJ_VALUES TRUE
76     #define DEBUG_ITERATION TRUE
77     #define DEBUG_CENTERING TRUE
78    
79     #define SLV5(s) ((slv5_system_t)(s))
80     #define SERVER (sys->slv)
81     #define slv5_PA_SIZE 27 /* MUST INCREMENT WHEN ADDING PARAMETERS */
82     #define slv5_RA_SIZE 4
83    
84     /* do not delete (or extend) this array definition.
85     */
86     #define IEX(n) slv5_iaexpln[(n)]
87     #define slv5_IA_SIZE 10
88     static char *slv5_iaexpln[slv5_IA_SIZE] = {
89     "If lifds != 0 and showlessimportant is TRUE, show direct solve details",
90     "If savlin != 0, write out matrix data file at each iteration to SlvLinsol.dat",
91     "Cutoff is the block size cutoff for MODEL-based reordering of partitions",
92     "Update jacobian every this many major iterations",
93     "Update row scalings every this many major iterations",
94     "Update column scalings every this many major iterations",
95     "Check jacobian for poorly scaled columns and whine if found",
96     "Reorder option. 0 = MODEL based, 1 = MODEL based2, 2 = simple spk1",
97     "Use safe calculation routines",
98     "scaleopt = 0: 2norm"
99     };
100    
101     /* change slv5_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
102     #define UPDATE_JACOBIAN_PTR (sys->parm_array[0])
103     #define UPDATE_JACOBIAN ((*(int *)UPDATE_JACOBIAN_PTR))
104     #define UPDATE_WEIGHTS_PTR (sys->parm_array[1])
105     #define UPDATE_WEIGHTS ((*(int *)UPDATE_WEIGHTS_PTR))
106     #define UPDATE_NOMINALS_PTR (sys->parm_array[2])
107     #define UPDATE_NOMINALS ((*(int *)UPDATE_NOMINALS_PTR))
108     #define DUMPCNORM_PTR (sys->parm_array[3])
109     #define DUMPCNORM ((*(int *)DUMPCNORM_PTR))
110     #define SAFE_CALC_PTR (sys->parm_array[4])
111     #define SAFE_CALC ((*(int *)SAFE_CALC_PTR))
112     #define SCALEOPT_PTR (sys->parm_array[5])
113     #define SCALEOPT ((*(char **)SCALEOPT_PTR))
114     #define TOO_SMALL_PTR (sys->parm_array[6])
115     #define TOO_SMALL ((*(real64 *)TOO_SMALL_PTR))
116     #define CNLOW_PTR (sys->parm_array[7])
117     #define CNLOW ((*(real64 *)CNLOW_PTR))
118     #define CNHIGH_PTR (sys->parm_array[8])
119     #define CNHIGH ((*(real64 *)CNHIGH_PTR))
120     #define TOWARD_BOUNDS_PTR (sys->parm_array[9])
121     #define TOWARD_BOUNDS ((*(real64 *)TOWARD_BOUNDS_PTR))
122     #define IGNORE_BOUNDS_PTR (sys->parm_array[10])
123     #define IGNORE_BOUNDS ((*(int32 *)IGNORE_BOUNDS_PTR))
124     #define RHO_PTR (sys->parm_array[11])
125     #define RHO ((*(real64 *)RHO_PTR))
126     #define PARTITION_PTR (sys->parm_array[12])
127     #define PARTITION ((*(int32 *)PARTITION_PTR))
128     #define SHOW_LESS_IMPT_PTR (sys->parm_array[13])
129     #define SHOW_LESS_IMPT ((*(int32 *)SHOW_LESS_IMPT_PTR))
130     #define TIME_LIMIT_PTR (sys->parm_array[14])
131     #define TIME_LIMIT ((*(int32 *)TIME_LIMIT_PTR))
132     #define ITER_LIMIT_PTR (sys->parm_array[15])
133     #define ITER_LIMIT ((*(int32 *)ITER_LIMIT_PTR))
134     #define ALPHA_PTR (sys->parm_array[16])
135     #define ALPHA ((*(real64 *)ALPHA_PTR))
136     #define SING_TOL_PTR (sys->parm_array[17])
137     #define SING_TOL ((*(real64 *)SING_TOL_PTR))
138     #define PIVOT_TOL_PTR (sys->parm_array[18])
139     #define PIVOT_TOL ((*(real64 *)PIVOT_TOL_PTR))
140     #define FEAS_TOL_PTR (sys->parm_array[19])
141     #define FEAS_TOL ((*(real64 *)FEAS_TOL_PTR))
142     #define LIFDS_PTR (sys->parm_array[20])
143     #define LIFDS ((*(int32 *)LIFDS_PTR))
144     #define SAVLIN_PTR (sys->parm_array[21])
145     #define SAVLIN ((*(int32 *)SAVLIN_PTR))
146     #define REORDER_OPTION_PTR (sys->parm_array[22])
147     #define REORDER_OPTION ((*(char **)REORDER_OPTION_PTR))
148     #define CUTOFF_PTR (sys->parm_array[23])
149     #define CUTOFF ((*(int32 *)CUTOFF_PTR))
150     #define FACTOR_OPTION_PTR (sys->parm_array[24])
151     #define FACTOR_OPTION ((*(char **)FACTOR_OPTION_PTR))
152     #define CONVOPT_PTR (sys->parm_array[25])
153     #define CONVOPT ((*(char **)CONVOPT_PTR))
154     #define LINTIME_PTR (sys->parm_array[26])
155     #define LINTIME ((*(int *)LINTIME_PTR))
156     /* change slv5_PA_SIZE above (MUST INCREMENT) WHEN ADDING PARAMETERS */
157    
158    
159     #define REX(n) slv5_raexpln[(n)]
160     static
161     char *slv5_raexpln[slv5_RA_SIZE] = {
162     "Var nominal to use if user specifies 0.0",
163     "Smallest column norm we won't complain about if checking",
164     "Largest column norm we won't complain about if checking",
165     "If bound is in the way, we go this fraction toward it"};
166    
167     struct update_data {
168     int jacobian; /* Countdown on jacobian updating */
169     int weights; /* Countdown on weights updating */
170     int nominals; /* Countdown on nominals updating */
171     };
172    
173    
174     struct jacobian_data {
175     linsolqr_system_t sys; /* Linear system */
176     mtx_matrix_t mtx; /* Transpose gradient of residuals */
177     real64 *rhs; /* RHS from linear system */
178     dof_t *dofdata; /* dof data pointer from server */
179     mtx_region_t reg; /* Current block region */
180     int32 rank; /* Numerical rank of the jacobian */
181     enum factor_method fm; /* Linear factorization method */
182     boolean accurate; /* ? Recalculate matrix */
183     boolean singular; /* ? Can matrix be inverted */
184     boolean old_partition;/* old value of partition flag */
185     };
186    
187     struct slv5_system_structure {
188    
189     /*
190     * Problem definition
191     */
192     slv_system_t slv; /* slv_system_t back-link */
193     struct rel_relation *obj; /* Objective function: NULL = none */
194     struct var_variable **vlist; /* Variable list (NULL terminated) */
195     struct rel_relation **rlist; /* Relation list (NULL terminated) */
196    
197     /*
198     * Solver information
199     */
200     int integrity; /* ? Has the system been created */
201     int32 presolved; /* ? Has the system been presolved */
202     slv_parameters_t p; /* Parameters */
203     slv_status_t s; /* Status (as of iteration end) */
204     struct update_data update; /* Jacobian frequency counters */
205     int32 cap; /* Order of matrix/vectors */
206     int32 rank; /* Symbolic rank of problem */
207     int32 vused; /* Free and incident variables */
208     int32 vtot; /* length of varlist */
209     int32 rused; /* Included relations */
210     int32 rtot; /* length of rellist */
211     double clock; /* CPU time */
212     void *parm_array[slv5_PA_SIZE]; /* array of pointers to param values */
213     struct slv_parameter pa[slv5_PA_SIZE];/* &pa[0] => sys->p.parms */
214    
215     /*
216     * Calculated data (scaled)
217     */
218     struct jacobian_data J; /* linearized system */
219     struct vector_data nominals; /* Variable nominals */
220     struct vector_data weights; /* Relation weights */
221     struct vector_data variables; /* Variable values */
222     struct vector_data residuals; /* Relation residuals */
223     struct vector_data newton_residuals; /* Newton Relation residuals */
224     struct vector_data perturbed_residuals;/* Perturbed Relation residuals */
225     struct vector_data newton; /* Dependent variables */
226     struct vector_data perturbed_newton; /* Perturbed Newton direction */
227     struct vector_data varnewstep;/* hypothetical newton step in variables */
228     struct vector_data varstep; /* Step in variables */
229    
230     real64 progress; /* expected progress */
231     real64 sigma; /* penalty parameter */
232     real64 mu; /* Complementary gap */
233     real64 muaff; /* Complementary gap after newton */
234     real64 sigmamu; /* Complementary gap times penalty */
235     real64 eta; /* penalty of objective */
236     real64 nu; /* objective residuals */
237     real64 psi; /* Objective */
238     real64 comp; /* no. of complementary eqns. */
239     };
240    
241    
242     /*
243     * Integrity checks
244     * ----------------
245     * check_system(sys)
246     */
247    
248     #define OK ((int)813029392)
249     #define DESTROYED ((int)103289182)
250     static int check_system(slv5_system_t sys)
251     /*
252     * Checks sys for NULL and for integrity.
253     */
254     {
255     if( sys == NULL ) {
256     FPRINTF(ASCERR,"ERROR: (slv5) check_system\n");
257     FPRINTF(ASCERR," NULL system handle.\n");
258     return 1;
259     }
260    
261     switch( sys->integrity ) {
262     case OK:
263     return 0;
264     case DESTROYED:
265     FPRINTF(ASCERR,"ERROR: (slv5) check_system\n");
266     FPRINTF(ASCERR," System was recently destroyed.\n");
267     return 1;
268     default:
269     FPRINTF(ASCERR,"ERROR: (slv5) check_system\n");
270     FPRINTF(ASCERR," System reused or never allocated.\n");
271     return 1;
272     }
273     }
274    
275     /*
276     * General input/output routines
277     * -----------------------------
278     * print_var_name(out,sys,var)
279     * print_rel_name(out,sys,rel)
280     */
281    
282     #define print_var_name(a,b,c) slv_print_var_name((a),(b)->slv,(c))
283     #define print_rel_name(a,b,c) slv_print_rel_name((a),(b)->slv,(c))
284    
285     /*
286     * Debug output routines
287     * ---------------------
288     * debug_delimiter(fp)
289     * debug_out_vector(fp,vec)
290     * debug_out_var_values(fp,sys)
291     * debug_out_rel_residuals(fp,sys)
292     * debug_out_jacobian(fp,sys)
293     * debug_write_array(fp,real64 *,length)
294     * debug_out_parameters(fp)
295     */
296    
297     /*
298     * Outputs a hyphenated line.
299     */
300     static void debug_delimiter( FILE *fp)
301     {
302     int i;
303     for( i=0; i<60; i++ ) PUTC('-',fp);
304     PUTC('\n',fp);
305     }
306    
307     #if DEBUG
308     /*
309     * Outputs a vector.
310     */
311     static void debug_out_vector( FILE *fp, struct vector_data *vec)
312     {
313     int32 ndx;
314     FPRINTF(fp,"Norm = %g, Accurate = %s, Vector range = %d to %d\n",
315     calc_sqrt_D0(vec->norm2), vec->accurate?"TRUE":"FALSE",
316     vec->rng->low,vec->rng->high);
317     FPRINTF(fp,"Vector --> ");
318     for( ndx=vec->rng->low ; ndx<=vec->rng->high ; ++ndx )
319     FPRINTF(fp, "%g ", vec->vec[ndx]);
320     PUTC('\n',fp);
321     }
322    
323     /*
324     * Outputs all variable values in current block.
325     */
326     static void debug_out_var_values( FILE *fp, slv5_system_t sys)
327     {
328     int32 col;
329     struct var_variable *var;
330    
331     FPRINTF(fp,"Var values --> \n");
332     for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high ; col++ ) {
333     var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
334     print_var_name(fp,sys,var);
335     FPRINTF(fp, "\nI Lb Value Ub Scale Col INom\n");
336     FPRINTF(fp,"%d\t%.4g\t%.4g\t%.4g\t%.4g\t%d\t%.4g\n",
337     var_sindex(var),var_lower_bound(var),var_value(var),
338     var_upper_bound(var),var_nominal(var),
339     col,sys->nominals.vec[col]);
340     }
341     }
342    
343     /*
344     * Outputs all relation residuals in current block.
345     */
346     static void debug_out_rel_residuals( FILE *fp, slv5_system_t sys)
347     {
348     int32 row;
349    
350     FPRINTF(fp,"Rel residuals --> \n");
351     for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high ; row++ ) {
352     struct rel_relation *rel;
353     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
354     FPRINTF(fp," %g : ",rel_residual(rel));
355     print_rel_name(fp,sys,rel);
356     PUTC('\n',fp);
357     }
358     PUTC('\n',fp);
359     }
360    
361     /*
362     * Outputs permutation and values of the nonzero elements in the
363     * the jacobian matrix.
364     */
365     static void debug_out_jacobian( FILE *fp, slv5_system_t sys)
366     {
367     mtx_coord_t nz;
368     real64 value;
369    
370     nz.row = sys->J.reg.row.low;
371     for( ; nz.row <= sys->J.reg.row.high; ++(nz.row) ) {
372     FPRINTF(fp," Row %d (rel %d)\n", nz.row,
373     mtx_row_to_org(sys->J.mtx,nz.row));
374     nz.col = mtx_FIRST;
375     while( value = mtx_next_in_row(sys->J.mtx,&nz,&(sys->J.reg.col)),
376     nz.col != mtx_LAST ) {
377     FPRINTF(fp," Col %d (var %d) has value %g\n", nz.col,
378     mtx_col_to_org(sys->J.mtx,nz.col), value);
379     }
380     }
381     }
382    
383     #endif
384    
385     static void debug_write_array(FILE *fp,real64 *vec, int32 length)
386     {
387     int32 i;
388     for (i=0; i< length;i++)
389     FPRINTF(fp,"%.20g\n",vec[i]);
390     }
391    
392     static char savlinfilename[]="SlvLinsol.dat.\0";
393     static char savlinfilebase[]="SlvLinsol.dat.\0";
394     static int savlinnum=0;
395    
396     /* The number to postfix to savlinfilebase. increases with file accesses. */
397    
398     /*
399     * Array/vector operations
400     * ----------------------------
401     * destroy_array(p)
402     * create_array(len,type)
403     * zero_vector(vec)
404     * copy_vector(vec1,vec2)
405     * prod = inner_product(vec1,vec2)
406     * norm2 = square_norm(vec)
407     * matrix_product(mtx,vec,prod,scale,transpose)
408     */
409    
410     #define destroy_array(p) \
411     if( (p) != NULL ) ascfree((p))
412     #define create_array(len,type) \
413     ((len) > 0 ? (type *)ascmalloc((len)*sizeof(type)) : NULL)
414     #define create_zero_array(len,type) \
415     ((len) > 0 ? (type *)asccalloc((len),sizeof(type)) : NULL)
416    
417     #define zero_vector(v) slv_zero_vector(v)
418     #define copy_vector(v,t) slv_copy_vector((v),(t))
419     #define inner_product(v,u) slv_inner_product((v),(u))
420     #define square_norm(v) slv_square_norm(v)
421     #define matrix_product(m,v,p,s,t) slv_matrix_product((m),(v),(p),(s),(t))
422    
423     /*
424     * Calculation routines
425     * --------------------
426     * calc_eta(server)
427     * ok = calc_residuals(sys)
428     * ok = calc_mu(sys)
429     * ok = calc_newton_residuals(sys)
430     * ok = calc_muaff(sys)
431     * ok = calc_sigma(sys)
432     * ok = calc_sigmamu(sys)
433     * ok = calc_perturbed_residuals(sys)
434     * ok = calc_J(sys)
435     * calc_nominals(sys)
436     * calc_weights(sys)
437     * scale_J(sys)
438     * scale_variables(sys)
439     * scale_residuals(sys)
440     * scale_perturbed_residuals(sys)
441     * calc_pivots(sys)
442     * calc_rhs(sys)
443     * calc_newton(sys)
444     * calc_perturbed_newton(sys)
445     * calc_varnewstep(sys)
446     * calc_varstep(sys)
447     * calc_nu(sys)
448     * calc_psi(sys)
449     */
450    
451    
452     /*
453     * Calculates the penalty of objective function
454     */
455     static void calc_comp( slv5_system_t sys )
456     {
457     int32 row;
458     struct rel_relation *rel;
459     real64 comp;
460    
461     comp = 0.0;
462     row = sys->residuals.rng->low;
463     for( ; row <= sys->residuals.rng->high; row++ ) {
464     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
465     #if DEBUG
466     if (!rel) {
467     int r;
468     r=mtx_row_to_org(sys->J.mtx,row);
469     FPRINTF(ASCERR,"NULL relation found !!\n");
470     FPRINTF(ASCERR,"at row %d rel %d in calc_comp\n",(int)row,r);
471     FFLUSH(ASCERR);
472     }
473     #endif
474     if (rel_active(rel) && rel_included (rel) && rel_complementary(rel)) {
475     #if DEBUG
476     FPRINTF(ASCERR,"Complementary equation in slv5 \n");
477     #endif /* DEBUG */
478     comp = comp + 1.0;
479     }
480     }
481    
482     #if DEBUG_ITERATION
483     FPRINTF(ASCERR," No. of complementary eqns = %g \n",comp);
484     #endif /* DEBUG_ITERATION */
485     sys->comp = comp;
486     }
487    
488     /*
489     * Calculates the penalty of objective function
490     */
491     static void calc_eta( slv5_system_t sys)
492     {
493    
494     real64 eta;
495    
496     if (sys->comp == 0.0) {
497     eta = 0.0;
498     } else {
499     eta = sys->comp * sqrt(sys->comp);
500     }
501     #if DEBUG_ITERATION
502     FPRINTF(ASCERR," eta = %g \n",eta);
503     #endif /* DEBUG_ITERATION */
504     sys->eta = eta;
505     }
506    
507    
508     /*
509     * Calculates all of the residuals in the current block and computes
510     * the residual norm for block status. Returns true iff calculations
511     * preceded without error.
512     */
513     static boolean calc_residuals( slv5_system_t sys)
514     {
515     int32 row;
516     struct rel_relation *rel;
517     double time0;
518    
519     if( sys->residuals.accurate ) return TRUE;
520    
521     calc_ok = TRUE;
522     row = sys->residuals.rng->low;
523     time0=tm_cpu_time();
524     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
525     for( ; row <= sys->residuals.rng->high; row++ ) {
526     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
527     #if DEBUG
528     if (!rel) {
529     int r;
530     r=mtx_row_to_org(sys->J.mtx,row);
531     FPRINTF(ASCERR,"NULL relation found !!\n");
532     FPRINTF(ASCERR,"at row %d rel %d in calc_residuals\n",(int)row,r);
533     FFLUSH(ASCERR);
534     }
535     #endif
536     sys->residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
537    
538     if (strcmp(CONVOPT,"ABSOLUTE") == 0) {
539     relman_calc_satisfied(rel,FEAS_TOL);
540     } else if (strcmp(CONVOPT,"RELNOM_SCALE") == 0) {
541     relman_calc_satisfied_scaled(rel,FEAS_TOL);
542     }
543     }
544     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
545     sys->s.block.functime += (tm_cpu_time() -time0);
546     sys->s.block.funcs++;
547     square_norm( &(sys->residuals) );
548     sys->s.block.residual = calc_sqrt_D0(sys->residuals.norm2);
549     return(calc_ok);
550     }
551    
552     /*
553     * Calculates the complementary gap in the current point
554     */
555     static boolean calc_mu( slv5_system_t sys)
556     {
557    
558     int32 row;
559     struct rel_relation *rel;
560     real64 muk;
561    
562     muk = 0.0;
563     row = sys->residuals.rng->low;
564    
565     #if DEBUG_CENTERING
566     FPRINTF(ASCERR,"row low is = %d \n",row);
567     FPRINTF(ASCERR,"row high is = %d \n",sys->residuals.rng->high);
568     #endif /* DEBUG_CENTERING */
569    
570     for(row ; row <= sys->residuals.rng->high; row++ ) {
571     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
572     #if DEBUG
573     if (!rel) {
574     int r;
575     r = mtx_row_to_org(sys->J.mtx,row);
576     FPRINTF(ASCERR,"NULL relation found !!\n");
577     FPRINTF(ASCERR,"at row %d rel %d in calc_mu \n",(int)row,r);
578     FFLUSH(ASCERR);
579     }
580     #endif
581     if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
582     muk = muk + rel_residual(rel);
583     #if DEBUG
584     FPRINTF(ASCERR,"Complementary equation in calc_mu \n");
585     FPRINTF(ASCERR,"row = %d\n",row);
586     FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
587     FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
588     FPRINTF(ASCERR,"Partial muk is = %g \n",muk);
589     #endif /* DEBUG */
590     }
591     }
592    
593     if (sys->comp > 0.0) {
594     muk = muk / sys->comp;
595     } else {
596     muk = 0.0;
597     }
598    
599     sys->mu = muk;
600    
601     #if DEBUG_CENTERING
602     FPRINTF(ASCERR,"muk is = %g \n",sys->mu);
603     #endif /* DEBUG_CENTERING */
604    
605     return(calc_ok);
606     }
607    
608    
609     /*
610     * Calculates all of the residuals in the current block and computes
611     * the residual norm for block status after a hypothetical Newton
612     * step. Returns true iff calculations
613     * preceded without error.
614     */
615     static boolean calc_newton_residuals( slv5_system_t sys)
616     {
617     int32 row;
618     struct rel_relation *rel;
619    
620     if( sys->newton_residuals.accurate ) return TRUE;
621    
622     calc_ok = TRUE;
623     row = sys->newton_residuals.rng->low;
624     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
625     for( ; row <= sys->newton_residuals.rng->high; row++ ) {
626     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
627     #if DEBUG
628     if (!rel) {
629     int r;
630     r=mtx_row_to_org(sys->J.mtx,row);
631     FPRINTF(ASCERR,"NULL relation found !!\n");
632     FPRINTF(ASCERR,"at row %d rel %d in calc_newton_residuals\n",(int)row,r);
633     FFLUSH(ASCERR);
634     }
635     #endif
636     sys->newton_residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
637     }
638     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
639     square_norm( &(sys->newton_residuals) );
640     return(calc_ok);
641     }
642    
643     /*
644     * Calculates the complementary gap after a hypothetical Newton Step
645     */
646     static boolean calc_muaff( slv5_system_t sys)
647     {
648     int32 row;
649     struct rel_relation *rel;
650     real64 muaff;
651    
652     muaff = 0.0;
653     row = sys->newton_residuals.rng->low;
654     for(row ; row <= sys->newton_residuals.rng->high; row++ ) {
655     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
656     #if DEBUG
657     if (!rel) {
658     int r;
659     r=mtx_row_to_org(sys->J.mtx,row);
660     FPRINTF(ASCERR,"NULL relation found !!\n");
661     FPRINTF(ASCERR,"at row %d rel %d in calc_muaff\n",(int)row,r);
662     FFLUSH(ASCERR);
663     }
664     #endif
665     if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
666     muaff = muaff + sys->newton_residuals.vec[row];
667     #if DEBUG
668     FPRINTF(ASCERR,"Complementary equation in calc_muaff \n");
669     FPRINTF(ASCERR,"row = %d\n",row);
670     FPRINTF(ASCERR,"residual vector = %g\n",sys->newton_residuals.vec[row]);
671     FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
672     FPRINTF(ASCERR,"Partial muaff is = %g \n",muaff);
673     #endif /* DEBUG */
674     }
675     }
676    
677     if (sys->comp > 0.0) {
678     muaff = muaff / sys->comp;
679     } else {
680     muaff = 0.0;
681     }
682    
683     sys->muaff = muaff;
684    
685     #if DEBUG_CENTERING
686     FPRINTF(ASCERR,"muaff is = %g \n",sys->muaff);
687     #endif /* DEBUG_CENTERING */
688    
689     return(calc_ok);
690     }
691    
692    
693     /*
694     * Calculates the penalty parameter sigma
695     */
696     static boolean calc_sigma( slv5_system_t sys)
697     {
698     real64 sigma, frac;
699    
700     if ((sys->mu) > 0.0) {
701     frac = (sys->muaff) / (sys->mu);
702     sigma = (frac) * (frac) * (frac);
703     } else {
704     frac = 0.0;
705     sigma = 0.0;
706     }
707    
708     sys->sigma = sigma;
709    
710     #if DEBUG_CENTERING
711     FPRINTF(ASCERR,"sigma is = %g \n",sys->sigma);
712     #endif /* DEBUG_CENTERING */
713    
714     return(calc_ok);
715     }
716    
717     /*
718     * Calculates the penalty parameter time the complementary gap
719     */
720     static boolean calc_sigmamu( slv5_system_t sys)
721     {
722     real64 sigmamu;
723    
724     if ((sys->mu) > 0.0 && (sys->sigma) > 0.0) {
725     sigmamu = (sys->mu) * (sys->sigma);
726     } else {
727     sigmamu = 0.0;
728     }
729    
730     sys->sigmamu = sigmamu;
731    
732     #if DEBUG_CENTERING
733     FPRINTF(ASCERR,"sigma mu is = %g \n",sys->sigmamu);
734     #endif /* DEBUG_CENTERING */
735    
736     return(calc_ok);
737     }
738    
739    
740     /*
741     * Calculates the perturbed residuals in the current block and computes
742     * the perturbed residual norm. Returns true iff calculations
743     * preceded without error.
744     */
745     static boolean calc_perturbed_residuals( slv5_system_t sys)
746     {
747     int32 row;
748     struct rel_relation *rel;
749    
750     if( sys->perturbed_residuals.accurate ) return TRUE;
751    
752     calc_ok = TRUE;
753     row = sys->perturbed_residuals.rng->low;
754     Asc_SignalHandlerPush(SIGFPE,SIG_IGN);
755     for( ; row <= sys->perturbed_residuals.rng->high; row++ ) {
756     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
757     #if DEBUG
758     if (!rel) {
759     int r;
760     r=mtx_row_to_org(sys->J.mtx,row);
761     FPRINTF(ASCERR,"NULL relation found !!\n");
762     FPRINTF(ASCERR,"at row %d rel %d in calc_perturbed _residuals\n",
763     (int)row,r);
764     FFLUSH(ASCERR);
765     }
766     #endif
767     sys->perturbed_residuals.vec[row] = relman_eval(rel,&calc_ok,SAFE_CALC);
768     #if DEBUG
769     FPRINTF(ASCERR,"calc_perturbed_residuals \n");
770     FPRINTF(ASCERR,"row = %d\n",row);
771     FPRINTF(ASCERR,"residual vector = %g\n",
772     sys->perturbed_residuals.vec[row]);
773     FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
774     #endif /* DEBUG */
775     if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
776     sys->perturbed_residuals.vec[row] = sys->perturbed_residuals.vec[row] -
777     sys->sigmamu ;
778     #if DEBUG
779     FPRINTF(ASCERR,"residual vector - sigmamu = %g\n",
780     sys->perturbed_residuals.vec[row]);
781     #endif /* DEBUG */
782     }
783     }
784     Asc_SignalHandlerPop(SIGFPE,SIG_IGN);
785     square_norm( &(sys->perturbed_residuals) );
786     return(calc_ok);
787     }
788    
789     /*
790     * Calculates the current block of the jacobian.
791     * It is initially unscaled.
792     */
793     static boolean calc_J( slv5_system_t sys)
794     {
795     int32 row;
796     var_filter_t vfilter;
797     double time0;
798     real64 resid;
799    
800     if( sys->J.accurate )
801     return TRUE;
802    
803     calc_ok = TRUE;
804     vfilter.matchbits = (VAR_INBLOCK | VAR_ACTIVE);
805     vfilter.matchvalue = (VAR_INBLOCK | VAR_ACTIVE);
806     time0=tm_cpu_time();
807     mtx_clear_region(sys->J.mtx,&(sys->J.reg));
808     for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
809     struct rel_relation *rel;
810     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
811     relman_diffs(rel,&vfilter,sys->J.mtx,&resid,SAFE_CALC);
812     }
813     sys->s.block.jactime += (tm_cpu_time() - time0);
814     sys->s.block.jacs++;
815    
816     if( --(sys->update.nominals) <= 0 ) sys->nominals.accurate = FALSE;
817     if( --(sys->update.weights) <= 0 ) sys->weights.accurate = FALSE;
818    
819     linsolqr_matrix_was_changed(sys->J.sys);
820     return(calc_ok);
821     }
822    
823     /*
824     * Retrieves the nominal values of all of the block variables,
825     * insuring that they are all strictly positive.
826     */
827     static void calc_nominals( slv5_system_t sys)
828     {
829     int32 col;
830     FILE *fp = MIF(sys);
831     if( sys->nominals.accurate ) return;
832     fp = MIF(sys);
833     col = sys->nominals.rng->low;
834     if(strcmp(SCALEOPT,"NONE") == 0){
835     for( ; col <= sys->nominals.rng->high; col++ ) {
836     sys->nominals.vec[col] = 1;
837     }
838     } else {
839     for( ; col <= sys->nominals.rng->high; col++ ) {
840     struct var_variable *var;
841     real64 n;
842     var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
843     n = var_nominal(var);
844     if( n <= 0.0 ) {
845     if( n == 0.0 ) {
846     n = TOO_SMALL;
847     FPRINTF(fp,"ERROR: (slv5) calc_nominals\n");
848     FPRINTF(fp," Variable ");
849     print_var_name(fp,sys,var);
850     FPRINTF(fp," \nhas nominal value of zero.\n");
851     FPRINTF(fp," Resetting to %g.\n",n);
852     var_set_nominal(var,n);
853     } else {
854     n = -n;
855     FPRINTF(fp,"ERROR: (slv5) calc_nominals\n");
856     FPRINTF(fp," Variable ");
857     print_var_name(fp,sys,var);
858     FPRINTF(fp," \nhas negative nominal value.\n");
859     FPRINTF(fp," Resetting to %g.\n",n);
860     var_set_nominal(var,n);
861     }
862     }
863     #if DEBUG
864     FPRINTF(fp,"Column %d is ",col);
865     print_var_name(fp,sys,var);
866     FPRINTF(fp,"\nScaling of column %d is %g\n",col,n);
867     #endif
868     sys->nominals.vec[col] = n;
869     }
870     }
871     square_norm( &(sys->nominals) );
872     sys->update.nominals = UPDATE_NOMINALS;
873     sys->nominals.accurate = TRUE;
874     }
875    
876     /*
877     * Calculates the weights of all of the block relations
878     * to scale the rows of the Jacobian.
879     */
880     static void calc_weights( slv5_system_t sys)
881     {
882     mtx_coord_t nz;
883     real64 sum;
884    
885     if( sys->weights.accurate )
886     return;
887    
888     nz.row = sys->weights.rng->low;
889     if(strcmp(SCALEOPT,"NONE") == 0) {
890     for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
891     sys->weights.vec[nz.row] = 1;
892     }
893     } else if (strcmp(SCALEOPT,"ROW_2NORM") == 0 ) {
894     for( ; nz.row <= sys->weights.rng->high; (nz.row)++ ) {
895     sum=mtx_sum_sqrs_in_row(sys->J.mtx,nz.row,&(sys->J.reg.col));
896     sys->weights.vec[nz.row] = (sum>0.0) ? 1.0/calc_sqrt_D0(sum) : 1.0;
897     }
898     }
899     square_norm( &(sys->weights) );
900     sys->update.weights = UPDATE_WEIGHTS;
901     sys->residuals.accurate = FALSE;
902     sys->weights.accurate = TRUE;
903     }
904    
905     /*
906     * Scales the jacobian.
907     */
908     static void scale_J( slv5_system_t sys)
909     {
910     int32 row;
911     int32 col;
912    
913     if( sys->J.accurate ) return;
914    
915     calc_nominals(sys);
916     for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ )
917     mtx_mult_col(sys->J.mtx,col,sys->nominals.vec[col],&(sys->J.reg.row));
918    
919     calc_weights(sys);
920     for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ )
921     mtx_mult_row(sys->J.mtx,row,sys->weights.vec[row],&(sys->J.reg.col));
922     }
923    
924     static void jacobian_scaled(slv5_system_t sys)
925     {
926     int32 col;
927     if (DUMPCNORM) {
928     for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
929     real64 cnorm;
930     cnorm =
931     calc_sqrt_D0(mtx_sum_sqrs_in_col(sys->J.mtx,col,&(sys->J.reg.row)));
932     if (cnorm >CNHIGH || cnorm <CNLOW) {
933     FPRINTF(ASCERR,"[col %d org %d] %g\n", col,
934     mtx_col_to_org(sys->J.mtx,col), cnorm);
935     }
936     }
937     }
938    
939     sys->update.jacobian = UPDATE_JACOBIAN;
940     sys->J.accurate = TRUE;
941     sys->J.singular = FALSE; /* yet to be determined */
942     #if DEBUG
943     FPRINTF(LIF(sys),"\nJacobian: \n");
944     debug_out_jacobian(LIF(sys),sys);
945     #endif
946     }
947    
948     static void scale_variables( slv5_system_t sys)
949     {
950     int32 col;
951    
952     if( sys->variables.accurate ) return;
953    
954     col = sys->variables.rng->low;
955     for( ; col <= sys->variables.rng->high; col++ ) {
956     struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
957     sys->variables.vec[col] = var_value(var)/sys->nominals.vec[col];
958     }
959     square_norm( &(sys->variables) );
960     sys->variables.accurate = TRUE;
961     #if DEBUG
962     FPRINTF(LIF(sys),"Variables: ");
963     debug_out_vector(LIF(sys),&(sys->variables));
964     #endif
965     }
966    
967     /*
968     * Scales the previously calculated residuals.
969     */
970     static void scale_residuals( slv5_system_t sys)
971     {
972     int32 row;
973    
974     if( sys->residuals.accurate ) return;
975    
976     row = sys->residuals.rng->low;
977     for( ; row <= sys->residuals.rng->high; row++ ) {
978     struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
979     sys->residuals.vec[row] = rel_residual(rel)*sys->weights.vec[row];
980     }
981     square_norm( &(sys->residuals) );
982     sys->residuals.accurate = TRUE;
983     #if DEBUG
984     FPRINTF(LIF(sys),"Residuals: ");
985     debug_out_vector(LIF(sys),&(sys->residuals));
986     #endif
987     }
988    
989     /*
990     * Scales the previously calculated residuals.
991     */
992     static void scale_perturbed_residuals( slv5_system_t sys)
993     {
994     int32 row;
995    
996     if( sys->perturbed_residuals.accurate ) return;
997    
998     row = sys->perturbed_residuals.rng->low;
999     for( ; row <= sys->perturbed_residuals.rng->high; row++ ) {
1000     struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1001     #if DEBUG
1002     FPRINTF(ASCERR,"scale_perturbed_residuals \n");
1003     FPRINTF(ASCERR,"row = %d\n",row);
1004     FPRINTF(ASCERR,"residual vector = %g\n",
1005     sys->perturbed_residuals.vec[row]);
1006     FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1007     #endif /* DEBUG */
1008     sys->perturbed_residuals.vec[row] = sys->perturbed_residuals.vec[row]
1009     * sys->weights.vec[row];
1010     }
1011     square_norm( &(sys->perturbed_residuals) );
1012     sys->perturbed_residuals.accurate = TRUE;
1013     #if DEBUG
1014     FPRINTF(LIF(sys),"Perturbed Residuals: ");
1015     debug_out_vector(LIF(sys),&(sys->perturbed_residuals));
1016     #endif
1017     }
1018    
1019    
1020     /*
1021     * Scale system dependent on interface parameters
1022     */
1023     static void scale_perturbed_system( slv5_system_t sys )
1024     {
1025     if(strcmp(SCALEOPT,"NONE") == 0){
1026     scale_perturbed_residuals(sys);
1027     return;
1028     }
1029     if(strcmp(SCALEOPT,"ROW_2NORM") == 0){
1030     scale_perturbed_residuals(sys);
1031     return;
1032     }
1033     }
1034    
1035     /*
1036     * Scale system dependent on interface parameters
1037     */
1038     static void scale_system( slv5_system_t sys )
1039     {
1040     if(strcmp(SCALEOPT,"NONE") == 0){
1041     if(sys->J.accurate == FALSE){
1042     calc_nominals(sys);
1043     calc_weights(sys);
1044     jacobian_scaled(sys);
1045     }
1046     scale_variables(sys);
1047     scale_residuals(sys);
1048     return;
1049     }
1050     if(strcmp(SCALEOPT,"ROW_2NORM") == 0 ){
1051     if(sys->J.accurate == FALSE){
1052     scale_J(sys);
1053     jacobian_scaled(sys);
1054     }
1055     scale_variables(sys);
1056     scale_residuals(sys);
1057     return;
1058     }
1059     }
1060    
1061    
1062     /*
1063     * Obtain the equations and variables which
1064     * are able to be pivoted.
1065     * return value is the row rank deficiency, which we hope is 0.
1066     */
1067     static int calc_pivots(slv5_system_t sys)
1068     {
1069     int row_rank_defect=0, oldtiming;
1070     linsolqr_system_t lsys = sys->J.sys;
1071     FILE *fp = LIF(sys);
1072    
1073     oldtiming = g_linsolqr_timing;
1074     g_linsolqr_timing =LINTIME;
1075     linsolqr_factor(lsys,sys->J.fm); /* factor */
1076     g_linsolqr_timing = oldtiming;
1077    
1078     sys->J.rank = linsolqr_rank(lsys);
1079     sys->J.singular = FALSE;
1080     row_rank_defect = sys->J.reg.row.high -
1081     sys->J.reg.row.low+1 - sys->J.rank;
1082     if( row_rank_defect > 0 ) {
1083     int32 row,krow;
1084     mtx_sparse_t *uprows=NULL;
1085     sys->J.singular = TRUE;
1086     uprows = linsolqr_unpivoted_rows(lsys);
1087     if (uprows !=NULL) {
1088     for( krow=0; krow < uprows->len ; krow++ ) {
1089     int32 org_row;
1090     struct rel_relation *rel;
1091    
1092     org_row = uprows->idata[krow];
1093     row = mtx_org_to_row(sys->J.mtx,org_row);
1094     rel = sys->rlist[org_row];
1095     FPRINTF(fp,"%-40s ---> ","Relation not pivoted");
1096     print_rel_name(fp,sys,rel);
1097     PUTC('\n',fp);
1098    
1099     /*
1100     * assign zeros to the corresponding weights
1101     * so that subsequent calls to "scale_residuals"
1102     * will only measure the pivoted equations.
1103     */
1104     sys->weights.vec[row] = 0.0;
1105     sys->residuals.vec[row] = 0.0;
1106     sys->residuals.accurate = FALSE;
1107     sys->perturbed_residuals.vec[row] = 0.0;
1108     sys->perturbed_residuals.accurate = FALSE;
1109     sys->newton_residuals.vec[row] = 0.0;
1110     sys->newton_residuals.accurate = FALSE;
1111     mtx_mult_row(sys->J.mtx,row,0.0,&(sys->J.reg.col));
1112     }
1113     mtx_destroy_sparse(uprows);
1114     }
1115     if( !sys->residuals.accurate ) {
1116     square_norm( &(sys->residuals) );
1117     sys->residuals.accurate = TRUE;
1118     sys->update.weights = 0; /* re-compute weights next iteration. */
1119     }
1120     }
1121     if( sys->J.rank < sys->J.reg.col.high-sys->J.reg.col.low+1 ) {
1122     int32 col,kcol;
1123     mtx_sparse_t *upcols=NULL;
1124     if (NOTNULL(upcols)) {
1125     for( kcol=0; upcols != NULL && kcol < upcols->len ; kcol++ ) {
1126     int32 org_col;
1127     struct var_variable *var;
1128    
1129     org_col = upcols->idata[kcol];
1130     col = mtx_org_to_col(sys->J.mtx,org_col);
1131     var = sys->vlist[org_col];
1132     FPRINTF(fp,"%-40s ---> ","Variable not pivoted");
1133     print_var_name(fp,sys,var);
1134     PUTC('\n',fp);
1135     /*
1136     * If we're not optimizing (everything should be
1137     * pivotable) or this was one of the dependent variables,
1138     * consider this variable as if it were fixed.
1139     */
1140     if( col <= sys->J.reg.col.high ) {
1141     mtx_mult_col(sys->J.mtx,col,0.0,&(sys->J.reg.row));
1142     }
1143     }
1144     mtx_destroy_sparse(upcols);
1145     }
1146     }
1147     if (SHOW_LESS_IMPT) {
1148     FPRINTF(LIF(sys),"%-40s ---> %d (%s)\n","Jacobian rank", sys->J.rank,
1149     sys->J.singular ? "deficient":"full");
1150     FPRINTF(LIF(sys),"%-40s ---> %g\n","Smallest pivot",
1151     linsolqr_smallest_pivot(sys->J.sys));
1152     }
1153     return row_rank_defect;
1154     }
1155    
1156    
1157     /*
1158     * Calculates just the jacobian RHS. This function should be used to
1159     * supplement calculation of the jacobian. The vector vec must
1160     * already be calculated and scaled so as to simply be added to the
1161     * rhs. Caller is responsible for initially zeroing the rhs vector.
1162     */
1163     static void calc_rhs(slv5_system_t sys, struct vector_data *vec,
1164     real64 scalar, boolean transpose)
1165     {
1166     if( transpose ) { /* vec is indexed by col */
1167     int32 col;
1168     for( col=vec->rng->low; col<=vec->rng->high; col++ ) {
1169     sys->J.rhs[mtx_col_to_org(sys->J.mtx,col)] += scalar*vec->vec[col];
1170     }
1171     } else { /* vec is indexed by row */
1172     int32 row;
1173     for( row=vec->rng->low; row<=vec->rng->high; row++ ) {
1174     sys->J.rhs[mtx_row_to_org(sys->J.mtx,row)] += scalar*vec->vec[row];
1175     }
1176     }
1177     linsolqr_rhs_was_changed(sys->J.sys,sys->J.rhs);
1178     }
1179    
1180     /*
1181     * Computes a step to solve the linearized equations.
1182     */
1183     static void calc_newton( slv5_system_t sys)
1184     {
1185     linsolqr_system_t lsys = sys->J.sys;
1186     int32 col;
1187    
1188     if( sys->newton.accurate )
1189     return;
1190    
1191     sys->J.rhs = linsolqr_get_rhs(lsys,0);
1192     mtx_zero_real64(sys->J.rhs,sys->cap);
1193     calc_rhs(sys, &(sys->residuals), -1.0, FALSE);
1194     linsolqr_solve(lsys,sys->J.rhs);
1195     col = sys->newton.rng->low;
1196     for( ; col <= sys->newton.rng->high; col++ ) {
1197     sys->newton.vec[col] =
1198     linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1199     }
1200     if (SAVLIN) {
1201     FILE *ldat;
1202     int32 ov;
1203     sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1204     ldat=fopen(savlinfilename,"w");
1205     FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1206     sys->s.iteration);
1207     debug_write_array(ldat,sys->J.rhs,sys->cap);
1208     FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1209     for(ov=0 ; ov < sys->cap; ov++ )
1210     FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1211     fclose(ldat);
1212     }
1213     square_norm( &(sys->newton) );
1214     sys->newton.accurate = TRUE;
1215     #if DEBUG
1216     FPRINTF(LIF(sys),"Newton: ");
1217     debug_out_vector(LIF(sys),&(sys->newton));
1218     #endif
1219     }
1220    
1221     /*
1222     * Computes a step to solve the linearized equations.
1223     */
1224     static void calc_perturbed_newton( slv5_system_t sys)
1225     {
1226     linsolqr_system_t lsys = sys->J.sys;
1227     int32 col;
1228    
1229     if( sys->perturbed_newton.accurate )
1230     return;
1231    
1232     sys->J.rhs = linsolqr_get_rhs(lsys,1);
1233     mtx_zero_real64(sys->J.rhs,sys->cap);
1234     calc_rhs(sys, &(sys->perturbed_residuals), -1.0, FALSE);
1235     linsolqr_solve(lsys,sys->J.rhs);
1236     col = sys->perturbed_newton.rng->low;
1237     for( ; col <= sys->perturbed_newton.rng->high; col++ ) {
1238     sys->perturbed_newton.vec[col] =
1239     linsolqr_var_value(lsys,sys->J.rhs,mtx_col_to_org(sys->J.mtx,col));
1240     }
1241     if (SAVLIN) {
1242     FILE *ldat;
1243     int32 ov;
1244     sprintf(savlinfilename,"%s%d",savlinfilebase,savlinnum++);
1245     ldat=fopen(savlinfilename,"w");
1246     FPRINTF(ldat,"================= resids (orgrowed) itn %d =====\n",
1247     sys->s.iteration);
1248     debug_write_array(ldat,sys->J.rhs,sys->cap);
1249     FPRINTF(ldat,"================= vars (orgcoled) ============\n");
1250     for(ov=0 ; ov < sys->cap; ov++ )
1251     FPRINTF(ldat,"%.20g\n",linsolqr_var_value(lsys,sys->J.rhs,ov));
1252     fclose(ldat);
1253     }
1254     square_norm( &(sys->perturbed_newton) );
1255     sys->perturbed_newton.accurate = TRUE;
1256     #if DEBUG
1257     FPRINTF(LIF(sys),"Perturbed Newton: ");
1258     debug_out_vector(LIF(sys),&(sys->perturbed_newton));
1259     #endif
1260     }
1261    
1262    
1263    
1264     /*
1265     * Calculate the perturbed descent direction
1266     * in the variables.
1267     */
1268     static void calc_varstep( slv5_system_t sys)
1269     {
1270     if( sys->varstep.accurate )
1271     return;
1272     copy_vector(&(sys->perturbed_newton),&(sys->varstep));
1273     sys->varstep.norm2 = sys->perturbed_newton.norm2;
1274    
1275     sys->varstep.accurate = TRUE;
1276     #if DEBUG
1277     FPRINTF(LIF(sys),"Varstep: ");
1278     debug_out_vector(LIF(sys),&(sys->varstep));
1279     #endif
1280     }
1281    
1282     /*
1283     * Calculate the hypothetical netown direction
1284     * in the variables.
1285     */
1286     static void calc_varnewstep( slv5_system_t sys)
1287     {
1288     if( sys->varnewstep.accurate )
1289     return;
1290     copy_vector(&(sys->newton),&(sys->varnewstep));
1291     sys->varnewstep.norm2 = sys->newton.norm2;
1292    
1293     sys->varnewstep.accurate = TRUE;
1294     #if DEBUG
1295     FPRINTF(LIF(sys),"Varnewstep: ");
1296     debug_out_vector(LIF(sys),&(sys->varnewstep));
1297     #endif
1298     }
1299    
1300     /*
1301     * Computes the error nu.
1302     */
1303     static void calc_nu( slv5_system_t sys)
1304     {
1305     struct rel_relation *rel;
1306     int32 row, r;
1307     real64 error;
1308    
1309     error = 0.0;
1310     row = sys->residuals.rng->low;
1311     for(row = sys->residuals.rng->low; row <= sys->residuals.rng->high; row++) {
1312     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1313     if (!rel) {
1314     r=mtx_row_to_org(sys->J.mtx,row);
1315     FPRINTF(ASCERR,"NULL relation found !!\n");
1316     FPRINTF(ASCERR,"at row %d rel %d in calc_nu\n",(int)row,r);
1317     FFLUSH(ASCERR);
1318     }
1319     if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
1320     #if DEBUG
1321     FPRINTF(ASCERR,"Complementary equation in calc_nu \n");
1322     FPRINTF(ASCERR,"row = %d\n",row);
1323     FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
1324     FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1325     #endif /* DEBUG */
1326     error = error + rel_residual(rel);
1327     } else {
1328     if (rel_active(rel) && rel_included(rel)) {
1329     error = error + (rel_residual(rel) *rel_residual(rel) );
1330     #if DEBUG
1331     FPRINTF(ASCERR,"Non complementary equation calc_nu \n");
1332     FPRINTF(ASCERR,"row = %d\n",row);
1333     FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
1334     FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1335     #endif /* DEBUG */
1336     }
1337     }
1338     }
1339    
1340     sys->nu = error;
1341    
1342     #if DEBUG_OBJ_VALUES
1343     FPRINTF(ASCERR," Error nu = %g \n",sys->nu);
1344     #endif /* DEBUG_OBJ_VALUES */
1345     }
1346    
1347    
1348     /*
1349     * Computes the objective function Psi.
1350     */
1351     static void calc_psi( slv5_system_t sys)
1352     {
1353    
1354     struct rel_relation *rel;
1355     int32 row, r;
1356     real64 sum, sumt;
1357    
1358     calc_nu( sys);
1359    
1360     sum = 0.0;
1361     for(row = sys->residuals.rng->low; row <= sys->residuals.rng->high; row++) {
1362     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1363     if (!rel) {
1364     r = mtx_row_to_org(sys->J.mtx,row);
1365     FPRINTF(ASCERR,"NULL relation found !!\n");
1366     FPRINTF(ASCERR,"at row %d rel %d in calc_psi\n",(int)row,r);
1367     FFLUSH(ASCERR);
1368     }
1369     if (rel_complementary(rel) && rel_active(rel) && rel_included(rel)) {
1370     #if DEBUG
1371     FPRINTF(ASCERR,"Complementary equation in calc_psi\n");
1372     FPRINTF(ASCERR,"row = %d\n",row);
1373     FPRINTF(ASCERR,"residual vector = %g\n",sys->residuals.vec[row]);
1374     FPRINTF(ASCERR,"rel_residual = %g \n",rel_residual(rel));
1375     #endif /* DEBUG */
1376     sumt = log10(rel_residual(rel));
1377     sum = sum + sumt;
1378     }
1379     }
1380    
1381     if (sys->eta > 0.0) {
1382     sys->psi = ( sys->eta * log10(sys->nu) ) - sum;
1383     } else {
1384     sys->psi = 0.5*sys->residuals.norm2;
1385     }
1386    
1387     #if DEBUG_OBJ_VALUES
1388     FPRINTF(ASCERR," Psi = %g \n",sys->psi);
1389     #endif /* DEBUG_OBJ_VALUES */
1390     }
1391    
1392    
1393    
1394     /*
1395     * Variable values maintenance
1396     * ---------------------------
1397     * restore_variables(sys)
1398     * coef = required_coef_to_stay_inbounds(sys)
1399     * apply_step(sys,coef)
1400     * step_accepted(sys)
1401     */
1402    
1403     /*
1404     * Restores the values of the variables before applying
1405     * a step.
1406     */
1407     static void restore_variables( slv5_system_t sys)
1408     {
1409     int32 col;
1410     real64 *vec;
1411     vec = (sys->nominals.vec);
1412     for( col = sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
1413     struct var_variable *var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1414     var_set_value(var,sys->variables.vec[col]*vec[col]);
1415     }
1416     }
1417    
1418    
1419     /*
1420     * Calculates the maximum fraction of the step which can be
1421     * taken without making negative the complementary vars.
1422     * It is assumed that the current complementary variable
1423     * is positive. The step must be calculated.
1424     */
1425     static real64 factor_for_complementary_vars( slv5_system_t sys, int32 v)
1426     {
1427     real64 factor, minfactor;
1428     struct var_variable *var;
1429     real64 dx,val,bnd;
1430     int32 col;
1431     struct vector_data step;
1432     real64 *vec;
1433    
1434     vec = (sys->nominals.vec);
1435    
1436     if (v == 1) {
1437     step = sys->varstep;
1438     } else {
1439     step = sys->varnewstep;
1440     }
1441    
1442     minfactor = 1.0;
1443     factor = 1.0;
1444     for( col=step.rng->low; col <= step.rng->high; col++ ) {
1445     var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1446     val = var_value(var);
1447     if (var_complementary(var) && (!var_fixed(var))) {
1448     dx = step.vec[col] * vec[col];
1449     bnd = 0.0;
1450     if( val + dx < bnd ) {
1451     factor = MIN((bnd-val)/dx, 1.0);
1452     if (factor < 1.0) {
1453     #if DEBUG_COMPLEMENTARY_VAR
1454     FPRINTF(ASCERR,"Negative Complementary Variable : \n");
1455     print_var_name(ASCERR,sys,var);
1456     FPRINTF(ASCERR,"\n");
1457     FPRINTF(ASCERR,"factor = %g \n",factor);
1458     #endif /* DEBUG_COMPLEMENTARY_VAR */
1459     }
1460     if( factor < minfactor ) {
1461     minfactor = factor;
1462     }
1463     }
1464     }
1465     }
1466     #if DEBUG_COMPLEMENTARY_VAR
1467     FPRINTF(ASCERR,"Minimum complementary factor = %g \n",minfactor);
1468     #endif /* DEBUG_COMPLEMENTARY_VAR */
1469     return minfactor;
1470     }
1471    
1472    
1473     /*
1474     * Adds step to the variable values in block. It projects
1475     * non complementary varaibles near bounds.
1476     */
1477     static void apply_step( slv5_system_t sys, int32 v, real64 factor)
1478     {
1479     FILE *lif = LIF(sys);
1480     struct var_variable *var;
1481     real64 dx,val,bnd;
1482     real64 bounds_coef;
1483     struct vector_data step;
1484     int32 col;
1485     real64 *vec;
1486    
1487     vec = (sys->nominals.vec);
1488     bounds_coef = 1.0;
1489    
1490     if (v == 1) {
1491     step = sys->varstep;
1492     } else {
1493     step = sys->varnewstep;
1494     }
1495    
1496     for(col=step.rng->low; col<=step.rng->high;col++) {
1497     var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1498     dx = vec[col]*step.vec[col];
1499     val = var_value(var);
1500     dx = dx*factor;
1501     if( val + dx > (bnd=var_upper_bound(var)) ) {
1502     dx = TOWARD_BOUNDS*(bnd-val);
1503     if (SHOW_LESS_IMPT) {
1504     FPRINTF(lif,"%-40s ---> ",
1505     " Variable projected to upper bound");
1506     print_var_name(lif,sys,var); PUTC('\n',lif);
1507     }
1508     } else if( val + dx < (bnd=var_lower_bound(var)) ) {
1509     dx = TOWARD_BOUNDS*(bnd-val);
1510     if (SHOW_LESS_IMPT) {
1511     FPRINTF(lif,"%-40s ---> ",
1512     " Variable projected to lower bound");
1513     print_var_name(lif,sys,var); PUTC('\n',lif);
1514     }
1515     }
1516     var_set_value(var,val+dx);
1517     }
1518    
1519     /* Allow weighted residuals to be recalculated at new point */
1520     sys->residuals.accurate = FALSE;
1521     sys->newton_residuals.accurate = FALSE;
1522     sys->perturbed_residuals.accurate = FALSE;
1523    
1524     return;
1525     }
1526    
1527    
1528     /*
1529     * This function should be called when the step is accepted.
1530     */
1531     static void step_accepted( slv5_system_t sys)
1532     {
1533     /* Maintain update status on jacobian and weights */
1534     if (--(sys->update.jacobian) <= 0) {
1535     sys->J.accurate = FALSE;
1536     }
1537    
1538     sys->variables.accurate = FALSE;
1539     sys->newton_residuals.accurate = FALSE;
1540     sys->perturbed_residuals.accurate = FALSE;
1541     sys->newton.accurate = FALSE;
1542     sys->perturbed_newton.accurate = FALSE;
1543     sys->varnewstep.accurate = FALSE;
1544     sys->varstep.accurate = FALSE;
1545     }
1546    
1547    
1548     /*
1549     * Block routines
1550     * --------------
1551     * feas = block_feasible(sys)
1552     * move_to_next_block(sys)
1553     * find_next_unconverged_block(sys)
1554     */
1555    
1556     /*
1557     * Returns TRUE if the current block is feasible, FALSE otherwise.
1558     * It is assumed that the residuals have been computed.
1559     */
1560     static boolean block_feasible( slv5_system_t sys)
1561     {
1562     int32 row;
1563    
1564     if( !sys->s.calc_ok )
1565     return(FALSE);
1566    
1567     for( row = sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
1568     struct rel_relation *rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1569     if( !rel_satisfied(rel) ) return FALSE;
1570     }
1571     return TRUE;
1572     }
1573    
1574     /*
1575     * Moves on to the next block, updating all of the solver information.
1576     * To move to the first block, set sys->s.block.current_block to -1 before
1577     * calling. If already at the last block, then sys->s.block.current_block
1578     * will equal the number of blocks and the system will be declared
1579     * converged. Otherwise, the residuals for the new block will be computed
1580     * and sys->s.calc_ok set according.
1581     */
1582     static void move_to_next_block( slv5_system_t sys)
1583     {
1584     struct var_variable *var;
1585     struct rel_relation *rel;
1586     int32 row;
1587     int32 col;
1588     int32 ci;
1589    
1590     if( sys->s.block.current_block >= 0 ) {
1591    
1592     /* Record cost accounting info here. */
1593     ci=sys->s.block.current_block;
1594     sys->s.cost[ci].size = sys->s.block.current_size;
1595     sys->s.cost[ci].iterations = sys->s.block.iteration;
1596     sys->s.cost[ci].funcs = sys->s.block.funcs;
1597     sys->s.cost[ci].jacs = sys->s.block.jacs;
1598     sys->s.cost[ci].functime = sys->s.block.functime;
1599     sys->s.cost[ci].jactime = sys->s.block.jactime;
1600     sys->s.cost[ci].time = sys->s.block.cpu_elapsed;
1601     sys->s.cost[ci].resid = sys->s.block.residual;
1602    
1603     /* De-initialize previous block */
1604     if (SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
1605     LIFDS)) {
1606     FPRINTF(LIF(sys),"Block %d converged.\n",
1607     sys->s.block.current_block);
1608     }
1609     for( col=sys->J.reg.col.low; col <= sys->J.reg.col.high; col++ ) {
1610     var = sys->vlist[mtx_col_to_org(sys->J.mtx,col)];
1611     var_set_in_block(var,FALSE);
1612     }
1613     for( row=sys->J.reg.row.low; row <= sys->J.reg.row.high; row++ ) {
1614     rel = sys->rlist[mtx_row_to_org(sys->J.mtx,row)];
1615     rel_set_in_block(rel,FALSE);
1616     }
1617     sys->s.block.previous_total_size += sys->s.block.current_size;
1618     }
1619    
1620     sys->s.block.current_block++;
1621     if( sys->s.block.current_block < sys->s.block.number_of ) {
1622     boolean ok;
1623    
1624     /* Initialize next block */
1625    
1626     sys->J.reg =
1627     (slv_get_solvers_blocks(SERVER))->block[sys->s.block.current_block];
1628    
1629    
1630     row = sys->J.reg.row.high - sys->J.reg.row.low + 1;
1631     col = sys->J.reg.col.high - sys->J.reg.col.low + 1;
1632     sys->s.block.current_size = MAX(row,col);
1633    
1634     sys->s.block.iteration = 0;
1635     sys->s.block.cpu_elapsed = 0.0;
1636     sys->s.block.functime = 0.0;
1637     sys->s.block.jactime = 0.0;
1638     sys->s.block.funcs = 0;
1639     sys->s.block.jacs = 0;
1640    
1641     if(SHOW_LESS_IMPT && (LIFDS ||
1642     sys->s.block.current_size > 1)) {
1643     debug_delimiter(LIF(sys));
1644     debug_delimiter(LIF(sys));
1645     }
1646     if(SHOW_LESS_IMPT && LIFDS) {
1647     FPRINTF(LIF(sys),"\n%-40s ---> %d in [%d..%d]\n",
1648     "Current block number", sys->s.block.current_block,
1649     0, sys->s.block.number_of-1);
1650     FPRINTF(LIF(sys),"%-40s ---> %d\n", "Current block size",
1651     sys->s.block.current_size);
1652     }
1653     sys->s.calc_ok = TRUE;
1654    
1655     if (!(sys->p.ignore_bounds) ) {
1656     slv_insure_bounds(SERVER, sys->J.reg.col.low,
1657     sys->J.reg.col.high,MIF(sys));
1658     }
1659    
1660     sys->residuals.accurate = FALSE;
1661     if( !(ok = calc_residuals(sys)) ) {
1662     FPRINTF(MIF(sys),
1663     "Residual calculation errors detected in move_to_next_block.\n");
1664     }
1665    
1666     /*
1667     * Update number of complementary equations
1668     */
1669     calc_comp(sys);
1670     calc_eta(sys);
1671    
1672     if( SHOW_LESS_IMPT &&
1673     (sys->s.block.current_size >1 ||
1674     LIFDS) ) {
1675     FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
1676     sys->s.block.residual);
1677     }
1678     sys->s.calc_ok = sys->s.calc_ok && ok;
1679    
1680     /* Must be updated as soon as required */
1681     sys->J.accurate = FALSE;
1682     sys->update.weights = 0;
1683     sys->update.nominals = 0;
1684     sys->variables.accurate = FALSE;
1685     sys->newton_residuals.accurate = FALSE;
1686     sys->perturbed_residuals.accurate = FALSE;
1687     sys->newton.accurate = FALSE;
1688     sys->perturbed_newton.accurate = FALSE;
1689     sys->varnewstep.accurate = FALSE;
1690     sys->varstep.accurate = FALSE;
1691    
1692     } else {
1693     boolean ok;
1694     /*
1695     * Before we claim convergence, we must check if we left behind
1696     * some unassigned relations. If and only if they happen to be
1697     * satisfied at the current point, convergence has been obtained.
1698     *
1699     * Also insures that all included relations have valid residuals.
1700     * Included inequalities will have correct residuals.
1701     * Unsatisfied included inequalities cause inconsistency.
1702     *
1703     * This of course ignores that fact an objective function might
1704     * be present. Then, feasibility isn't enough, is it now.
1705     */
1706     if( sys->s.struct_singular ) {
1707     /* black box w/singletons provoking bug here, maybe */
1708     sys->s.block.current_size = sys->rused - sys->rank;
1709     if(SHOW_LESS_IMPT) {
1710     debug_delimiter(LIF(sys));
1711     FPRINTF(LIF(sys),"%-40s ---> %d\n", "Unassigned Relations",
1712     sys->s.block.current_size);
1713     }
1714     sys->J.reg.row.low = sys->J.reg.col.low = sys->rank;
1715     sys->J.reg.row.high = sys->J.reg.col.high = sys->rused - 1;
1716     sys->residuals.accurate = FALSE;
1717     if( !(ok=calc_residuals(sys)) ) {
1718     FPRINTF(MIF(sys),
1719     "Residual calculation errors detected in leftover equations.\n");
1720     }
1721     if(SHOW_LESS_IMPT) {
1722     FPRINTF(LIF(sys),"%-40s ---> %g\n", "Residual norm (unscaled)",
1723     sys->s.block.residual);
1724     }
1725     if( block_feasible(sys) ) {
1726     if(SHOW_LESS_IMPT) {
1727     FPRINTF(LIF(sys),"\nUnassigned relations ok. You lucked out.\n");
1728     }
1729     sys->s.converged = TRUE;
1730     } else {
1731     if(SHOW_LESS_IMPT) {
1732     FPRINTF(LIF(sys),"\nProblem inconsistent: %s.\n",
1733     "Unassigned relations not satisfied");
1734     }
1735     sys->s.inconsistent = TRUE;
1736     }
1737     if(SHOW_LESS_IMPT) {
1738     debug_delimiter(LIF(sys));
1739     }
1740     } else {
1741     sys->s.converged = TRUE;
1742     }
1743     }
1744     }
1745    
1746    
1747     /*
1748     * Calls the appropriate reorder function on a block
1749     */
1750     static void reorder_new_block(slv5_system_t sys)
1751     {
1752     int32 method;
1753     if( sys->s.block.current_block < sys->s.block.number_of ) {
1754     if (strcmp(REORDER_OPTION,"SPK1") == 0) {
1755     method = 2;
1756     } else {
1757     method = 1;
1758     }
1759    
1760     if( sys->s.block.current_block <= sys->s.block.current_reordered_block &&
1761     sys->s.cost[sys->s.block.current_block].reorder_method == method &&
1762     sys->s.block.current_block >= 0 ) {
1763     #if DEBUG
1764     FPRINTF(ASCERR,"YOU JUST AVOIDED A REORDERING\n");
1765     #endif
1766     slv_set_up_block(SERVER,sys->s.block.current_block);
1767     /* tell linsol to bless it and get on with things */
1768     linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
1769     return; /*must have been reordered since last system build*/
1770     }
1771    
1772     /* Let the slv client function take care of reordering things
1773     * and setting in block flags.
1774     */
1775     if (strcmp(REORDER_OPTION,"SPK1") == 0) {
1776     sys->s.cost[sys->s.block.current_block].reorder_method = 2;
1777     slv_spk1_reorder_block(SERVER,sys->s.block.current_block,1);
1778     } else if (strcmp(REORDER_OPTION,"TEAR_DROP") == 0) {
1779     sys->s.cost[sys->s.block.current_block].reorder_method = 1;
1780     slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
1781     CUTOFF,
1782     0,mtx_SPK1);
1783     /* khack: try tspk1 for transpose case */
1784     } else if (strcmp(REORDER_OPTION,"OVER_TEAR") == 0) {
1785     sys->s.cost[sys->s.block.current_block].reorder_method = 1;
1786     slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
1787     CUTOFF,
1788     1,mtx_SPK1);
1789     } else {
1790     sys->s.cost[sys->s.block.current_block].reorder_method = 1;
1791     FPRINTF(MIF(sys),"IPSlv called with unknown reorder option\n");
1792     FPRINTF(MIF(sys),"IPSlv using single edge tear drop (TEAR_DROP).\n");
1793     slv_tear_drop_reorder_block(SERVER,sys->s.block.current_block,
1794     CUTOFF,0,mtx_SPK1);
1795     }
1796     /* tell linsol to bless it and get on with things */
1797     linsolqr_reorder(sys->J.sys,&(sys->J.reg),natural);
1798     if (sys->s.block.current_block > sys->s.block.current_reordered_block) {
1799     sys->s.block.current_reordered_block = sys->s.block.current_block;
1800     }
1801     }
1802     }
1803    
1804     /*
1805     * Moves to next unconverged block, assuming that the current block has
1806     * converged (or is -1, to start).
1807     */
1808     static void find_next_unconverged_block( slv5_system_t sys)
1809     {
1810     do {
1811     move_to_next_block(sys);
1812     #if DEBUG
1813     debug_out_var_values(ASCERR,sys);
1814     debug_out_rel_residuals(ASCERR,sys);
1815     #endif
1816     } while( !sys->s.converged && block_feasible(sys));
1817     reorder_new_block(sys);
1818     }
1819    
1820    
1821     /*
1822     * Iteration begin/end routines
1823     * ----------------------------
1824     * iteration_begins(sys)
1825     * iteration_ends(sys)
1826     */
1827    
1828     /*
1829     * Prepares sys for entering an iteration, increasing the iteration counts
1830     * and starting the clock.
1831     */
1832     static void iteration_begins( slv5_system_t sys)
1833     {
1834     sys->clock = tm_cpu_time();
1835     ++(sys->s.block.iteration);
1836     ++(sys->s.iteration);
1837     if(SHOW_LESS_IMPT&& (sys->s.block.current_size >1 ||
1838     LIFDS)) {
1839     FPRINTF(LIF(sys),"\n%-40s ---> %d\n",
1840     "Iteration", sys->s.block.iteration);
1841     FPRINTF(LIF(sys),"%-40s ---> %d\n",
1842     "Total iteration", sys->s.iteration);
1843     }
1844     }
1845    
1846     /*
1847     * Prepares sys for exiting an iteration, stopping the clock and recording
1848     * the cpu time.
1849     */
1850     static void iteration_ends( slv5_system_t sys)
1851     {
1852     double cpu_elapsed; /* elapsed this iteration */
1853    
1854     cpu_elapsed = (double)(tm_cpu_time() - sys->clock);
1855     sys->s.block.cpu_elapsed += cpu_elapsed;
1856     sys->s.cpu_elapsed += cpu_elapsed;
1857     if(SHOW_LESS_IMPT && (sys->s.block.current_size >1 ||
1858     LIFDS)) {
1859     FPRINTF(LIF(sys),"%-40s ---> %g\n",
1860     "Elapsed time", sys->s.block.cpu_elapsed);
1861     FPRINTF(LIF(sys),"%-40s ---> %g\n",
1862     "Total elapsed time", sys->s.cpu_elapsed);
1863     }
1864     }
1865    
1866     /*
1867     * Updates the solver status.
1868     */
1869     static void update_status( slv5_system_t sys)
1870     {
1871     boolean unsuccessful;
1872    
1873     if( !sys->s.converged ) {
1874     sys->s.time_limit_exceeded =
1875     (sys->s.block.cpu_elapsed >= TIME_LIMIT);
1876     sys->s.iteration_limit_exceeded =
1877     (sys->s.block.iteration >= ITER_LIMIT);
1878     }
1879    
1880     unsuccessful = sys->s.diverged || sys->s.inconsistent ||
1881     sys->s.iteration_limit_exceeded || sys->s.time_limit_exceeded;
1882    
1883     sys->s.ready_to_solve = !unsuccessful && !sys->s.converged;
1884     sys->s.ok = !unsuccessful && sys->s.calc_ok && !sys->s.struct_singular;
1885     }
1886    
1887     static
1888     int32 slv5_get_default_parameters(slv_system_t server, SlvClientToken asys,
1889     slv_parameters_t *parameters)
1890     {
1891     slv5_system_t sys;
1892     union parm_arg lo,hi,val;
1893     struct slv_parameter *new_parms = NULL;
1894     int32 make_macros = 0;
1895    
1896     static char *factor_names[] = {
1897     "SPK1/RANKI","SPK1/RANKI+ROW",
1898     "Fast-SPK1/RANKI","Fast-SPK1/RANKI+ROW",
1899     "Fastest-SPK1/MR-RANKI","CondQR","CPQR"
1900     /* ,"GAUSS","GAUSS_EASY" currently only works for ken */
1901     };
1902     static char *reorder_names[] = {
1903     "SPK1","TEAR_DROP","OVER_TEAR"
1904     };
1905     static char *converge_names[] = {
1906     "ABSOLUTE","RELNOM_SCALE"
1907     };
1908     static char *scaling_names[] = {
1909     "NONE","ROW_2NORM"
1910     };
1911    
1912     if (server != NULL && asys != NULL) {
1913     sys = SLV5(asys);
1914     make_macros = 1;
1915     }
1916    
1917     #ifndef NDEBUG /* keep purify from whining on UMR */
1918     lo.argr = hi.argr = val.argr = 0.0;
1919     #endif
1920    
1921     if (parameters->parms == NULL) {
1922     /* an external client wants our parameter list.
1923     * an instance of slv5_system_structure has this pointer
1924     * already set in slv5_create
1925     */
1926     new_parms = (struct slv_parameter *)
1927     ascmalloc((slv5_PA_SIZE)*sizeof(struct slv_parameter));
1928     if (new_parms == NULL) {
1929     return -1;
1930     }
1931     parameters->parms = new_parms;
1932     parameters->dynamic_parms = 1;
1933     }
1934     parameters->num_parms = 0;
1935    
1936     /* begin defining parameters */
1937    
1938     slv_define_parm(parameters, bool_parm,
1939     "ignorebounds","ignore bounds?","ignore bounds?",
1940     U_p_bool(val,0),U_p_bool(lo,0),U_p_bool(hi,1),-1);
1941     SLV_BPARM_MACRO(IGNORE_BOUNDS_PTR,parameters);
1942    
1943     slv_define_parm(parameters, real_parm,
1944     "rho", "search parameter", "search parameter",
1945     U_p_real(val,0.5),U_p_real(lo, 0),U_p_real(hi,1.0), 1);
1946     SLV_RPARM_MACRO(RHO_PTR,parameters);
1947    
1948     slv_define_parm(parameters, bool_parm,
1949     "partition", "partitioning enabled", "partitioning enabled",
1950     U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1951     SLV_BPARM_MACRO(PARTITION_PTR,parameters);
1952    
1953     slv_define_parm(parameters, bool_parm,
1954     "showlessimportant", "detailed solving info",
1955     "detailed solving info",
1956     U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1957     SLV_BPARM_MACRO(SHOW_LESS_IMPT_PTR,parameters);
1958    
1959     slv_define_parm(parameters, int_parm,
1960     "timelimit", "time limit (CPU sec/block)",
1961     "time limit (CPU sec/block)",
1962     U_p_int(val,1500),U_p_int(lo, 1),U_p_int(hi,20000),1);
1963     SLV_IPARM_MACRO(TIME_LIMIT_PTR,parameters);
1964    
1965     slv_define_parm(parameters, int_parm,
1966     "iterationlimit", "max iterations/block",
1967     "max iterations/block",
1968     U_p_int(val, 30),U_p_int(lo, 1),U_p_int(hi,20000),1);
1969     SLV_IPARM_MACRO(ITER_LIMIT_PTR,parameters);
1970    
1971     slv_define_parm(parameters,real_parm,
1972     "alpha","search coefficient" ,"search coefficient" ,
1973     U_p_real(val,0.5),U_p_real(lo,0),U_p_real(hi,1.0), 1);
1974     SLV_RPARM_MACRO(ALPHA_PTR,parameters);
1975    
1976     slv_define_parm(parameters, real_parm,
1977     "singtol", "epsilon (min pivot)", "epsilon (min pivot)",
1978     U_p_real(val, 1e-12),U_p_real(lo, 1e-12),U_p_real(hi,1.0),1);
1979     SLV_RPARM_MACRO(SING_TOL_PTR,parameters);
1980    
1981     slv_define_parm(parameters, real_parm,
1982     "pivottol", "condition tolerance", "condition tolerance",
1983     U_p_real(val, 0.5),U_p_real(lo, 0),U_p_real(hi, 1),1);
1984     SLV_RPARM_MACRO(PIVOT_TOL_PTR,parameters);
1985    
1986     slv_define_parm(parameters, real_parm,
1987     "feastol", "max residual (absolute)",
1988     "max residual (absolute)",
1989     U_p_real(val, 1e-8),U_p_real(lo, 1e-13),
1990     U_p_real(hi,100000.5),1);
1991     SLV_RPARM_MACRO(FEAS_TOL_PTR,parameters);
1992    
1993     slv_define_parm(parameters, bool_parm,
1994     "lifds", "show singletons details", IEX(0),
1995     U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
1996     SLV_BPARM_MACRO(LIFDS_PTR,parameters);
1997    
1998     slv_define_parm(parameters, bool_parm,
1999     "savlin", "write to file SlvLinsol.dat", IEX(1),
2000     U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2001     SLV_BPARM_MACRO(SAVLIN_PTR,parameters);
2002    
2003     slv_define_parm(parameters, bool_parm,
2004     "safe_calc", "safe calculations", IEX(8),
2005     U_p_bool(val, 1),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2006    
2007     SLV_BPARM_MACRO(SAFE_CALC_PTR,parameters);
2008    
2009    
2010     slv_define_parm(parameters, int_parm,
2011     "cutoff", "block size cutoff (MODEL-based)", IEX(2),
2012     U_p_int(val, 500),U_p_int(lo,0),U_p_int(hi,20000), 2);
2013     SLV_IPARM_MACRO(CUTOFF_PTR,parameters);
2014    
2015    
2016     slv_define_parm(parameters, int_parm,
2017     "upjac", "Jacobian update frequency", IEX(3),
2018     U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
2019     SLV_IPARM_MACRO(UPDATE_JACOBIAN_PTR,parameters);
2020    
2021     slv_define_parm(parameters, int_parm,
2022     "upwts", "Row scaling update frequency", IEX(4),
2023     U_p_int(val, 1),U_p_int(lo,0),U_p_int(hi,20000), 3);
2024     SLV_IPARM_MACRO(UPDATE_WEIGHTS_PTR,parameters);
2025    
2026     slv_define_parm(parameters, int_parm,
2027     "upnom", "Column scaling update frequency", IEX(5),
2028     U_p_int(val, 1000),U_p_int(lo,0),U_p_int(hi,20000), 3);
2029     SLV_IPARM_MACRO(UPDATE_NOMINALS_PTR,parameters);
2030    
2031     slv_define_parm(parameters, char_parm,
2032     "convopt", "convergence test", "convergence test",
2033     U_p_string(val,converge_names[0]),
2034     U_p_strings(lo,converge_names),
2035     U_p_int(hi,sizeof(converge_names)/sizeof(char *)),1);
2036     SLV_CPARM_MACRO(CONVOPT_PTR,parameters);
2037    
2038     slv_define_parm(parameters, char_parm,
2039     "scaleopt", "jacobian scaling option", IEX(9),
2040     U_p_string(val,scaling_names[1]),
2041     U_p_strings(lo,scaling_names),
2042     U_p_int(hi,sizeof(scaling_names)/sizeof(char *)),1);
2043     SLV_CPARM_MACRO(SCALEOPT_PTR,parameters);
2044    
2045     slv_define_parm(parameters, bool_parm,
2046     "cncols", "Check poorly scaled columns", IEX(6),
2047     U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2048     SLV_BPARM_MACRO(DUMPCNORM_PTR,parameters);
2049    
2050     slv_define_parm(parameters, bool_parm,
2051     "lintime", "Enable linsolqr timing",
2052     "Enable linsolqr factor timing",
2053     U_p_bool(val, 0),U_p_bool(lo,0),U_p_bool(hi,1), 2);
2054     SLV_BPARM_MACRO(LINTIME_PTR,parameters);
2055    
2056    
2057     slv_define_parm(parameters, char_parm,
2058     "reorder", "reorder method", IEX(7),
2059     U_p_string(val,reorder_names[0]),
2060     U_p_strings(lo,reorder_names),
2061     U_p_int(hi,sizeof(reorder_names)/sizeof(char *)),1);
2062     SLV_CPARM_MACRO(REORDER_OPTION_PTR,parameters);
2063    
2064    
2065     slv_define_parm(parameters, real_parm,
2066     "toosmall", "default for zero nominal", REX(0),
2067     U_p_real(val, 1e-8),U_p_real(lo, 1e-12),U_p_real(hi,1.5), 3);
2068     SLV_RPARM_MACRO(TOO_SMALL_PTR,parameters);
2069    
2070     slv_define_parm(parameters, real_parm,
2071     "cnlow", "smallest allowable column norm", REX(1),
2072     U_p_real(val, 0.01),U_p_real(lo, 0),
2073     U_p_real(hi,100000000.5), 3);
2074     SLV_RPARM_MACRO(CNLOW_PTR,parameters);
2075    
2076     slv_define_parm(parameters, real_parm,
2077     "cnhigh", "largest allowable column norm", REX(2),
2078     U_p_real(val, 100.0),U_p_real(lo,0),
2079     U_p_real(hi,100000000.5), 3);
2080     SLV_RPARM_MACRO(CNHIGH_PTR,parameters);
2081    
2082     slv_define_parm(parameters, real_parm,
2083     "tobnds", "fraction move to bounds", REX(3),
2084     U_p_real(val, 0.995),U_p_real(lo, 0),U_p_real(hi,1.0), 3);
2085     SLV_RPARM_MACRO(TOWARD_BOUNDS_PTR,parameters);
2086    
2087     slv_define_parm(parameters, char_parm,
2088     "bppivoting","linear method","linear method choice",
2089     U_p_string(val,factor_names[4]),
2090     U_p_strings(lo,factor_names),
2091     U_p_int(hi,sizeof(factor_names)/sizeof(char *)),1);
2092     SLV_CPARM_MACRO(FACTOR_OPTION_PTR,parameters);
2093    
2094     return 1;
2095     }
2096    
2097     /*
2098     * External routines
2099     * -----------------
2100     * See slv_client.h
2101     */
2102    
2103     static SlvClientToken slv5_create(slv_system_t server, int *statusindex)
2104     {
2105     slv5_system_t sys;
2106    
2107     sys = (slv5_system_t)asccalloc(1, sizeof(struct slv5_system_structure) );
2108     if (sys==NULL) {
2109     *statusindex = 1;
2110     return sys;
2111     }
2112     SERVER = server;
2113     sys->p.parms = sys->pa;
2114     sys->p.dynamic_parms = 0;
2115     slv5_get_default_parameters(server,(SlvClientToken)sys,&(sys->p));
2116     sys->integrity = OK;
2117     sys->presolved = 0;
2118     sys->p.output.more_important = stdout;
2119     sys->p.output.less_important = stdout;
2120     sys->J.old_partition = TRUE;
2121     sys->p.whose = (*statusindex);
2122     sys->s.ok = TRUE;
2123     sys->s.calc_ok = TRUE;
2124     sys->s.costsize = 0;
2125     sys->s.cost = NULL; /*redundant, but sanity preserving */
2126     sys->vlist = slv_get_solvers_var_list(server);
2127     sys->rlist = slv_get_solvers_rel_list(server);
2128     sys->obj = slv_get_obj_relation(server);
2129     if (sys->vlist == NULL) {
2130     ascfree(sys);
2131     FPRINTF(ASCERR,"IPSlv called with no variables.\n");
2132     *statusindex = -2;
2133     return NULL; /*_prolly leak here */
2134     }
2135     if (sys->rlist == NULL && sys->obj == NULL) {
2136     ascfree(sys);
2137     FPRINTF(ASCERR,"IPSlv called with no relations or objective.\n");
2138     *statusindex = -1;
2139     return NULL; /*_prolly leak here */
2140     }
2141     /* we don't give a damn about the objective list or the pars or
2142     * bounds or extrels or any of the other crap.
2143     */
2144     slv_check_var_initialization(server);
2145     *statusindex = 0;
2146     return((SlvClientToken)sys);
2147    
2148     }
2149    
2150     static void destroy_matrices( slv5_system_t sys)
2151     {
2152     if( sys->J.sys ) {
2153     int count = linsolqr_number_of_rhs(sys->J.sys)-1;
2154     for( ; count >= 0; count-- ) {
2155     destroy_array(linsolqr_get_rhs(sys->J.sys,count));
2156     }
2157     mtx_destroy(linsolqr_get_matrix(sys->J.sys));
2158     linsolqr_set_matrix(sys->J.sys,NULL);
2159     linsolqr_destroy(sys->J.sys);
2160     sys->J.sys = NULL;
2161     }
2162     }
2163    
2164     static void destroy_vectors( slv5_system_t sys)
2165     {
2166     destroy_array(sys->nominals.vec);
2167     destroy_array(sys->weights.vec);
2168     destroy_array(sys->variables.vec);
2169     destroy_array(sys->residuals.vec);
2170     destroy_array(sys->newton_residuals.vec);
2171     destroy_array(sys->perturbed_residuals.vec);
2172     destroy_array(sys->newton.vec);
2173     destroy_array(sys->perturbed_newton.vec);
2174     destroy_array(sys->varnewstep.vec);
2175     destroy_array(sys->varstep.vec);
2176     }
2177    
2178     static int slv5_eligible_solver(slv_system_t server)
2179     {
2180     struct rel_relation **rp;
2181     rel_filter_t rfilter;
2182    
2183     rfilter.matchbits = (REL_INCLUDED | REL_ACTIVE);
2184     rfilter.matchvalue = (REL_INCLUDED | REL_ACTIVE);
2185    
2186     if (!slv_count_solvers_rels(server,&rfilter)) {
2187     return (FALSE);
2188     }
2189    
2190     for( rp=slv_get_solvers_rel_list(server); *rp != NULL ; ++rp ) {
2191     if( rel_less(*rp) || rel_greater(*rp) ) return(FALSE);
2192     }
2193     return(TRUE);
2194     }
2195    
2196     static
2197     void slv5_get_parameters(slv_system_t server, SlvClientToken asys,
2198     slv_parameters_t *parameters)
2199     {
2200     slv5_system_t sys;
2201     (void) server;
2202     sys = SLV5(asys);
2203     if (check_system(sys)) return;
2204     mem_copy_cast(&(sys->p),parameters,sizeof(slv_parameters_t));
2205     }
2206    
2207     static void slv5_set_parameters(slv_system_t server, SlvClientToken asys,
2208     slv_parameters_t *parameters)
2209     {
2210     slv5_system_t sys;
2211     (void) server;
2212     sys = SLV5(asys);
2213     if (check_system(sys)) return;
2214     mem_copy_cast(parameters,&(sys->p),sizeof(slv_parameters_t));
2215     }
2216    
2217     static void slv5_get_status(slv_system_t server, SlvClientToken asys,
2218     slv_status_t *status)
2219     {
2220     slv5_system_t sys;
2221     (void) server;
2222     sys = SLV5(asys);
2223     if (check_system(sys)) return;
2224     mem_copy_cast(&(sys->s),status,sizeof(slv_status_t));
2225     }
2226    
2227     static linsolqr_system_t slv5_get_linsolqr_sys(slv_system_t server,
2228     SlvClientToken asys)
2229     {
2230     slv5_system_t sys;
2231     (void) server;
2232     sys = SLV5(asys);
2233     if (check_system(sys)) return NULL;
2234     return(sys->J.sys);
2235     }
2236    
2237     static linsol_system_t slv5_get_linsol_sys(slv_system_t server,
2238     SlvClientToken asys)
2239     {
2240     (void) server;
2241     (void) asys;
2242     return( NULL );
2243     }
2244    
2245     /*
2246     * Performs structural analysis on the system, setting the flags in
2247     * status. The problem must be set up, the relation/variable list
2248     * must be non-NULL. The
2249     * jacobian (linear) system must be created and have the correct order
2250     * (stored in sys->cap). Everything else will be determined here.
2251     * On entry there isn't yet a correspondence between var_sindex and
2252     * jacobian column. Here we establish that.
2253     */
2254     static void structural_analysis(slv_system_t server, slv5_system_t sys)
2255     {
2256     var_filter_t vfilter;
2257     rel_filter_t rfilter;
2258    
2259     /*
2260     * The server has marked incidence flags already.
2261     */
2262     /* count included equalities */
2263     rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
2264     rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
2265     sys->rused = slv_count_solvers_rels(server,&rfilter);
2266    
2267     /* count free and incident vars */
2268     vfilter.matchbits = (VAR_FIXED | VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
2269     vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
2270     sys->vused = slv_count_solvers_vars(server,&vfilter);
2271    
2272     /* Symbolic analysis */
2273     sys->rtot = slv_get_num_solvers_rels(server);
2274     sys->vtot = slv_get_num_solvers_vars(server);
2275     if (sys->rtot) {
2276     slv_block_partition(server);
2277     }
2278     sys->J.dofdata = slv_get_dofdata(server);
2279     sys->rank = sys->J.dofdata->structural_rank;
2280    
2281     if( !(PARTITION) ) {
2282     /* maybe we should reorder blocks here? maybe not */
2283     slv_block_unify(server);
2284     }
2285    
2286     slv_check_bounds(SERVER,sys->vused,sys->vtot-1,MIF(sys),"fixed");
2287    
2288     /* Initialize Status */
2289     sys->s.over_defined = (sys->rused > sys->vused);
2290     sys->s.under_defined = (sys->rused < sys->vused);
2291     sys->s.struct_singular = (sys->rank < sys->rused);
2292     sys->s.block.number_of = (slv_get_solvers_blocks(SERVER))->nblocks;
2293     }
2294    
2295    
2296    
2297     static void set_factor_options (slv5_system_t sys)
2298     {
2299     if (strcmp(FACTOR_OPTION,"SPK1/RANKI") == 0) {
2300     sys->J.fm = ranki_kw;
2301     } else if (strcmp(FACTOR_OPTION,"SPK1/RANKI+ROW") == 0) {
2302     sys->J.fm = ranki_jz;
2303     } else if (strcmp(FACTOR_OPTION,"Fast-SPK1/RANKI") == 0) {
2304     sys->J.fm = ranki_kw2;
2305     } else if (strcmp(FACTOR_OPTION,"Fast-SPK1/RANKI+ROW") == 0) {
2306     sys->J.fm = ranki_jz2;
2307     } else if (strcmp(FACTOR_OPTION,"Fastest-SPK1/MR-RANKI") == 0) {
2308     sys->J.fm = ranki_ba2;
2309     } else {
2310     sys->J.fm = ranki_ba2;
2311     }
2312     mtx_set_order(sys->J.mtx,sys->cap);
2313    
2314     linsolqr_set_matrix(sys->J.sys,sys->J.mtx);
2315     linsolqr_prep(sys->J.sys,linsolqr_fmethod_to_fclass(sys->J.fm));
2316     linsolqr_set_pivot_zero(sys->J.sys, SING_TOL);
2317    
2318     /* KHACK NOTE: looks like drop tol never set on interface */
2319     linsolqr_set_drop_tolerance(sys->J.sys, sys->p.tolerance.drop);
2320     linsolqr_set_pivot_tolerance(sys->J.sys, PIVOT_TOL);
2321     /* this next one is fishy, but we don't use qr so not panicking */
2322     linsolqr_set_condition_tolerance(sys->J.sys, PIVOT_TOL);
2323     }
2324    
2325    
2326     /*
2327     configures linsolqr system, among other things.
2328     sets type to be ranki_kw or ranki_jz as determined by
2329     parameters.
2330     */
2331     static void create_matrices(slv_system_t server, slv5_system_t sys)
2332     {
2333     sys->J.sys = linsolqr_create();
2334     sys->J.mtx = mtx_create();
2335    
2336     set_factor_options(sys);
2337    
2338     /* rhs 0 for sys->newton */
2339     sys->J.rhs = create_array(sys->cap,real64);
2340     linsolqr_add_rhs(sys->J.sys,sys->J.rhs,FALSE);
2341     /* rhs 1 for sys->perturbed_newton */
2342     sys->J.rhs = create_array(sys->cap,real64);
2343     linsolqr_add_rhs(sys->J.sys,sys->J.rhs,FALSE);
2344     structural_analysis(server,sys);
2345    
2346     }
2347    
2348     static void create_vectors(sys)
2349     slv5_system_t sys;
2350     {
2351     sys->nominals.vec = create_array(sys->cap,real64);
2352     sys->nominals.rng = &(sys->J.reg.col);
2353     sys->weights.vec = create_array(sys->cap,real64);
2354     sys->weights.rng = &(sys->J.reg.row);
2355     sys->variables.vec = create_array(sys->cap,real64);
2356     sys->variables.rng = &(sys->J.reg.col);
2357     sys->residuals.vec = create_array(sys->cap,real64);
2358     sys->residuals.rng = &(sys->J.reg.row);
2359     sys->newton_residuals.vec = create_array(sys->cap,real64);
2360     sys->newton_residuals.rng = &(sys->J.reg.row);
2361     sys->perturbed_residuals.vec = create_array(sys->cap,real64);
2362     sys->perturbed_residuals.rng = &(sys->J.reg.row);
2363     sys->newton.vec = create_array(sys->cap,real64);
2364     sys->newton.rng = &(sys->J.reg.col);
2365     sys->perturbed_newton.vec = create_array(sys->cap,real64);
2366     sys->perturbed_newton.rng = &(sys->J.reg.col);
2367     sys->varnewstep.vec = create_array(sys->cap,real64);
2368     sys->varnewstep.rng = &(sys->J.reg.col);
2369     sys->varstep.vec = create_array(sys->cap,real64);
2370     sys->varstep.rng = &(sys->J.reg.col);
2371    
2372     }
2373    
2374     static
2375     void slv5_dump_internals(slv_system_t server, SlvClientToken sys,int level)
2376     {
2377     (void) server;
2378     check_system(sys);
2379     if (level > 0) {
2380     FPRINTF(ASCERR,"ERROR: (slv5) slv5_dump_internals\n");
2381     FPRINTF(ASCERR," slv5 does not dump its internals.\n");
2382     }
2383     }
2384    
2385     /*
2386     * Here we will check if any fixed or included flags have
2387     * changed since the last presolve.
2388     */
2389     static
2390     int32 slv5_dof_changed(slv5_system_t sys)
2391     {
2392     int32 ind, result = 0;
2393     /* Currently we have two copies of the fixed and included flags
2394     which must be kept in sync. The var_fixed and rel_included
2395     functions perform the syncronization and hence must be called
2396     over the whole var list and rel list respectively. When we move
2397     to using only one set of flags (bit flags) this function can
2398     be changed to return 1 at the first indication of a change
2399     in the dof. */
2400    
2401     /* search for vars that were fixed and are now free */
2402     for( ind = sys->vused; ind < sys->vtot; ++ind ) {
2403     if( !var_fixed(sys->vlist[ind]) && var_active(sys->vlist[ind]) ) {
2404     ++result;
2405     }
2406     }
2407     /* search for rels that were unincluded and are now included */
2408     for( ind = sys->rused; ind < sys->rtot; ++ind ) {
2409     if( rel_included(sys->rlist[ind]) && rel_active(sys->rlist[ind])) {
2410     ++result;
2411     }
2412     }
2413     /* search for vars that were free and are now fixed */
2414     for( ind = sys->vused -