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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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