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

Annotation of /trunk/ascend4/solver/model_reorder.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, 9 months ago) by aw0a
File MIME type: text/x-csrc
File size: 17160 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 /*
2     * Model-based Reordering Routines
3     * by Benjamin Andrew Allan
4     * 6/22/96
5     * Version: $Revision: 1.11 $
6     * Version control file: $RCSfile: model_reorder.c,v $
7     * Date last modified: $Date: 1997/07/18 12:14:43 $
8     * Last modified by: $Author: mthomas $
9     * Author $$
10     * Copyright(C) 1996 Benjamin Andrew Allan
11     *
12     * This file is part of the ASCEND IV math programming system.
13     * These functions are part of a new design for feeding
14     * solvers from the ASCEND compiler.
15     *
16     * File to play MODEL-relation based reordering games.
17     * We're starting with a basic RBBD implementation.
18     * Assumptions:
19     * - If the MODELs are hierarchical, then they have been indexed in a
20     * tree bottom-up fashion.
21     * - Input is a square block that needs tearing and not a rectangle.
22     *
23     * The ASCEND IV math programming system is free software; you can redistribute
24     * it and/or modify it under the terms of the GNU General Public License as
25     * published by the Free Software Foundation; either version 2 of the
26     * License, or (at your option) any later version.
27     *
28     * The ASCEND IV math programming system is distributed in hope that it will be
29     * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
30     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
31     * General Public License for more details.
32     *
33     * You should have received a copy of the GNU General Public License along with
34     * the program; if not, write to the Free Software Foundation, Inc., 675
35     * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
36     * COPYING is found in ../compiler.
37     */
38    
39     #include "utilities/ascConfig.h"
40     #include "compiler/compiler.h"
41     #include "utilities/ascMalloc.h"
42     #include "general/list.h"
43     #include "solver/slv_types.h"
44     #include "solver/var.h"
45     #include "solver/rel.h"
46     #include "solver/discrete.h"
47     #include "solver/conditional.h"
48     #include "solver/logrel.h"
49     #include "solver/bnd.h"
50     #include "solver/mtx.h"
51     #include "solver/slv_common.h"
52     #include "solver/linsol.h"
53     #include "solver/linsolqr.h"
54     #include "solver/slv_client.h"
55     #include "solver/model_reorder.h"
56    
57    
58     mr_reorder_t *mr_reorder_create(slv_system_t slv, mtx_matrix_t mtx, int32 size)
59     {
60     mr_reorder_t *res;
61     if(slv==NULL || mtx == NULL || size < 1) {
62     return NULL;
63     }
64     res = (mr_reorder_t *)asccalloc(1,sizeof(mr_reorder_t));
65     if (res == NULL) return res;
66     res->slv = slv;
67     res->mtx = mtx;
68     res->cutoff = CUTOFFDEFAULT;
69     res->activelen = size;
70     res->local = (int32 *)asccalloc(2*size,sizeof(int32));
71     res->tmpmodlist = res->local+size; /* 2nd half of allocation is tmpmodlist */
72     res->active = (char *)asccalloc(size,sizeof(char));
73     if (res->local == NULL || res->active == NULL) {
74     if(res->local !=NULL) ascfree(res->local);
75     if(res->active !=NULL) ascfree(res->active);
76     ascfree(res);
77     res = NULL;
78     }
79     return res;
80     }
81    
82     void mr_reorder_destroy(mr_reorder_t *sys)
83     {
84     if (sys==NULL) return;
85     sys->slv=NULL;
86     sys->mtx=NULL;
87     ascfree(sys->local); /*tmpmodlist freed too, hereby */
88     ascfree(sys->active);
89     ascfree(sys);
90     }
91    
92     /* comparing distinct integers with qsort this yields an increasing
93     integer list (or else integer overflow) */
94     static int intcompare(int *i, int *j)
95     {
96     return(*i - *j);
97     }
98    
99     int mr_bisect_partition(mr_reorder_t *sys, mtx_region_t *reg,
100     int top, MRBlockReorderF rfunc)
101     {
102     #if MRDEBUG
103     FILE *fptr;
104     static mtx_region_t graph_reg;
105     static int fnum;
106     static char fname[80];
107     #endif
108     mtx_block_t mb;
109     mtx_block_t *mbptr;
110     mtx_region_t *block;
111     struct rel_relation **rp;
112     mtx_coord_t coord;
113     int partition, newpart,size,threshold, status = 0;
114     int32 k,row,mod,rmax,stop,nrows;
115    
116     /* check all our pointers */
117    
118     if (sys==NULL || sys->slv == NULL || sys->mtx == NULL ||
119     sys->local == NULL || sys->active == NULL ||
120     sys->tmpmodlist == NULL || sys->activelen == 0 ||
121     sys->cutoff < 10 || rfunc == NULL
122     ) {
123     #if MRDEBUG
124     FPRINTF(stderr,"mr_bisect_partition miscalled.\n");
125     #endif
126     return 1;
127     }
128    
129     rp = slv_get_solvers_rel_list(sys->slv);
130    
131     /* CREATE block list as needed */
132     if (top==1) {
133     mbptr = &mb;
134     mb.nblocks = 1;
135     mb.block = reg;
136     #if MRDEBUG
137     graph_reg = *reg;
138     fnum=0;
139     #endif
140     } else {
141     mbptr = mtx_block_partition(sys->mtx,reg);
142     if (mbptr == NULL || (mbptr->block == NULL && mbptr->nblocks) ) {
143     if (mbptr != NULL) ascfree(mbptr);
144     return 1;
145     }
146     sys->nblts++;
147     }
148    
149     #if MRDEBUG
150     FPRINTF(stderr,"mtx_bisect_partition2 (%d,%d - %d,%d).\n",
151     reg->row.low,reg->col.low,reg->row.high,reg->col.high);
152     sprintf(fname,"/tmp/bisect1.plot.%d",fnum);
153     fptr = fopen(fname,"w+");
154     mtx_write_region_plot(fptr,sys->mtx,&graph_reg);
155     fclose(fptr);
156     fnum++;
157     #endif
158    
159     /* tear (or skip) the k remaining partitions */
160     for (k=0; k < mbptr->nblocks; k++) {
161     block = &(mbptr->block[k]);
162     size = block->row.high - block->row.low +1;
163    
164     /* skip 1x1, 2x2 blocks which must be full */
165     if (block->row.high - block->row.low < 3) {
166     continue; /* next k */
167     }
168    
169     /* blocks below cutoff the user reorders instead */
170     if (size < sys->cutoff) {
171     #if MRDEBUG
172     FPRINTF(stderr,"Calling user reorder on rows %d - %d.\n",
173     block->row.low,block->row.high);
174     #endif
175     rfunc(sys->slv,sys->mtx,block);
176     continue; /* next k */
177     }
178    
179     /* ok, here's where we do the work */
180     threshold = size/2; /* we want to divide in approximately half */
181     rmax = block->row.high;
182    
183     /* tag in block models and count in block relations */
184     for (row = block->row.low; row <= rmax; row++) {
185     mod = rel_model(rp[mtx_row_to_org(sys->mtx,row)]);
186     if (!(sys->active[mod])) {
187     sys->active[mod] = 1;
188     sys->tmpmodlist[sys->tmpmodused++] = mod;
189     }
190     sys->local[mod]++;
191     }
192    
193     /* if block is all one MODEL, tell user to eat it */
194     if (sys->tmpmodused < 2) {
195     #if MRDEBUG
196     FPRINTF(stderr,"Found oversized MODEL including relation\n");
197     rel_write_name(sys->slv,rp[mtx_row_to_org(sys->mtx,rmax)],stderr);
198     FPRINTF(stderr,"\nCalling user reorder on rows %d - %d.\n",
199     block->row.low,block->row.high);
200     #endif
201     /* rezero active and local arrays, just 1! */
202     sys->local[0] = 0;
203     sys->active[0] = 0;
204     rfunc(sys->slv,sys->mtx,block);
205     continue; /* next k */
206     }
207    
208     /* sort the MODEL indices in tmpmodlist. using this fact later
209     * enables us to punt the rel_partition flag in favor of just
210     * checking MODEL index of a row against the last MODEL in the
211     * first partition when hunting in incidence for complicating cols.
212     * We might want to reverse things if the tree was numbered top
213     * down rather than bottom up.
214     * It would be nice if we could prove this sort is not needed,
215     * but I haven't had the time yet to try.
216     * A slightly different, probably bigger, data structure
217     * might also allow the sort to be avoided.
218     * Probably we should create a doubly linked list version
219     * of tmpmodlist and insert sorted as it is created,
220     * then map the final thing down to the vector needed.
221     * Who knows? It's a small list in any case.
222     */
223     qsort(sys->tmpmodlist,sys->tmpmodused,sizeof(int),
224     (int (*)(const void *, const void *))intcompare);
225    
226     /* now try to split tmpmodlist approximately in relation count half */
227     stop = 0;
228     nrows = 0;
229     /* local[mod] should be the size of the last MODEL in the first half
230     * on loop exit.
231     */
232     while(stop < sys->tmpmodused && nrows < threshold) {
233     /* stop at the first MODEL division after the halfway point */
234     mod = sys->tmpmodlist[stop++];
235     nrows += sys->local[mod];
236     }
237     /* check if the division BEFORE the halfway point is closer to center
238     * but not at beginning.
239     * This should help with large, centered blocks that OTHERWISE
240     * would put us way past half and result in overtearing.
241     * |------d1---------|t|---------------------d2------------|
242     * d1: threshold-(nrows-local[mod]) d2: nrows-threshold
243     * the abs should not be needed, but to be safe.
244     */
245     if ( stop > 1 &&
246     abs(nrows - threshold) > abs(threshold - (nrows - sys->local[mod])) ) {
247     nrows -= sys->local[mod];
248     stop--;
249     mod = sys->tmpmodlist[stop-1];
250     }
251    
252     /* now all relations with MODEL number > mod are in the
253     second partition with MODEL number <= mod are in the first. */
254    
255     /* now we're going to identify and permute out the tears,
256     diddling with the block. rmax preserves block->row.high */
257     coord.col = block->row.low;
258     while (coord.col <= block->row.high) {
259     coord.row = mtx_FIRST;
260     partition = -1; /* neither partition seen */
261     while ( mtx_next_in_col(sys->mtx,&coord,&(block->row)),
262     coord.row != mtx_LAST) {
263    
264     /* check if this column crossed partitions */
265     newpart = (rel_model(rp[mtx_row_to_org(sys->mtx,coord.row)]) > mod);
266     if (partition >= 0 && partition != newpart) {
267     /* this is a tear. */
268     sys->ntears++;
269     /* symmetrically permute tear to block outer edge */
270     mtx_swap_cols(sys->mtx,coord.col,block->row.high);
271     mtx_swap_rows(sys->mtx,coord.col,block->row.high);
272     rel_set_torn(rp[mtx_row_to_org(sys->mtx,block->row.high)],1);
273     /* reduce the block row range, break next_in loop.
274     * with any luck, this stops us from double tearing, though
275     * in a completely arbitrary way. Another covert tie-breaking.
276     */
277     block->row.high--;
278     break; /* next coord.col. */
279     }
280     /* it doesn't cross a partition yet, keep looking */
281     partition = newpart; /* slightly redundant. cheaper than checking */
282     }
283     coord.col++;
284     }
285     block->col.high = block->row.high;
286    
287     /* rezero active and local arrays so as not to confuse recursion */
288     while (sys->tmpmodused > 0) {
289     mod = sys->tmpmodlist[--(sys->tmpmodused)];
290     sys->local[mod] = 0;
291     sys->active[mod] = 0;
292     }
293     /* there are no flags on rels to rezero */
294    
295     /* attack block left after tearing. */
296     status += mr_bisect_partition(sys,block,0,rfunc);
297     }
298    
299     /* if block list is locally created, destroy it */
300     if (!top) {
301     mtx_destroy_blocklist(mbptr);
302     }
303     return status;
304     }
305    
306    
307     int mr_bisect_partition2(mr_reorder_t *sys, mtx_region_t *reg,
308     int top, MRBlockReorderF rfunc)
309     {
310     #if MRDEBUG
311     FILE *fptr;
312     static mtx_region_t graph_reg;
313     static int fnum;
314     static char fname[80];
315     #endif
316     mtx_block_t mb;
317     mtx_block_t *mbptr;
318     mtx_region_t *block;
319     struct rel_relation **rp;
320     mtx_coord_t coord;
321     int partition, newpart,size,threshold, status = 0;
322     int32 k,row,mod,rmax,stop,nrows,nexttear;
323    
324     /* check all our pointers */
325    
326     if (sys==NULL || sys->slv == NULL || sys->mtx == NULL ||
327     sys->local == NULL || sys->active == NULL ||
328     sys->tmpmodlist == NULL || sys->activelen == 0 ||
329     sys->cutoff < 10 || rfunc == NULL
330     ) {
331     #if MRDEBUG
332     FPRINTF(stderr,"mr_bisect_partition2 miscalled.\n");
333     #endif
334     return 1;
335     }
336    
337     rp = slv_get_solvers_rel_list(sys->slv);
338    
339     /* CREATE block list as needed */
340     if (top==1) {
341     mbptr = &mb;
342     mb.nblocks = 1;
343     mb.block = reg;
344     #if MRDEBUG
345     graph_reg = *reg;
346     fnum=0;
347     #endif
348     } else {
349     #if MRDEBUG
350     FPRINTF(stderr,"mtx_bisect_partition2 (%d,%d - %d,%d).\n",
351     reg->row.low,reg->col.low,reg->row.high,reg->col.high);
352     #endif
353     mbptr = mtx_block_partition(sys->mtx,reg);
354     if (mbptr == NULL || (mbptr->block == NULL && mbptr->nblocks) ) {
355     if (mbptr != NULL) ascfree(mbptr);
356     return 1;
357     }
358    
359     #if MRDEBUG
360     FPRINTF(stderr,"BLT List:\n");
361     for (k=0;k < mbptr->nblocks; k++) {
362     FPRINTF(stderr,"B: %d,%d - %d,%d\n",
363     mbptr->block[k].row.low,mbptr->block[k].col.low,
364     mbptr->block[k].row.high,mbptr->block[k].col.high);
365     }
366     #endif
367     sys->nblts++;
368     }
369    
370     #if MRDEBUG
371     FPRINTF(stderr,"mtx_bisect_partition2 (%d,%d - %d,%d).\n",
372     reg->row.low,reg->col.low,reg->row.high,reg->col.high);
373     sprintf(fname,"/tmp/bisect2.plot.%d",fnum);
374     fptr = fopen(fname,"w+");
375     mtx_write_region_plot(fptr,sys->mtx,&graph_reg);
376     fclose(fptr);
377     fnum++;
378     #endif
379    
380     /* tear (or skip) the k remaining partitions */
381     for (k=0; k < mbptr->nblocks; k++) {
382     block = &(mbptr->block[k]);
383     size = block->row.high - block->row.low +1;
384    
385     /* skip 1x1, 2x2 blocks which must be full */
386     if (block->row.high - block->row.low < 3) {
387     continue; /* next k */
388     }
389    
390     /* blocks below cutoff the user reorders instead */
391     if (size < sys->cutoff) {
392     #if MRDEBUG
393     FPRINTF(stderr,"Calling user reorder on rows %d - %d.\n",
394     block->row.low,block->row.high);
395     #endif
396     rfunc(sys->slv,sys->mtx,block);
397     continue; /* next k */
398     }
399    
400     /* ok, here's where we do the work */
401     threshold = size/2; /* we want to divide in approximately half */
402     rmax = block->row.high;
403    
404     #if MRDEBUG
405     FPRINTF(stderr,"Splitting block with threshold %d.\n", threshold);
406     #endif
407    
408     /* tag in block models and count in block relations */
409     for (row = block->row.low; row <= rmax; row++) {
410     mod = rel_model(rp[mtx_row_to_org(sys->mtx,row)]);
411     if (!(sys->active[mod])) {
412     sys->active[mod] = 1;
413     sys->tmpmodlist[sys->tmpmodused++] = mod;
414     }
415     sys->local[mod]++;
416     }
417    
418     /* if block is all one MODEL, tell user to eat it */
419     if (sys->tmpmodused < 2) {
420     #if MRDEBUG
421     FPRINTF(stderr,"Found oversized MODEL including relation\n");
422     rel_write_name(sys->slv,rp[mtx_row_to_org(sys->mtx,rmax)],stderr);
423     FPRINTF(stderr,"\nCalling user reorder on rows %d - %d.\n",
424     block->row.low,block->row.high);
425     #endif
426     /* rezero active and local arrays, just 1! */
427     sys->local[0] = 0;
428     sys->active[0] = 0;
429     rfunc(sys->slv,sys->mtx,block);
430     continue; /* next k */
431     }
432    
433     #if MRDEBUG
434     FPRINTF(stderr,"Number of models in block %d\n",sys->tmpmodused);
435     FPRINTF(stderr,"tearing block %d-%d\n",block->row.low,block->row.high);
436     #endif
437     /* sort the MODEL indices in tmpmodlist. using this fact later
438     * enables us to punt the rel_partition flag in favor of just
439     * checking MODEL index of a row against the last MODEL in the
440     * first partition.
441     * We might want to reverse things if the tree was numbered top
442     * down rather than bottom up.
443     */
444     qsort(sys->tmpmodlist,sys->tmpmodused,sizeof(int),
445     (int (*)(const void *, const void *))intcompare);
446    
447     /* now try to split tmpmodlist approximately in relation count half */
448     stop = 0;
449     nrows = 0;
450     while(stop < sys->tmpmodused && nrows < threshold) {
451     mod = sys->tmpmodlist[stop++];
452     nrows += sys->local[mod];
453     }
454     /* Now all relations with MODEL number > mod are in the
455     * second partition with MODEL number <= mod are in the first.
456     */
457     #if MRDEBUG
458     FPRINTF(stderr,"Number of rows in first half %d\n",nrows);
459     FPRINTF(stderr,"Number of models in first half %d\n",stop);
460     #endif
461    
462     /* now we're going to identify and permute out the tears,
463     * diddling with the block. rmax preserves block->row.high
464     */
465     coord.col = block->row.low;
466     nexttear = block->row.high;
467     while (coord.col <= nexttear) {
468     coord.row = mtx_FIRST;
469     partition = -1; /* neither partition seen */
470     while ( mtx_next_in_col(sys->mtx,&coord,&(block->row)),
471     coord.row != mtx_LAST) {
472    
473     /* check if this column crossed partitions */
474     newpart = (rel_model(rp[mtx_row_to_org(sys->mtx,coord.row)]) > mod);
475     #if MRDEBUG
476     if (newpart <0) {
477     FPRINTF(stderr,"unexpect newpart < 0 (%d)\n",coord.row);
478     }
479     #endif
480     if (partition >= 0 && partition != newpart) {
481     /* this is a tear. */
482     sys->ntears++;
483     /* symmetrically permute tear to block outer edge */
484     mtx_swap_cols(sys->mtx,coord.col,nexttear);
485     mtx_swap_rows(sys->mtx,coord.col,nexttear);
486     rel_set_torn(rp[mtx_row_to_org(sys->mtx,nexttear)],1);
487     /* reduce the block row range, break next_in loop.
488     * with any luck, this stops us from double tearing, though
489     * in a completely arbitrary way. Another covert tie-breaking.
490     */
491     #if MRDEBUG
492     FPRINTF(stderr,"tear row %d\n",nexttear);
493     #endif
494     nexttear--;
495     break; /* next coord.col. */
496     }
497     /* it doesn't cross a partition yet, keep looking */
498     partition = newpart; /* slightly redundant. cheaper than checking */
499     }
500     coord.col++;
501     }
502     #if MRDEBUG
503     FPRINTF(stderr,"TORE %d rows\n", block->col.high-nexttear);
504     #endif
505     block->col.high = block->row.high = nexttear;
506    
507     /* rezero active and local arrays so as not to confuse recursion */
508     while (sys->tmpmodused > 0) {
509     mod = sys->tmpmodlist[--(sys->tmpmodused)];
510     sys->local[mod] = 0;
511     sys->active[mod] = 0;
512     }
513     /* there are no flags on rels to rezero */
514    
515     /* attack block left after tearing. */
516     status += mr_bisect_partition2(sys,block,0,rfunc);
517     }
518    
519     /* if block list is locally created, destroy it */
520     if (!top) {
521     mtx_destroy_blocklist(mbptr);
522     }
523     return status;
524     }

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