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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 48803 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 /*
2     * Solver Standard Clients
3     * by Benjamin Andrew Allan
4     * 5/19/96
5     * Version: $Revision: 1.28 $
6     * Version control file: $RCSfile: slv_stdcalls.c,v $
7     * Date last modified: $Date: 1998/06/16 16:53:04 $
8     * Last modified by: $Author: mthomas $
9     *
10     * Copyright(C) 1996 Benjamin Andrew Allan
11     * Copyright(C) 1998 Carnegie Mellon University
12     *
13     * This file is part of the ASCEND IV math programming system.
14     * Here is where we register our normal solvers and keep things
15     * like reordering clients and so forth. This file should probably
16     * be broken up into one file-one client layout for the unregistered
17     * clients so that we can load them selectively.
18     *
19     * The Ascend Math Programming System is free software; you can
20     * redistribute it and/or modify it under the terms of the GNU
21     * General Public License as published by the Free Software
22     * Foundation; either version 2 of the License, or (at your option)
23     * any later version.
24     *
25     * The Ascend Math Programming System is distributed in hope that it
26     * will be useful, but WITHOUT ANY WARRANTY; without even the implied
27     * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
28     * See the GNU General Public License for more details.
29     *
30     * You should have received a copy of the GNU General Public License
31     * along with the program; if not, write to the Free Software
32     * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check
33     * the file named COPYING.
34     */
35    
36     #include "utilities/ascConfig.h"
37     #include "compiler/compiler.h"
38     #include "utilities/ascMalloc.h"
39     #include "general/list.h"
40     #include "solver/mtx.h"
41     #include "solver/slv_types.h"
42     #include "solver/var.h"
43     #include "solver/rel.h"
44     #include "solver/discrete.h"
45     #include "solver/conditional.h"
46     #include "solver/logrel.h"
47     #include "solver/relman.h"
48     #include "solver/logrelman.h"
49     #include "solver/bnd.h"
50     #include "solver/slv_common.h"
51     #include "solver/linsol.h"
52     #include "solver/linsolqr.h"
53     #include "solver/slv_client.h"
54     /* header of registered clients */
55     #include "solver/slv0.h"
56     #include "solver/slv1.h"
57     #include "solver/slv2.h"
58     #include "solver/slv3.h"
59     #include "solver/slv6.h"
60     #include "solver/slv7.h"
61     #include "solver/slv8.h"
62     #include "solver/slv9.h"
63     #include "solver/slv9a.h"
64     #include "solver/model_reorder.h"
65     #include "solver/slv_stdcalls.h"
66    
67     #define KILL 0 /* deleteme code. compile old code if kill = 1 */
68     #define NEEDSTOBEDONE 0 /* indicates code/comments that are not yet ready */
69     #define USECODE 0 /* Code in good shape but not used currently */
70    
71     /*
72     * Here we play some hefty jacobian reordering games.
73     * What we want to happen eventually is as follows:
74     *
75     * Get the free & incident pattern for include relations.
76     * Output assign jacobian.
77     * BLT permute the jacobian. if underspecified, fake rows
78     * to make things appear square.
79     * For all leading square blocks apply reordering (kirk, etc)
80     * If trailing rectangular block, apply clever kirkbased scheme not known.
81     *
82     * At present, we aren't quite so clever. We:
83     *
84     * Get the free & incident pattern for include relations.
85     * Output assign jacobian.
86     * BLT permute the square output assigned region of the jacobian.
87     * For all leading square blocks apply reordering (kirk, etc)
88     *
89     * Collect the block list as part of the master data structure.
90     * Set sindices for rels, vars as current rows/cols in matrix so
91     * that the jacobian is 'naturally' preordered for all solvers.
92     *
93     * Solvers are still free to reorder their own matrices any way they like;
94     * it's probably a dumb idea, though.
95     *
96     * returns 0 if ok, 1 if out of memory.
97     */
98    
99     #define MIMDEBUG 0
100     /*
101     * Populates a matrix according to the sys solvers_vars, solvers_rels
102     * lists and the filters given. The filter should have at least
103     * SVAR = 1 bit on. mtx given must be created (not null) and with
104     * order >= MAX( slv_get_num_solvers_rels(sys),slv_get_num_solvers_vars(sys));
105     * returns 0 if went ok.
106     */
107     int slv_std_make_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
108     var_filter_t *vf,rel_filter_t *rf)
109     {
110     #if MIMDEBUG
111     FILE *fp;
112     #endif
113     int32 r, rlen,ord;
114     struct rel_relation **rp;
115    
116     if (sys==NULL || mtx == NULL || vf == NULL || rf == NULL) {
117     FPRINTF(stderr,"make_incidence called with null\n");
118     return 1;
119     }
120     rp = slv_get_solvers_rel_list(sys);
121     assert(rp!=NULL);
122     rlen = slv_get_num_solvers_rels(sys);
123     ord = MAX(rlen,slv_get_num_solvers_vars(sys));
124     if (ord > mtx_order(mtx)) {
125     FPRINTF(stderr,"make_incidence called with undersized matrix\n");
126     return 2;
127     }
128     for (r=0; r < rlen; r++) {
129     if (rel_apply_filter(rp[r],rf)) {
130     relman_get_incidence(rp[r],vf,mtx);
131     }
132     }
133     #if MIMDEBUG
134     fp = fopen("/tmp/mim.plot","w+");
135     if (fp !=NULL) {
136     mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
137     fclose(fp);
138     }
139     #endif
140     return 0;
141     }
142    
143     #define RIDEBUG 0
144     /* returns 0 if successful, 1 if insufficient memory. Does not change
145     * the data in mtx. Orders the solvers_var list of the system to
146     * match the permutation on the given mtx. It is assumed that
147     * the org cols of mtx == var list position. Only the
148     * vars in range lo to hi of the var list are affected and
149     * only these vars should appear in the same org column range of
150     * the mtx. This should not be called on blocks less than 3x3.
151     */
152     static int reindex_vars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
153     const mtx_matrix_t mtx)
154     {
155     struct var_variable **vtmp, **vp;
156     int32 c,v,vlen;
157    
158     if (lo >= hi +1) return 0; /* job too small */
159     vp = slv_get_solvers_var_list(sys);
160     vlen = slv_get_num_solvers_vars(sys);
161     /* on vtmp we DONT have the terminating null */
162     vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
163     if (vtmp == NULL) {
164     return 1;
165     }
166     /* copy pointers to vtmp in order desired and copy back rather than sort */
167     for (c=lo;c<=hi;c++) {
168     v = mtx_col_to_org(mtx,c);
169     vtmp[c] = vp[v];
170     }
171     /* copying back and re-sindexing */
172     for (c=lo;c<=hi;c++) {
173     vp[c] = vtmp[c];
174     var_set_sindex(vp[c],c);
175     }
176     ascfree(vtmp);
177     return 0;
178     }
179     /* returns 0 if successful, 1 if insufficient memory. Does not change
180     * the data in mtx. Orders the solvers_rel list of the system to
181     * match the permutation on the given mtx. It is assumed that
182     * the org rows of mtx == rel list position. Only the
183     * rels in range lo to hi of the rel list are affected and
184     * only these rels should appear in the same org row range of
185     * the input mtx. This should not be called on blocks less than 3x3.
186     */
187     static int reindex_rels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
188     const mtx_matrix_t mtx)
189     {
190     struct rel_relation **rtmp, **rp;
191     int32 c,v,rlen;
192    
193     rp = slv_get_solvers_rel_list(sys);
194     rlen = slv_get_num_solvers_rels(sys);
195     /* on rtmp we DONT have the terminating null */
196     rtmp = (struct rel_relation **)ascmalloc(rlen*sizeof(struct rel_relation *));
197     if (rtmp == NULL) {
198     return 1;
199     }
200     /* copy pointers to rtmp in order desired and copy back rather than sort.
201     * do this only in the row range of interest. */
202     for (c=lo;c<=hi;c++) {
203     v = mtx_row_to_org(mtx,c);
204     #if RIDEBUG
205     FPRINTF(stderr,"Old rel sindex (org) %d becoming sindex (cur) %d\n",v,c);
206     #endif
207     rtmp[c] = rp[v];
208     }
209     /* copying back and re-sindexing */
210     for (c=lo;c<=hi;c++) {
211     rp[c] = rtmp[c];
212     rel_set_sindex(rp[c],c);
213     }
214     ascfree(rtmp);
215     return 0;
216     }
217    
218    
219     /* returns 0 if ok, OTHERWISE if madness detected.
220     * doesn't grok inequalities.
221     */
222     #define SBPDEBUG 0
223     int slv_block_partition_real(slv_system_t sys,int uppertriangular)
224     {
225     #if SBPDEBUG
226     FILE *fp;
227     #endif
228     struct rel_relation **rp;
229     struct var_variable **vp;
230     mtx_region_t *newblocks;
231     dof_t *d;
232     int32 nrow,ncol;
233     int32 ncolpfix, nrowpun, vlenmnv;
234     mtx_matrix_t mtx;
235     int32 order, rank;
236     int32 c,len,vlen,r,rlen;
237     var_filter_t vf;
238     rel_filter_t rf;
239    
240     rp = slv_get_solvers_rel_list(sys);
241     vp = slv_get_solvers_var_list(sys);
242     rlen = slv_get_num_solvers_rels(sys);
243     vlen = slv_get_num_solvers_vars(sys);
244     if (rlen ==0 || vlen == 0) return 1;
245     order = MAX(rlen,vlen);
246    
247     rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
248     rf.matchvalue = (REL_ACTIVE);
249     vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
250     vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
251    
252     /* nrow plus unincluded rels and ineqs */
253     nrowpun = slv_count_solvers_rels(sys,&rf);
254     ncolpfix = slv_count_solvers_vars(sys,&vf); /* ncol plus fixed vars */
255    
256     vf.matchbits = (VAR_SVAR);
257     vf.matchvalue = 0;
258     vlenmnv = vlen - slv_count_solvers_vars(sys,&vf); /* vlen minus non vars */
259    
260     rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
261     rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
262     vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
263     vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
264    
265     nrow = slv_count_solvers_rels(sys,&rf);
266     ncol = slv_count_solvers_vars(sys,&vf);
267    
268     mtx = mtx_create();
269     mtx_set_order(mtx,order);
270    
271     if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
272     FPRINTF(stderr,
273     "slv_block_partition: failure in creating incidence matrix.\n");
274     mtx_destroy(mtx);
275     return 1;
276     }
277    
278     mtx_output_assign(mtx,rlen,vlen);
279     rank = mtx_symbolic_rank(mtx);
280     if (rank == 0 ) return 1; /* nothing to do, eh? */
281     /* lot of whining about dof */
282     if (rank < nrow) {
283     FPRINTF(stderr,"System is row rank deficient (%d dependent equations)\n",
284     nrow - rank);
285     }
286     if (rank < ncol) {
287     if ( nrow != rank) {
288     FPRINTF(stderr,"System is row rank deficient with %d excess columns.\n",
289     ncol - rank);
290     } else {
291     FPRINTF(stderr,"System has %d degrees of freedom.\n", ncol - rank);
292     }
293     }
294     if (ncol == nrow) {
295     FPRINTF(stderr,"System is (%d) square",ncol);
296     if (ncol != rank) {
297     FPRINTF(stderr,"but rank deficient.\n");
298     } else {
299     FPRINTF(stderr,".\n");
300     }
301     }
302     if (uppertriangular) {
303     mtx_ut_partition(mtx);
304     } else {
305     mtx_partition(mtx);
306     }
307     /* copy the block list. there is at least one if rank >=1 as above */
308     len = mtx_number_of_blocks(mtx);
309     newblocks = (mtx_region_t *)ascmalloc(len*sizeof(mtx_region_t));
310     if (newblocks == NULL) {
311     mtx_destroy(mtx);
312     return 2;
313     }
314     for (c = 0 ; c < len; c++) {
315     mtx_block(mtx,c,&(newblocks[c]));
316     }
317     slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
318     d = slv_get_dofdata(sys);
319     d->structural_rank = rank;
320     d->n_rows = nrow;
321     d->n_cols = ncol;
322     /* the next two lines assume inequalities are unincluded.
323     */
324     #define MINE 1
325     #if MINE
326     d->n_fixed = ncolpfix - ncol;
327     d->n_unincluded = nrowpun - nrow;
328     #else
329     d->n_fixed = vlen - ncol;
330     d->n_unincluded = rlen - nrow;
331     #endif
332     d->reorder.partition = 1; /* yes */
333     d->reorder.basis_selection = 0; /* none yet */
334     d->reorder.block_reordering = 0; /* none */
335    
336     #if SBPDEBUG
337     fp = fopen("/tmp/sbp1.plot","w+");
338     if (fp !=NULL) {
339     mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
340     fclose(fp);
341     }
342     #endif
343    
344    
345     #if MINE
346     len = vlenmnv - 1; /* row of last possible solver variable */
347     for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
348     r = mtx_col_to_org(mtx,c);
349     if ( var_fixed(vp[r]) && var_active(vp[r])) {
350     c++;
351     } else {
352     mtx_swap_cols(mtx,c,len);
353     len--;
354     }
355     }
356     #endif
357    
358     /* Now, leftover columns are where we want them.
359     * That is, unassigned incidence is next to assigned and fixed on right with
360     * no columns which do not correspond to variables in the range 0..vlen-1.
361     */
362     /* there aren't any nonincident vars in the solvers varlist, so
363     * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
364    
365     /* rows 0->rlen-1 should be actual rels, though, and not phantoms.
366     * At this point we should have
367     * rows - assigned
368     * - unassigned w/incidence
369     * - unincluded and included w/out free incidence and inactive mixed.
370     */
371     /* sort unincluded relations to bottom of good region */
372     /* this is needed because the idiot user may have fixed everything */
373     /* in a relation, but we want those rows not to be confused with */
374     /* unincluded rows. we're sort of dropping inequalities on the floor.*/
375     len = rlen - 1; /* row of last possible relation */
376     /* sort the inactive out of the unassigned/unincluded */
377     for (c=rank; len > c ; ) {
378     r = mtx_row_to_org(mtx,c);
379     if ( rel_active(rp[r])) {
380     c++;
381     } else {
382     mtx_swap_rows(mtx,c,len);
383     len--;
384     }
385     }
386    
387     #if MINE
388     /* sort the unincluded out of the unassigned */
389     len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
390     for (c=rank; len > c ; ) {
391     r = mtx_row_to_org(mtx,c);
392     if ( rel_included(rp[r]) ) {
393     c++;
394     } else {
395     mtx_swap_rows(mtx,c,len);
396     len--;
397     }
398     }
399     #endif
400    
401    
402     #if SBPDEBUG
403     fp = fopen("/tmp/sbp2.plot","w+");
404     if (fp !=NULL) {
405     mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
406     fclose(fp);
407     }
408     #endif
409     /* now, at last we have cols jacobian in the order we want the lists to
410     * be handed to the solvers. So, we need to reset the var sindex values
411     * and reorder the lists.
412     */
413     if (reindex_vars_from_mtx(sys,0,vlen-1,mtx)) {
414     mtx_destroy(mtx);
415     return 2;
416     }
417     /* now, at last we have rows jacobian in the order we want the lists to
418     * be handed to the solvers. So, we need to reset the rel sindex values.
419     */
420     if (reindex_rels_from_mtx(sys,0,rlen-1,mtx)) {
421     mtx_destroy(mtx);
422     return 2;
423     }
424     mtx_destroy(mtx);
425     return 0;
426     }
427    
428    
429     #if 0 /* code not currently used */
430     #ifdef STATIC_HARWELL
431     extern void mc21b();
432     extern void mc13emod();
433     /* returns 0 if ok, OTHERWISE if madness detected.
434     * doesn't grok inequalities.
435     * CURRENTLY DOESN'T DO WELL WHEN NCOL<NROW
436     */
437     #define SBPDEBUG 0
438     int slv_block_partition_harwell(slv_system_t sys)
439     {
440     #if SBPDEBUG
441     FILE *fp;
442     #endif
443     struct rel_relation **rp;
444     struct rel_relation **rtmp;
445     struct rel_relation *rel;
446     struct var_variable **vp;
447     struct var_variable **vtmp;
448     struct var_variable *var;
449     const struct var_variable **list;
450     mtx_region_t *newblocks;
451     dof_t *d;
452     int32 nrow,ncol,size;
453     int32 order, rank;
454     int32 c,len,vlen,r,v,rlen,col,row;
455     int32 licn,rel_tmp_end,rel_count,var_perm_count,var_tmp_end,row_perm_count;
456     int32 col_incidence,numnz,row_high,col_high,dummy_ptr,row_ptr,num;
457     int32 *dummy_array_ptr;
458     var_filter_t vf;
459     rel_filter_t rf;
460     int32 *icn, *ip, *lenr, *iperm, *iw1, *iw2, *iw3, *arp, *ipnew, *ib;
461     int32 *row_perm;
462    
463     /* get rel and var info */
464     rp = slv_get_solvers_rel_list(sys);
465     vp = slv_get_solvers_var_list(sys);
466     rlen = slv_get_num_solvers_rels(sys);
467     vlen = slv_get_num_solvers_vars(sys);
468     if (rlen ==0 || vlen == 0) return 1;
469     order = MAX(rlen,vlen);
470    
471     /* allocate temp arrays */
472     vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
473     var = (struct var_variable *)ascmalloc(sizeof(struct var_variable *));
474     if (vtmp == NULL || var ==NULL) {
475     return 1;
476     }
477     rtmp = (struct rel_relation **)ascmalloc(rlen*sizeof(struct rel_relation *));
478     rel = (struct rel_relation *)ascmalloc(sizeof(struct rel_relation *));
479     if (rtmp == NULL || rel ==NULL) {
480     return 1;
481     }
482    
483     /* set up filters */
484     rf.matchbits = (REL_INCLUDED | REL_EQUALITY);
485     rf.matchvalue = (REL_INCLUDED | REL_EQUALITY);
486     vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
487     vf.matchvalue = (VAR_INCIDENT | VAR_SVAR);
488    
489     /* count rows and cols */
490     nrow = slv_count_solvers_rels(sys,&rf);
491     ncol = slv_count_solvers_vars(sys,&vf);
492    
493    
494     /* count incidence for equality relations that are included
495     and active. Sorts out unincluded and inactive rels to
496     the end of the temporary relation list */
497     licn = 0;
498     rel_tmp_end = rlen -1;
499     for (r=0; r < rlen; r++) {
500     rel = rp[r];
501     if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
502     len = rel_n_incidences(rel);
503     list = rel_incidence_list(rel);
504     for (c=0; c < len; c++) {
505     if( var_apply_filter(list[c],&vf) ) {
506     licn++;
507     }
508     }
509     } else {
510     rtmp[rel_tmp_end] = rp[r];
511     rel_tmp_end--;
512     }
513     }
514    
515     /* sort the inactive out of the unincluded and move to end */
516     len = rlen -1;
517     for(c = rel_tmp_end + 1; len > c; ) {
518     rel = rtmp[c];
519     if (rel_active(rel)) {
520     c++;
521     } else {
522     rtmp[c] = rtmp[len];
523     rtmp[len] = rel;
524     len--;
525     }
526     }
527    
528     /* Sort variable list */
529     var_perm_count = 0;
530     var_tmp_end = vlen - 1;
531     for (c = 0; c < vlen; c++) {
532     var = vp[c];
533     if (var_apply_filter(var,&vf)) {
534     vtmp[var_perm_count] = var;
535     var_perm_count++;
536     } else {
537     vtmp[var_tmp_end] = var;
538     var_tmp_end--;
539     }
540     }
541     /* sort the inactive out of the unincluded and move to end */
542     len = vlen -1;
543     for(c = var_tmp_end + 1; len > c; ) {
544     var = vtmp[c];
545     if (var_active(var) && var_active(var)) {
546     c++;
547     } else {
548     vtmp[c] = vtmp[len];
549     vtmp[len] = var;
550     len--;
551     }
552     }
553    
554     /* reset solver indicies to make life easier */
555     for (c = 0; c < vlen; c++) {
556     vp[c] = vtmp[c];
557     var_set_sindex(vp[c],c);
558     }
559    
560     size = MAX(nrow,ncol);
561     /* Create vectors for fortran calls */
562     icn = (int32 *)ascmalloc(licn*sizeof(int32));
563     ip = (int32 *)ascmalloc(size*sizeof(int32));
564     lenr = (int32 *)ascmalloc(size*sizeof(int32));
565     iperm = (int32 *)ascmalloc(size*sizeof(int32));
566     iw1 = (int32 *)ascmalloc(size*sizeof(int32));
567     iw2 = (int32 *)ascmalloc(size*sizeof(int32));
568     iw3 = (int32 *)ascmalloc(size*sizeof(int32));
569     arp = (int32 *)ascmalloc(size*sizeof(int32));
570     ipnew = (int32 *)ascmalloc(size*sizeof(int32));
571     ib = (int32 *)ascmalloc(size*sizeof(int32));
572     row_perm = (int32 *)ascmalloc(size*sizeof(int32));
573    
574    
575     /* Fill incidence vectors and place included relations w/o
576     * incidence in temporary relation list before the unincluded
577     * and inactive relations
578     */
579     col_incidence = 0;
580     row_perm_count = 0;
581     for (r=0; r < rlen; r++) {
582     rel = rp[r];
583     if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
584     len = rel_n_incidences(rel);
585     if (len > 0) {
586     ip[row_perm_count] = col_incidence + 1; /*FORTRAN*/
587     lenr[row_perm_count] = len;
588     row_perm[row_perm_count] = r;
589     row_perm_count++;
590     list = rel_incidence_list(rel);
591     for (c=0; c < len; c++) {
592     if( var_apply_filter(list[c],&vf) ) {
593     col = var_sindex(list[c]);
594     icn[col_incidence] = col + 1; /*FORTRAN*/
595     col_incidence++;
596     }
597     }
598     } else {
599     rtmp[rel_tmp_end] = rel;
600     rel_tmp_end--;
601     }
602     }
603     }
604    
605     licn = col_incidence;
606     order = row_perm_count;
607     mc21b(&order,icn,&licn,ip,lenr,iperm,&numnz,iw1,arp,iw2,iw3);
608     if (order == numnz){
609     FPRINTF(stderr,"they are equal\n");
610     } else if(order == numnz + 1){
611     FPRINTF(stderr,"ADD ONE\n");
612     } else {
613     FPRINTF(stderr,"no relationship\n");
614     }
615    
616     if (rank == 0 ) {
617     return 1; /* nothing to do, eh? */
618     }
619     row_high = nrow - 1;
620     col_high = ncol -1;
621    
622     /* lot of whining about dof */
623     if (rank < nrow) {
624     FPRINTF(stderr,"System is row rank deficient (dependent equations)\n");
625     row_high = rank - 1;
626     /* KHACK: have found that this is executed erroneously */
627     }
628     if (rank < ncol) {
629     if ( nrow != rank) {
630     FPRINTF(stderr,"System is row rank deficient with excess columns.\n");
631     } else {
632     FPRINTF(stderr,"System has degrees of freedom.\n");
633     }
634     }
635     if (ncol == nrow) {
636     FPRINTF(stderr,"System is (%d) square",ncol);
637     if (ncol != rank) {
638     FPRINTF(stderr,"but rank deficient.\n");
639     } else {
640     FPRINTF(stderr,".\n");
641     }
642     }
643    
644     dummy_array_ptr = arp;
645     arp = lenr;
646     lenr = dummy_array_ptr;
647    
648     /* Fix row pointers for assigned relations. Put unassigned
649     onto temporary relation list (before included w/o incidence)
650     and put dummy dense rows into ipnew */
651     dummy_ptr = licn + 1; /* FORTRAN */
652     for (row = 0; row < order; row++) {
653     row_ptr = iperm[row] - 1; /* FORTRAN */
654     if (row_ptr != -1) { /* indicates unassigned row */
655     ipnew[row_ptr] = ip[row];
656     lenr[row_ptr] = arp[row];
657     } else {
658     ipnew[row_ptr] = dummy_ptr;
659     dummy_ptr += ncol;
660     lenr[row_ptr] = ncol;
661     /* unassigned before fixed*/
662     rtmp[rel_tmp_end] = rp[row_perm[row_ptr]];
663     rel_tmp_end--;
664     }
665     }
666     for (row = order; row < ncol; row++) {
667     ipnew[row] = dummy_ptr;
668     dummy_ptr += ncol;
669     lenr[row] = ncol;
670     }
671    
672     mc13emod(&size,icn,&licn,ipnew,lenr,arp,ib,&num,iw1,iw2,iw3);
673    
674    
675     /* copy the block list. there is at least one if rank >=1 as above */
676     len = num;
677     newblocks = (mtx_region_t *)ascmalloc(len*sizeof(mtx_region_t));
678     if (newblocks == NULL) {
679     return 2;
680     }
681     /* set up locations of block corners with the row and column
682     orderings which will be set soon */
683     for (c = 0 ; c < len; c++) {
684     newblocks[c].row.low = ib[c]-1;
685     newblocks[c].col.low = newblocks[c].row.low;
686     }
687     for (c = 0 ; c < len -1; c++) {
688     newblocks[c].row.high = newblocks[c+1].row.low -1;
689     newblocks[c].col.high = newblocks[c].row.high;
690     }
691     newblocks[len - 1].row.high = row_high;
692     newblocks[len - 1].col.high = col_high;
693    
694     slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
695     d = slv_get_dofdata(sys);
696     d->structural_rank = rank;
697     d->n_rows = nrow;
698     d->n_cols = ncol;
699     /* the next two lines assume only free and fixed incident solvervar
700     * appear in the solvers var list and inequalities are unincluded.
701     */
702     d->n_fixed = vlen - ncol;
703     d->n_unincluded = rlen - nrow;
704     d->reorder.partition = 1; /* yes */
705     d->reorder.basis_selection = 0; /* none yet */
706     d->reorder.block_reordering = 0; /* none */
707    
708     for (c = 0; c < ncol; c++) {
709     v = arp[c] -1; /* FORTRAN */
710     vtmp[c] = vp[v];
711     }
712     for (c = 0; c < vlen; c++) {
713     vp[c] = vtmp[c];
714     var_set_sindex(vp[c],c);
715     }
716    
717     rel_count = 0;
718     for (c = 0; c < size; c++) {
719     r = arp[c] - 1;
720     if (r < order){ /* weed out fake rows */
721     r = iperm[r] - 1;
722     if (r != -1) { /* unassigned already in rtmp */
723     rtmp[rel_count] = rp[row_perm[r]];
724     rel_count++;
725     }
726     }
727     }
728     for (c = 0; c < rlen; c++) {
729     rp[c] = rtmp[c];
730     rel_set_sindex(rp[c],c);
731     }
732     ascfree(vtmp);
733     ascfree(rtmp);
734     ascfree(icn);
735     ascfree(ip);
736     ascfree(lenr);
737     ascfree(iperm);
738     ascfree(iw1);
739     ascfree(iw2);
740     ascfree(iw3);
741     ascfree(arp);
742     ascfree(ipnew);
743     ascfree(ib);
744     return 0;
745     }
746     #endif /*STATIC_HARWELL*/
747     #endif /* 0 */
748    
749    
750     void slv_sort_rels_and_vars(slv_system_t sys,
751     int32 *rel_count,
752     int32 *var_count)
753     {
754     struct rel_relation **rp;
755     struct rel_relation **rtmp;
756     struct rel_relation *rel;
757     struct var_variable **vp;
758     struct var_variable **vtmp;
759     struct var_variable *var;
760     int32 nrow,ncol,rlen,vlen,order,rel_tmp_end;
761     int32 r,c,var_tmp_end,len;
762     var_filter_t vf;
763     rel_filter_t rf;
764    
765     /* get rel and var info */
766     rp = slv_get_solvers_rel_list(sys);
767     vp = slv_get_solvers_var_list(sys);
768     rlen = slv_get_num_solvers_rels(sys);
769     vlen = slv_get_num_solvers_vars(sys);
770     if (rlen ==0 || vlen == 0) return;
771     order = MAX(rlen,vlen);
772    
773     assert(var_count != NULL && rel_count != NULL);
774     *var_count = *rel_count = -1;
775    
776     /* allocate temp arrays */
777     vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
778     if (vtmp == NULL) {
779     return;
780     }
781     rtmp = (struct rel_relation **)ascmalloc(rlen*sizeof(struct rel_relation *));
782     if (rtmp == NULL) {
783     ascfree(vtmp);
784     return;
785     }
786    
787     /* set up filters */
788     rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
789     rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
790     vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
791     vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
792    
793     /* count rows and cols */
794     nrow = slv_count_solvers_rels(sys,&rf);
795     ncol = slv_count_solvers_vars(sys,&vf);
796    
797    
798     /* Sort out unincluded and inactive rels to
799     * the end of the temporary relation list.
800     */
801     *rel_count = 0;
802     rel_tmp_end = rlen -1;
803     for (r=0; r < rlen; r++) {
804     rel = rp[r];
805     if (rel_apply_filter(rel,&rf)) {
806     rtmp[*rel_count] = rel;
807     (*rel_count)++;
808     } else {
809     rtmp[rel_tmp_end] = rel;
810     rel_tmp_end--;
811     }
812     }
813    
814     /* sort the inactive out of the unincluded and move to end */
815     len = rlen -1;
816     for(c = rel_tmp_end + 1; len > c; ) {
817     rel = rtmp[c];
818     if (rel_active(rel)) {
819     c++;
820     } else {
821     rtmp[c] = rtmp[len];
822     rtmp[len] = rel;
823     len--;
824     }
825     }
826    
827     /* Sort variable list */
828     *var_count = 0;
829     var_tmp_end = vlen - 1;
830     for (c = 0; c < vlen; c++) {
831     var = vp[c];
832     if (var_apply_filter(var,&vf)) {
833     vtmp[*var_count] = var;
834     (*var_count)++;
835     } else {
836     vtmp[var_tmp_end] = var;
837     var_tmp_end--;
838     }
839     }
840     /* sort the inactive out of the unincluded and move to end */
841     len = vlen -1;
842     for(c = var_tmp_end + 1; len > c; ) {
843     var = vtmp[c];
844     if (var_active(var) && var_active(var)) {
845     c++;
846     } else {
847     vtmp[c] = vtmp[len];
848     vtmp[len] = var;
849     len--;
850     }
851     }
852    
853     /* reset solver indicies */
854     for (c = 0; c < vlen; c++) {
855     vp[c] = vtmp[c];
856     var_set_sindex(vp[c],c);
857     }
858    
859     for (c = 0; c < rlen; c++) {
860     rp[c] = rtmp[c];
861     rel_set_sindex(rp[c],c);
862     }
863     ascfree(vtmp);
864     ascfree(rtmp);
865     return;
866     }
867    
868    
869    
870    
871    
872     int slv_block_unify(slv_system_t sys)
873     {
874     dof_t *d;
875     const mtx_block_t *mbt;
876     mtx_region_t *newblocks;
877    
878     if (sys==NULL) return 1;
879    
880     d = slv_get_dofdata(sys);
881     mbt = slv_get_solvers_blocks(sys);
882     assert(d!=NULL && mbt!=NULL);
883     if (d->structural_rank && mbt->nblocks >1) {
884     newblocks = (mtx_region_t *)ascmalloc(sizeof(mtx_region_t));
885     if (newblocks == NULL) return 2;
886     newblocks->row.low = newblocks->col.low = 0;
887     newblocks->row.high = d->structural_rank - 1;
888     newblocks->col.high = d->n_cols - 1;
889     if ( newblocks->col.high < 0 || newblocks->row.high < 0) {
890     ascfree(newblocks);
891     return 1;
892     }
893     slv_set_solvers_blocks(sys,1,newblocks);
894     }
895     return 0;
896     }
897    
898     int slv_set_up_block(slv_system_t sys,int bnum)
899     {
900     struct rel_relation **rp;
901     struct var_variable **vp;
902     mtx_region_t reg;
903     const mtx_block_t *b;
904     int32 c,vlen,rlen;
905     var_filter_t vf;
906     rel_filter_t rf;
907    
908     if (sys==NULL) return 1;
909     rlen = slv_get_num_solvers_rels(sys);
910     vlen = slv_get_num_solvers_vars(sys);
911     if (rlen ==0 || vlen == 0) return 1;
912    
913     rp = slv_get_solvers_rel_list(sys);
914     vp = slv_get_solvers_var_list(sys);
915     assert(rp!=NULL);
916     assert(vp!=NULL);
917    
918     rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
919     rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
920     vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
921     vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
922    
923     b = slv_get_solvers_blocks(sys);
924     assert(b!=NULL);
925     if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
926     reg = b->block[bnum];
927     for (c=reg.col.low; c<=reg.col.high; c++) {
928     var_set_in_block(vp[c],1);
929     }
930     for (c=reg.row.low; c<=reg.row.high; c++) {
931     rel_set_in_block(rp[c],1);
932     }
933     return 0;
934     }
935    
936     /* returns 0 if ok, OTHERWISE if madness detected.
937     */
938     int slv_spk1_reorder_block(slv_system_t sys,int bnum,int transpose)
939     {
940     struct rel_relation **rp;
941     struct var_variable **vp;
942     mtx_region_t reg;
943     const mtx_block_t *b;
944     mtx_matrix_t mtx;
945     mtx_coord_t coord;
946     int32 c,vlen,rlen;
947     var_filter_t vf;
948     rel_filter_t rf;
949     dof_t *d;
950    
951     if (sys==NULL) return 1;
952     rlen = slv_get_num_solvers_rels(sys);
953     vlen = slv_get_num_solvers_vars(sys);
954     if (rlen ==0 || vlen == 0) return 1;
955    
956     rp = slv_get_solvers_rel_list(sys);
957     vp = slv_get_solvers_var_list(sys);
958     assert(rp!=NULL);
959     assert(vp!=NULL);
960     rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
961     rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
962     vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
963     vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
964    
965     mtx = mtx_create();
966     mtx_set_order(mtx,MAX(rlen,vlen));
967     b = slv_get_solvers_blocks(sys);
968     assert(b!=NULL);
969     if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
970     reg = b->block[bnum];
971     for (c=reg.col.low; c<=reg.col.high; c++) {
972     var_set_in_block(vp[c],1);
973     }
974     for (c=reg.row.low; c<=reg.row.high; c++) {
975     rel_set_in_block(rp[c],1);
976     }
977     if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
978     return 1; /* must be square */
979     /* could also enforce minsize 3x3, but someone my call with a
980     * partitionable region, so don't want to.
981     */
982     }
983    
984     if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
985     FPRINTF(stderr,
986     "slv_spk1_reorder_block: failure in creating incidence matrix.\n");
987     mtx_destroy(mtx);
988     return 1;
989     }
990     /* verify that block has no empty columns, though not checking diagonal */
991     for (c = reg.row.low; c <= reg.row.high; c++) {
992     coord.col = mtx_FIRST;
993     coord.row = c;
994     if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
995     mtx_destroy(mtx);
996     FPRINTF(stderr, "slv_spk1_reorder_block: empty row (%d) found.\n",c);
997     return 1;
998     }
999     coord.row = mtx_FIRST;
1000     coord.col = c;
1001     if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
1002     FPRINTF(stderr, "slv_spk1_reorder_block: empty col (%d) found.\n",c);
1003     mtx_destroy(mtx);
1004     return 1;
1005     }
1006     }
1007     if (transpose) {
1008     mtx_reorder(mtx,&reg,mtx_TSPK1);
1009     } else {
1010     mtx_reorder(mtx,&reg,mtx_SPK1);
1011     }
1012     if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
1013     mtx_destroy(mtx);
1014     return 2;
1015     }
1016     if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
1017     mtx_destroy(mtx);
1018     return 2;
1019     }
1020     d = slv_get_dofdata(sys);
1021     d->reorder.block_reordering = 1; /* spk1 */
1022    
1023     mtx_destroy(mtx);
1024     return 0;
1025     }
1026    
1027     /* global to get around the mr header */
1028     static
1029     enum mtx_reorder_method g_blockmethod = mtx_UNKNOWN;
1030    
1031     /*
1032     A function to be called on the leftover diagonal blocks of mr_bisect.
1033     This function should be thoroughly investigated, which it has not been.
1034     */
1035     static int tear_subreorder(slv_system_t server,
1036     mtx_matrix_t mtx, mtx_region_t *reg)
1037     {
1038     assert(mtx != NULL && reg !=NULL);
1039     (void)server;
1040     if (g_blockmethod==mtx_UNKNOWN) {
1041     mtx_reorder(mtx,reg,mtx_SPK1);
1042     } else {
1043     mtx_reorder(mtx,reg,g_blockmethod);
1044     }
1045     return 0;
1046     }
1047    
1048     int slv_tear_drop_reorder_block(slv_system_t sys, int32 bnum,
1049     int32 cutoff, int two,
1050     enum mtx_reorder_method blockmethod)
1051     {
1052     struct rel_relation **rp;
1053     struct var_variable **vp;
1054     const mtx_block_t *b;
1055     dof_t *d;
1056     mr_reorder_t *mrsys;
1057     mtx_matrix_t mtx;
1058     mtx_region_t reg;
1059     mtx_coord_t coord;
1060     int32 c,vlen,rlen,modcount;
1061     var_filter_t vf;
1062     rel_filter_t rf;
1063    
1064     if (sys==NULL) return 1;
1065     rlen = slv_get_num_solvers_rels(sys);
1066     vlen = slv_get_num_solvers_vars(sys);
1067     if (rlen ==0 || vlen == 0) return 1;
1068    
1069     rp = slv_get_solvers_rel_list(sys);
1070     vp = slv_get_solvers_var_list(sys);
1071     assert(rp!=NULL);
1072     assert(vp!=NULL);
1073     rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
1074     rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK| REL_ACTIVE);
1075     vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
1076     vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
1077    
1078     mtx = mtx_create();
1079     mtx_set_order(mtx,MAX(rlen,vlen));
1080     b = slv_get_solvers_blocks(sys);
1081     assert(b!=NULL); /* probably shouldn't be an assert here ... */
1082     if (bnum <0 || bnum >= b->nblocks || b->block == NULL) {
1083     mtx_destroy(mtx);
1084     return 1;
1085     }
1086     reg = b->block[bnum];
1087     /* do the flag setting even if we don't do the reorder */
1088     for (c=reg.col.low; c<=reg.col.high; c++) {
1089     var_set_in_block(vp[c],1);
1090     }
1091     for (c=reg.row.low; c<=reg.row.high; c++) {
1092     rel_set_in_block(rp[c],1);
1093     }
1094     if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
1095     mtx_destroy(mtx);
1096     return 1; /* must be square */
1097     }
1098     if (reg.row.high - reg.row.low < 3) {
1099     mtx_destroy(mtx);
1100     return 0; /* must be 3x3 or bigger to have any effect */
1101     }
1102    
1103     if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
1104     FPRINTF(stderr,
1105     "slv_tear_drop_reorder_block: failure in creating incidence matrix.\n");
1106     mtx_destroy(mtx);
1107     return 1;
1108     }
1109     /* verify that block has no empty columns, though not checking diagonal */
1110     for (c = reg.row.low; c <= reg.row.high; c++) {
1111     coord.col = mtx_FIRST;
1112     coord.row = c;
1113     if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
1114     mtx_destroy(mtx);
1115     FPRINTF(stderr, "slv_tear_drop_reorder_block: empty row (%d) found.\n",c);
1116     return 1;
1117     }
1118     coord.row = mtx_FIRST;
1119     coord.col = c;
1120     if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
1121     FPRINTF(stderr, "slv_tear_drop_reorder_block: empty col (%d) found.\n",c);
1122     mtx_destroy(mtx);
1123     return 1;
1124     }
1125     }
1126    
1127     modcount = slv_get_num_models(sys);
1128     mrsys = mr_reorder_create(sys,mtx,modcount);
1129     if (mrsys==NULL) return 1;
1130     mrsys->cutoff = cutoff;
1131     g_blockmethod = blockmethod;
1132     if (!two) {
1133     mr_bisect_partition(mrsys,&reg,1,tear_subreorder);
1134     } else {
1135     mr_bisect_partition2(mrsys,&reg,1,tear_subreorder);
1136     }
1137     #if 1
1138     FPRINTF(stderr,"Model-based reordering: (tear-style = %d)\n",two);
1139     FPRINTF(stderr,"Min. tearable size:\t%d\n",mrsys->cutoff);
1140     FPRINTF(stderr,"Tears:\t\t%d\n",mrsys->ntears);
1141     FPRINTF(stderr,"Partitionings:\t%d\n",mrsys->nblts);
1142     #endif
1143     reg = b->block[bnum]; /* bisect likely munged reg */
1144     if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
1145     mtx_destroy(mtx);
1146     mr_reorder_destroy(mrsys);
1147     return 2;
1148     }
1149     if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
1150     mtx_destroy(mtx);
1151     mr_reorder_destroy(mrsys);
1152     return 2;
1153     }
1154     d = slv_get_dofdata(sys);
1155     d->reorder.block_reordering = 2; /* tear_drop_baa */
1156    
1157     mtx_destroy(mtx);
1158     mr_reorder_destroy(mrsys);
1159     return 0;
1160     }
1161    
1162     int slv_insure_bounds(slv_system_t sys,int32 lo,int32 hi, FILE *mif)
1163     {
1164     real64 val,low,high;
1165     int32 c,nchange=0;
1166     struct var_variable *var, **vp;
1167    
1168     vp = slv_get_solvers_var_list(sys);
1169     if (vp==NULL) return -1;
1170     for (c= lo; c <= hi; c++) {
1171     var = vp[c];
1172     low = var_lower_bound(var);
1173     high = var_upper_bound(var);
1174     val = var_value(var);
1175     if( low > high ) {
1176     if (mif!=NULL) {
1177     FPRINTF(mif,"Bounds for variable ");
1178     var_write_name(sys,var,mif);
1179     FPRINTF(mif," are inconsistent [%g,%g].\n",low,high);
1180     FPRINTF(mif,"Bounds will be swapped.\n");
1181     }
1182     var_set_upper_bound(var, low);
1183     var_set_lower_bound(var, high);
1184     low = var_lower_bound(var);
1185     high = var_upper_bound(var);
1186     nchange++;
1187     }
1188    
1189     if( low > val ) {
1190     if (mif!=NULL) {
1191     FPRINTF(mif,"Variable ");
1192     var_write_name(sys,var,mif);
1193     FPRINTF(mif," was initialized below its lower bound.\n");
1194     FPRINTF(mif,"It will be moved to its lower bound.\n");
1195     }
1196     var_set_value(var, low);
1197     } else {
1198     if( val > high ) {
1199     if (mif!=NULL) {
1200     FPRINTF(mif,"Variable ");
1201     var_write_name(sys,var,mif);
1202     FPRINTF(mif," was initialized above its upper bound.\n");
1203     FPRINTF(mif,"It will be moved to its upper bound.\n");
1204     }
1205     var_set_value(var, high);
1206     nchange++;
1207     }
1208     }
1209     }
1210     return nchange;
1211     }
1212    
1213     void slv_check_bounds(const slv_system_t sys,int32 lo,int32 hi,
1214     FILE *mif,const char *label)
1215     {
1216     real64 val,low,high;
1217     int32 c,len;
1218     struct var_variable *var, **vp;
1219     static char defaultlabel[] = "";
1220    
1221     if (label==NULL) label = defaultlabel;
1222     if (sys==NULL || mif == NULL) return;
1223     vp = slv_get_solvers_var_list(sys);
1224     if (vp==NULL) return;
1225     len = slv_get_num_solvers_vars(sys);
1226     if (lo > len || hi > len) {
1227     FPRINTF(stderr,"slv_check_bounds miscalled\n");
1228     return;
1229     }
1230     for (c= lo; c <= hi; c++) {
1231     var = vp[c];
1232     low = var_lower_bound(var);
1233     high = var_upper_bound(var);
1234     val = var_value(var);
1235     if( low > high ) {
1236     FPRINTF(mif,"Bounds for %s variable ",label);
1237     var_write_name(sys,var,mif);
1238     FPRINTF(mif,"\nare inconsistent [%g,%g].\n",low,high);
1239     }
1240    
1241     if( low > val ) {
1242     FPRINTF(mif,"%s variable ",label);
1243     var_write_name(sys,var,mif);
1244     FPRINTF(mif,"\nwas initialized below its lower bound.\n");
1245     } else {
1246     if( val > high ) {
1247     FPRINTF(mif,"%s variable ",label);
1248     var_write_name(sys,var,mif);
1249     FPRINTF(mif," was initialized above its upper bound.\n");
1250     }
1251     }
1252     }
1253     return;
1254     }
1255    
1256     int SlvRegisterStandardClients(void)
1257     {
1258     int nclients = 0;
1259     int status;
1260     #ifdef STATIC_SLV
1261     status = slv_register_client(slv0_register,NULL,NULL);
1262     if (status) {
1263     FPRINTF(stderr,"Unable to register slv0.\n");
1264     } else {
1265     nclients++;
1266     }
1267     #endif
1268     #ifdef STATIC_MINOS
1269     status = slv_register_client(slv1_register,NULL,NULL);
1270     if (status) {
1271     FPRINTF(stderr,"Unable to register slv1.\n");
1272     } else {
1273     nclients++;
1274     }
1275     #endif
1276     #ifdef STATIC_QRSLV
1277     status = slv_register_client(slv3_register,NULL,NULL);
1278     if (status) {
1279     FPRINTF(stderr,"Unable to register slv3.\n");
1280     } else {
1281     nclients++;
1282     }
1283     #endif
1284     #ifdef STATIC_CSLV
1285     status = slv_register_client(slv4_register,NULL,NULL);
1286     if (status) {
1287     FPRINTF(stderr,"Unable to register slv4.\n");
1288     } else {
1289     nclients++;
1290     }
1291     #endif
1292     #ifdef STATIC_LSSLV
1293     status = slv_register_client(slv5_register,NULL,NULL);
1294     if (status) {
1295     FPRINTF(stderr,"Unable to register slv5.\n");
1296     } else {
1297     nclients++;
1298     }
1299     #endif
1300     #ifdef STATIC_MPS
1301     status = slv_register_client(slv6_register,NULL,NULL);
1302     if (status) {
1303     FPRINTF(stderr,"Unable to register slv6.\n");
1304     } else {
1305     nclients++;
1306     }
1307     #endif
1308     #ifdef STATIC_NGSLV
1309     status = slv_register_client(slv7_register,NULL,NULL);
1310     if (status) {
1311     FPRINTF(stderr,"Unable to register slv7.\n");
1312     } else {
1313     nclients++;
1314     }
1315     #endif
1316     #ifdef STATIC_OPTSQP
1317     status = slv_register_client(slv2_register,NULL,NULL);
1318     if (status) {
1319     FPRINTF(stderr,"Unable to register slv2.\n");
1320     } else {
1321     nclients++;
1322     }
1323     #endif
1324     #if (defined(STATIC_CONOPT) || defined(DYNAMIC_CONOPT))
1325     status = slv_register_client(slv8_register,NULL,NULL);
1326     if (status) {
1327     FPRINTF(stderr,"Unable to register slv8.\n");
1328     } else {
1329     nclients++;
1330     }
1331     #endif
1332     #ifdef STATIC_CMSLV
1333     status = slv_register_client(slv9_register,NULL,NULL);
1334     if (status) {
1335     FPRINTF(stderr,"Unable to register slv9.\n");
1336     } else {
1337     nclients++;
1338     }
1339     #endif
1340     #ifdef STATIC_LRSLV
1341     status = slv_register_client(slv9a_register,NULL,NULL);
1342     if (status) {
1343     FPRINTF(stderr,"Unable to register slv9a.\n");
1344     } else {
1345     nclients++;
1346     }
1347     #endif
1348     return nclients;
1349     }
1350    
1351    
1352    
1353     /************************************************************************\
1354     Output Assignment and partitiong in Logical Relations
1355     \************************************************************************/
1356    
1357     #define MLIMDEBUG 0
1358     /*
1359     * Populates a matrix according to the sys solvers_dvars, solvers_logrels
1360     * lists and the filters given. mtx given must be created (not null)
1361     * and with order >= MAX( slv_get_num_solvers_logrels(sys),
1362     * slv_get_num_solvers_dvars(sys));
1363     * returns 0 if went ok.
1364     */
1365     int slv_std_make_log_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
1366     dis_filter_t *dvf,logrel_filter_t *lrf)
1367     {
1368     #if MLIMDEBUG
1369     FILE *fp;
1370     #endif
1371     int32 lr, lrlen,ord;
1372     struct logrel_relation **lrp;
1373    
1374     if (sys==NULL || mtx == NULL || dvf == NULL || lrf == NULL) {
1375     FPRINTF(stderr,"make_log_incidence called with null\n");
1376     return 1;
1377     }
1378     lrp = slv_get_solvers_logrel_list(sys);
1379     assert(lrp!=NULL);
1380     lrlen = slv_get_num_solvers_logrels(sys);
1381     ord = MAX(lrlen,slv_get_num_solvers_dvars(sys));
1382     if (ord > mtx_order(mtx)) {
1383     FPRINTF(stderr,"make_incidence called with undersized matrix\n");
1384     return 2;
1385     }
1386     for (lr=0; lr < lrlen; lr++) {
1387     if (logrel_apply_filter(lrp[lr],lrf)) {
1388     logrelman_get_incidence(lrp[lr],dvf,mtx);
1389     }
1390     }
1391     #if MLIMDEBUG
1392     fp = fopen("/tmp/mim.plot","w+");
1393     if (fp !=NULL) {
1394     mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1395     fclose(fp);
1396     }
1397     #endif
1398     return 0;
1399     }
1400    
1401     #if USEDCODE
1402    
1403     /* returns 0 if successful, 1 if insufficient memory. Does not change
1404     * the data in mtx. Orders the solvers_dvar list of the system to
1405     * match the permutation on the given mtx. It is assumed that
1406     * the org cols of mtx == dvar list position. Only the
1407     * dvars in range lo to hi of the dvar list are affected and
1408     * only these dvars should appear in the same org column range of
1409     * the mtx. This should not be called on blocks less than 3x3.
1410     */
1411     static int reindex_dvars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
1412     const mtx_matrix_t mtx)
1413     {
1414     struct dis_discrete **dvtmp, **dvp;
1415     int32 c,v,vlen;
1416    
1417     if (lo >= hi +1) return 0; /* job too small */
1418     dvp = slv_get_solvers_dvar_list(sys);
1419     vlen = slv_get_num_solvers_dvars(sys);
1420     /* on dvtmp we DONT have the terminating null */
1421     dvtmp = (struct dis_discrete **)
1422     ascmalloc(vlen*sizeof(struct dis_discrete *));
1423     if (dvtmp == NULL) {
1424     return 1;
1425     }
1426     /* copy pointers to dvtmp in order desired and copy back rather than sort */
1427     for (c=lo;c<=hi;c++) {
1428     v = mtx_col_to_org(mtx,c);
1429     dvtmp[c] = dvp[v];
1430     }
1431     /* copying back and re-sindexing */
1432     for (c=lo;c<=hi;c++) {
1433     dvp[c] = dvtmp[c];
1434     dis_set_sindex(dvp[c],c);
1435     }
1436     ascfree(dvtmp);
1437     return 0;
1438     }
1439    
1440     /* returns 0 if successful, 1 if insufficient memory. Does not change
1441     * the data in mtx. Orders the solvers_logrel list of the system to
1442     * match the permutation on the given mtx. It is assumed that
1443     * the org rows of mtx == logrel list position. Only the
1444     *logrels in range lo to hi of the logrel list are affected and
1445     * only these logrels should appear in the same org row range of
1446     * the input mtx. This should not be called on blocks less than 3x3.
1447     */
1448     static int reindex_logrels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
1449     const mtx_matrix_t mtx)
1450     {
1451     struct logrel_relation **lrtmp, **lrp;
1452     int32 c,v,rlen;
1453    
1454     lrp = slv_get_solvers_logrel_list(sys);
1455     rlen = slv_get_num_solvers_logrels(sys);
1456     /* on lrtmp we DONT have the terminating null */
1457     lrtmp = (struct logrel_relation **)
1458     ascmalloc(rlen*sizeof(struct logrel_relation *));
1459     if (lrtmp == NULL) {
1460     return 1;
1461     }
1462     /* copy pointers to lrtmp in order desired and copy back rather than sort.
1463     * do this only in the row range of interest. */
1464     for (c=lo;c<=hi;c++) {
1465     v = mtx_row_to_org(mtx,c);
1466     #if RIDEBUG
1467     FPRINTF(stderr,"Old logrel sindex(org) %d becoming sindex(cur) %d\n",v,c);
1468     #endif
1469     lrtmp[c] = lrp[v];
1470     }
1471     /* copying back and re-sindexing */
1472     for (c=lo;c<=hi;c++) {
1473     lrp[c] = lrtmp[c];
1474     logrel_set_sindex(lrp[c],c);
1475     }
1476     ascfree(lrtmp);
1477     return 0;
1478     }
1479    
1480     #endif /* USEDCODE */
1481    
1482     /*
1483     * returns 0 if ok, OTHERWISE if madness detected.
1484     */
1485     #define SLBPDEBUG 0
1486     int slv_log_block_partition(slv_system_t sys)
1487     {
1488     #if SLBPDEBUG
1489     FILE *fp;
1490     #endif
1491     struct logrel_relation **lrp;
1492     struct dis_discrete **dvp;
1493     mtx_region_t *newblocks;
1494     dof_t *d;
1495     int32 nrow,ncol;
1496     int32 ncolpfix, nrowpun, vlenmnv;
1497     mtx_matrix_t mtx;
1498     int32 order, rank;
1499     int32 c,len,dvlen,r,lrlen;
1500     dis_filter_t dvf;
1501     logrel_filter_t lrf;
1502    
1503     lrp = slv_get_solvers_logrel_list(sys);
1504     dvp = slv_get_solvers_dvar_list(sys);
1505     lrlen = slv_get_num_solvers_logrels(sys);
1506     dvlen = slv_get_num_solvers_dvars(sys);
1507     if (lrlen ==0 || dvlen == 0) return 1;
1508     order = MAX(lrlen,dvlen);
1509    
1510     lrf.matchbits = (LOGREL_ACTIVE);
1511     lrf.matchvalue = (LOGREL_ACTIVE);
1512     dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1513     dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1514    
1515     /* nrow plus unincluded logrels */
1516     nrowpun = slv_count_solvers_logrels(sys,&lrf);
1517     /* ncol plus fixed discrete vars */
1518     ncolpfix = slv_count_solvers_dvars(sys,&dvf);
1519    
1520     dvf.matchbits = (DIS_BVAR);
1521     dvf.matchvalue = 0;
1522    
1523     /* dvlen minus non boolean vars */
1524     vlenmnv = dvlen - slv_count_solvers_dvars(sys,&dvf);
1525    
1526     lrf.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
1527     lrf.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
1528     dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_FIXED | DIS_ACTIVE);
1529     dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);
1530    
1531     nrow = slv_count_solvers_logrels(sys,&lrf);
1532     ncol = slv_count_solvers_dvars(sys,&dvf);
1533    
1534     mtx = slv_get_sys_mtx(sys);
1535     mtx_set_order(mtx,order);
1536    
1537     if (slv_std_make_log_incidence_mtx(sys,mtx,&dvf,&lrf)) {
1538     FPRINTF(stderr,
1539     "slv_log_block_partition: failure in creating incidence matrix.\n");
1540     mtx_destroy(mtx);
1541     return 1;
1542     }
1543    
1544     mtx_output_assign(mtx,lrlen,dvlen);
1545     rank = mtx_symbolic_rank(mtx);
1546     if (rank == 0 ) return 1;
1547     if (rank < nrow) {
1548     FPRINTF(stderr,"rank<nrow in slv_log_block_partition\n");
1549     FPRINTF(stderr,"dependent logical relations ? \n");
1550     }
1551     if (rank < ncol) {
1552     if ( nrow != rank) {
1553     FPRINTF(stderr,"rank<ncol and nrow!=rank in slv_log_block_partition.\n");
1554     FPRINTF(stderr,"Excess of columns ? \n");
1555     } else {
1556     FPRINTF(stderr,"rank<ncol but nrow==rank in slv_log_block_partition.\n");
1557     FPRINTF(stderr,"Degrees of freedom ? \n");
1558     }
1559     }
1560     if (ncol == nrow) {
1561     if (ncol == rank) {
1562     FPRINTF(stderr,"\n");
1563     FPRINTF(stderr,
1564     "System of logical relations does not need Logic Inference \n");
1565     }
1566     if (ncol != rank) {
1567     FPRINTF(stderr,"but ncol!=rank.\n");
1568     FPRINTF(stderr,"Rank deficient ? \n");
1569     } else {
1570     FPRINTF(stderr,"\n");
1571     }
1572     }
1573     mtx_partition(mtx);
1574     /* copy the block list. there is at least one if rank >=1 as above */
1575     len = mtx_number_of_blocks(mtx);
1576     newblocks = (mtx_region_t *)ascmalloc(len*sizeof(mtx_region_t));
1577     if (newblocks == NULL) {
1578     mtx_destroy(mtx);
1579     return 2;
1580     }
1581     for (c = 0 ; c < len; c++) {
1582     mtx_block(mtx,c,&(newblocks[c]));
1583     }
1584     slv_set_solvers_log_blocks(sys,len,newblocks); /* give it to the system */
1585     d = slv_get_log_dofdata(sys);
1586     d->structural_rank = rank;
1587     d->n_rows = nrow;
1588     d->n_cols = ncol;
1589     d->n_fixed = ncolpfix - ncol;
1590     d->n_unincluded = nrowpun - nrow;
1591     d->reorder.partition = 1; /* yes */
1592     d->reorder.basis_selection = 0; /* none yet */
1593     d->reorder.block_reordering = 0; /* none */
1594    
1595     #if SLBPDEBUG
1596     fp = fopen("/tmp/sbp1.plot","w+");
1597     if (fp !=NULL) {
1598     mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1599     fclose(fp);
1600     }
1601     #endif
1602    
1603    
1604     len = vlenmnv - 1; /* row of last possible boolean variable */
1605     for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
1606     r = mtx_col_to_org(mtx,c);
1607     if ( dis_fixed(dvp[r]) && dis_active(dvp[r])) {
1608     c++;
1609     } else {
1610     mtx_swap_cols(mtx,c,len);
1611     len--;
1612     }
1613     }
1614    
1615     /* Now, leftover columns are where we want them.
1616     * That is, unassigned incidence is next to assigned and fixed on right with
1617     * no columns which do not correspond to dis vars in the range 0..dvlen-1.
1618     */
1619     /* there aren't any nonincident dvars in the solvers dvarlist, so
1620     * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
1621    
1622     /* rows 0->lrlen-1 should be actual logrels, though, and not phantoms.
1623     * At this point we should have
1624     * rows - assigned
1625     * - unassigned w/incidence
1626     * - unincluded and included w/out free incidence and inactive mixed.
1627     */
1628    
1629     /* sort unincluded logrelations to bottom of good region */
1630     /* this is needed because the idiot user may have fixed everything */
1631     /* in a logrelation, but we want those rows not to be confused with */
1632     /* unincluded rows.*/
1633     len = lrlen - 1; /* row of last possible relation */
1634     /* sort the inactive out of the unassigned/unincluded */
1635     for (c=rank; len > c ; ) {
1636     r = mtx_row_to_org(mtx,c);
1637     if ( logrel_active(lrp[r])) {
1638     c++;
1639     } else {
1640     mtx_swap_rows(mtx,c,len);
1641     len--;
1642     }
1643     }
1644    
1645     /* sort the unincluded out of the unassigned */
1646     len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
1647     for (c=rank; len > c ; ) {
1648     r = mtx_row_to_org(mtx,c);
1649     if ( logrel_included(lrp[r]) ) {
1650     c++;
1651     } else {
1652     mtx_swap_rows(mtx,c,len);
1653     len--;
1654     }
1655     }
1656    
1657     #if SBPDEBUG
1658     fp = fopen("/tmp/sbp2.plot","w+");
1659     if (fp !=NULL) {
1660     mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
1661     fclose(fp);
1662     }
1663     #endif
1664    
1665     /* These functions are called for slv_block_partitioning, but I do not
1666     * know if I need them here:
1667     *
1668     * now, at last we have cols in the order we want the lists to
1669     * be handed to the solver. So, we need to reset the dvar sindex values
1670     * and reorder the lists.
1671     */
1672    
1673     #if USEDCODE
1674     if (reindex_dvars_from_mtx(sys,0,dvlen-1,mtx)) {
1675     mtx_destroy(mtx);
1676     return 2;
1677     }
1678     #endif /* USEDCODE */
1679    
1680     /*
1681     * now, at last we have rows jacobian in the order we want the lists to
1682     * be handed to the solvers. So, we need to reset the rel sindex values.
1683     */
1684    
1685     #if USEDCODE
1686     if (reindex_logrels_from_mtx(sys,0,lrlen-1,mtx)) {
1687     mtx_destroy(mtx);
1688     return 2;
1689     }
1690     #endif /* USEDCODE */
1691    
1692     return 0;
1693     }
1694    

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22