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

Contents of /trunk/ascend4/solver/model_reorder.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 17160 byte(s)
Setting up web subdirectory in repository
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