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