/[ascend]/trunk/solvers/ipslv/slv5.mehrotra_no_quad.c
ViewVC logotype

Annotation of /trunk/solvers/ipslv/slv5.mehrotra_no_quad.c

Parent Directory Parent Directory | Revision Log Revision Log


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