/[ascend]/branches/ksenija2/ascend/system/block.c
ViewVC logotype

Contents of /branches/ksenija2/ascend/system/block.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2867 - (show annotations) (download) (as text)
Mon Mar 23 02:31:01 2015 UTC (7 years, 3 months ago) by jpye
File MIME type: text/x-csrc
File size: 34428 byte(s)
minor changes to debug output.
comments relating to disused var flags.

1 /* ASCEND modelling environment
2 Copyright (C) 1998, 2006-2007 Carnegie Mellon University
3 Copyright (C) 1996 Benjamin Andrew Allan
4
5 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, see <http://www.gnu.org/licenses/>.
17 *//* @file
18 Block partitioning implementation (for real-valued variables)
19 *//*
20 slv_stdcalls.c by Benjamin Andrew Allan, May 1996.
21 chopped out block-paritioning into separate file, John Pye, Jan 2007.
22 Last in CVS: $Revision: 1.28 $ $Date: 1998/06/16 16:53:04 $ $Author: mthomas $
23 */
24
25 #include "block.h"
26 #include <ascend/general/ascMalloc.h>
27 #include <ascend/utilities/ascPrint.h>
28 #include <ascend/general/panic.h>
29 #include <ascend/general/mathmacros.h>
30 #include "slv_client.h"
31 #include "slv_stdcalls.h"
32 #include "model_reorder.h"
33
34 /* #define REINDEX_DEBUG */
35 /* #define BLOCKPARTITION_DEBUG */
36 /* #define CUT_DEBUG */
37
38 /* global to get around the mr header (for tear_subreorder) */
39 static
40 enum mtx_reorder_method g_blockmethod = mtx_UNKNOWN;
41
42 /*-----------------------------------------------------------------------------
43 VAR/REL ORDERING (TO MATCH MATRIX PERMUTATIONS)
44 */
45
46 /**
47 Orders the solvers_var list of the system to match the permutation
48 on the given mtx. Does not change the data in mtx.
49
50 @return 0 on success, 1 on out-of-memory
51
52 It is assumed that the org cols of mtx == var list position. Only the
53 vars in range lo to hi of the var list are affected and
54 only these vars should appear in the same org column range of
55 the mtx. This should not be called on blocks less than 3x3.
56 */
57 static int reindex_vars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
58 const mtx_matrix_t mtx)
59 {
60 struct var_variable **vtmp, **vp;
61 int32 c,v,vlen;
62
63 if (lo >= hi +1) return 0; /* job too small */
64 vp = slv_get_solvers_var_list(sys);
65 vlen = slv_get_num_solvers_vars(sys);
66 /* on vtmp we DONT have the terminating null */
67 vtmp = ASC_NEW_ARRAY(struct var_variable *,vlen);
68 if (vtmp == NULL) {
69 return 1;
70 }
71 /* copy pointers to vtmp in order desired and copy back rather than sort */
72 for (c=lo;c<=hi;c++) {
73 v = mtx_col_to_org(mtx,c);
74 vtmp[c] = vp[v];
75 }
76 /* copying back and re-sindexing */
77 for (c=lo;c<=hi;c++) {
78 vp[c] = vtmp[c];
79 var_set_sindex(vp[c],c);
80 }
81 ascfree(vtmp);
82 return 0;
83 }
84 /**
85 Orders the solvers_rel list of the system to match the permutation
86 on the given mtx. Does not change the data in mtx.
87
88 @return 0 on success, 1 on out-of-memory
89
90 It is assumed that the org rows of mtx == rel list position. Only the
91 rels in range lo to hi of the rel list are affected and
92 only these rels should appear in the same org row range of
93 the input mtx. This should not be called on blocks less than 3x3.
94 */
95 static int reindex_rels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
96 const mtx_matrix_t mtx)
97 {
98 struct rel_relation **rtmp, **rp;
99 int32 c,v,rlen;
100
101 rp = slv_get_solvers_rel_list(sys);
102 rlen = slv_get_num_solvers_rels(sys);
103 /* on rtmp we DONT have the terminating null */
104 rtmp = ASC_NEW_ARRAY(struct rel_relation*,rlen);
105 if (rtmp == NULL) {
106 return 1;
107 }
108 /* copy pointers to rtmp in order desired and copy back rather than sort.
109 * do this only in the row range of interest. */
110 for (c=lo;c<=hi;c++) {
111 v = mtx_row_to_org(mtx,c);
112 #ifdef REINDEX_DEBUG
113 if(c!=v){
114 CONSOLE_DEBUG("Old rel sindex (org) %d becoming sindex (cur) %d\n",v,c);
115 }
116 #endif
117 rtmp[c] = rp[v];
118 }
119 /* copying back and re-sindexing */
120 for (c=lo;c<=hi;c++) {
121 rp[c] = rtmp[c];
122 rel_set_sindex(rp[c],c);
123 }
124 ascfree(rtmp);
125 return 0;
126 }
127
128 /*------------------------------------------------------------------------------
129 PARTITIONING INTO BLOCK LOWER/UPPER TRIANGULAR FORM
130 */
131
132 /**
133 Perform var and rel reordering to achieve block form.
134
135 @param uppertrianguler if non-zero, BUT form. if 0, make BLT form.
136
137 @callergraph
138 */
139 int slv_block_partition_real(slv_system_t sys,int uppertriangular){
140 #ifdef BLOCKPARTITION_DEBUG
141 FILE *fp;
142 #endif
143 struct rel_relation **rp;
144 struct var_variable **vp;
145 mtx_region_t *newblocks;
146 dof_t *d;
147 int32 nrow,ncol;
148 int32 ncolpfix, nrowpun, vlenmnv;
149 mtx_matrix_t mtx;
150 int32 order, rank;
151 int32 c,len,vlen,r,rlen;
152 var_filter_t vf;
153 rel_filter_t rf;
154
155 /* CONSOLE_DEBUG("..."); */
156
157 rp = slv_get_solvers_rel_list(sys);
158 vp = slv_get_solvers_var_list(sys);
159 rlen = slv_get_num_solvers_rels(sys);
160 vlen = slv_get_num_solvers_vars(sys);
161 if (rlen ==0 || vlen == 0) return 1;
162 order = MAX(rlen,vlen);
163
164 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
165 rf.matchvalue = (REL_ACTIVE);
166 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
167 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
168
169 /* nrow plus unincluded rels and ineqs */
170 nrowpun = slv_count_solvers_rels(sys,&rf);
171 ncolpfix = slv_count_solvers_vars(sys,&vf); /* ncol plus fixed vars */
172
173 vf.matchbits = (VAR_SVAR);
174 vf.matchvalue = 0;
175 vlenmnv = vlen - slv_count_solvers_vars(sys,&vf); /* vlen minus non vars */
176
177 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
178 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
179 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
180 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
181
182 nrow = slv_count_solvers_rels(sys,&rf);
183 ncol = slv_count_solvers_vars(sys,&vf);
184
185 mtx = mtx_create();
186 mtx_set_order(mtx,order);
187
188 if (slv_make_incidence_mtx(sys,mtx,&vf,&rf)) {
189 ERROR_REPORTER_HERE(ASC_PROG_ERR,"failure in creating incidence matrix.");
190 mtx_destroy(mtx);
191 return 1;
192 }
193
194 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
195
196 mtx_output_assign(mtx,rlen,vlen);
197 rank = mtx_symbolic_rank(mtx);
198 if (rank == 0 ) return 1; /* nothing to do, eh? */
199
200 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
201
202 /* lot of whining about dof */
203 if (rank < nrow) {
204 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is row rank deficient (%d dependent equations)",
205 nrow - rank);
206 }
207 if (rank < ncol) {
208 if ( nrow != rank) {
209 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is row rank deficient with %d excess columns.",
210 ncol - rank);
211 } else {
212 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System has %d degrees of freedom.", ncol - rank);
213 }
214 }
215 if (ncol == nrow) {
216 if (ncol != rank) {
217 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is (%d) square but rank deficient.",ncol);
218 } else {
219 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"System is (%d) square.",ncol);
220 }
221 }
222 if (uppertriangular) {
223 mtx_ut_partition(mtx);
224 } else {
225 mtx_partition(mtx);
226 }
227 /* copy the block list. there is at least one if rank >=1 as above */
228 len = mtx_number_of_blocks(mtx);
229 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
230 if (newblocks == NULL) {
231 mtx_destroy(mtx);
232 return 2;
233 }
234 for (c = 0 ; c < len; c++) {
235 mtx_block(mtx,c,&(newblocks[c]));
236 }
237 slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
238 d = slv_get_dofdata(sys);
239 d->structural_rank = rank;
240 d->n_rows = nrow;
241 d->n_cols = ncol;
242
243 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
244
245 /* the next two lines assume inequalities are un-included. */
246 #define MINE 1
247 #if MINE
248 d->n_fixed = ncolpfix - ncol;
249 d->n_unincluded = nrowpun - nrow;
250 #else
251 d->n_fixed = vlen - ncol;
252 d->n_unincluded = rlen - nrow;
253 #endif
254 d->reorder.partition = 1; /* yes */
255 d->reorder.basis_selection = 0; /* none yet */
256 d->reorder.block_reordering = 0; /* none */
257
258 #ifdef BLOCKPARTITION_DEBUG
259 fp = fopen("/tmp/sbp1.plot","w+");
260 if (fp !=NULL) {
261 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
262 fclose(fp);
263 }
264 #endif
265
266 #if MINE
267 len = vlenmnv - 1; /* row of last possible solver variable */
268 for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
269 r = mtx_col_to_org(mtx,c);
270 if ( var_fixed(vp[r]) && var_active(vp[r])) {
271 c++;
272 } else {
273 mtx_swap_cols(mtx,c,len);
274 len--;
275 }
276 }
277 #endif
278
279 /* Now, leftover columns are where we want them.
280 * That is, unassigned incidence is next to assigned and fixed on right with
281 * no columns which do not correspond to variables in the range 0..vlen-1.
282 */
283 /* there aren't any nonincident vars in the solvers varlist, so
284 * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */
285
286 /* rows 0->rlen-1 should be actual rels, though, and not phantoms.
287 * At this point we should have
288 * rows - assigned
289 * - unassigned w/incidence
290 * - unincluded and included w/out free incidence and inactive mixed.
291 */
292 /* sort unincluded relations to bottom of good region */
293 /* this is needed because the idiot user may have fixed everything */
294 /* in a relation, but we want those rows not to be confused with */
295 /* unincluded rows. we're sort of dropping inequalities on the floor.*/
296 len = rlen - 1; /* row of last possible relation */
297 /* sort the inactive out of the unassigned/unincluded */
298 for (c=rank; len > c ; ) {
299 r = mtx_row_to_org(mtx,c);
300 if ( rel_active(rp[r])) {
301 c++;
302 } else {
303 mtx_swap_rows(mtx,c,len);
304 len--;
305 }
306 }
307
308 #if MINE
309 /* sort the unincluded out of the unassigned */
310 len = nrowpun - 1; /* row of last possible unassigned/unincluded relation */
311 for (c=rank; len > c ; ) {
312 r = mtx_row_to_org(mtx,c);
313 if ( rel_included(rp[r]) ) {
314 c++;
315 } else {
316 mtx_swap_rows(mtx,c,len);
317 len--;
318 }
319 }
320 #endif
321
322
323 #ifdef BLOCKPARTITION_DEBUG
324 fp = fopen("/tmp/sbp2.plot","w+");
325 if (fp !=NULL) {
326 mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
327 fclose(fp);
328 }
329 #endif
330 /* now, at last we have cols jacobian in the order we want the lists to
331 * be handed to the solvers. So, we need to reset the var sindex values
332 * and reorder the lists.
333 */
334 if (reindex_vars_from_mtx(sys,0,vlen-1,mtx)) {
335 mtx_destroy(mtx);
336 return 2;
337 }
338 /* now, at last we have rows jacobian in the order we want the lists to
339 * be handed to the solvers. So, we need to reset the rel sindex values.
340 */
341 if (reindex_rels_from_mtx(sys,0,rlen-1,mtx)) {
342 mtx_destroy(mtx);
343 return 2;
344 }
345
346 /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */
347
348 mtx_destroy(mtx);
349 return 0;
350 }
351
352
353 #if 0 /* code not currently used */
354 #\ifdef STATIC_HARWELL /* added the backslash so that syntax highlighting behaves */
355 /* note that you can't statically link to Harwell routines under the terms of the GPL */
356 extern void mc21b();
357 extern void mc13emod();
358 /* returns 0 if ok, OTHERWISE if madness detected.
359 * doesn't grok inequalities.
360 * CURRENTLY DOESN'T DO WELL WHEN NCOL<NROW
361 */
362 #define BLOCKPARTITION_DEBUG 0
363 int slv_block_partition_harwell(slv_system_t sys)
364 {
365 #\if BLOCKPARTITION_DEBUG
366 FILE *fp;
367 #\endif
368 struct rel_relation **rp;
369 struct rel_relation **rtmp;
370 struct rel_relation *rel;
371 struct var_variable **vp;
372 struct var_variable **vtmp;
373 struct var_variable *var;
374 const struct var_variable **list;
375 mtx_region_t *newblocks;
376 dof_t *d;
377 int32 nrow,ncol,size;
378 int32 order, rank;
379 int32 c,len,vlen,r,v,rlen,col,row;
380 int32 licn,rel_tmp_end,rel_count,var_perm_count,var_tmp_end,row_perm_count;
381 int32 col_incidence,numnz,row_high,col_high,dummy_ptr,row_ptr,num;
382 int32 *dummy_array_ptr;
383 var_filter_t vf;
384 rel_filter_t rf;
385 int32 *icn, *ip, *lenr, *iperm, *iw1, *iw2, *iw3, *arp, *ipnew, *ib;
386 int32 *row_perm;
387
388 /* get rel and var info */
389 rp = slv_get_solvers_rel_list(sys);
390 vp = slv_get_solvers_var_list(sys);
391 rlen = slv_get_num_solvers_rels(sys);
392 vlen = slv_get_num_solvers_vars(sys);
393 if (rlen ==0 || vlen == 0) return 1;
394 order = MAX(rlen,vlen);
395
396 /* allocate temp arrays */
397 vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
398 var = (struct var_variable *)ascmalloc(sizeof(struct var_variable *));
399 if (vtmp == NULL || var ==NULL) {
400 return 1;
401 }
402 rtmp = ASC_NEW_ARRAY(struct rel_relation *,rlen);
403 rel = ASC_NEW(struct rel_relation);
404 if (rtmp == NULL || rel ==NULL) {
405 return 1;
406 }
407
408 /* set up filters */
409 rf.matchbits = (REL_INCLUDED | REL_EQUALITY);
410 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY);
411 vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
412 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR);
413
414 /* count rows and cols */
415 nrow = slv_count_solvers_rels(sys,&rf);
416 ncol = slv_count_solvers_vars(sys,&vf);
417
418
419 /* count incidence for equality relations that are included
420 and active. Sorts out unincluded and inactive rels to
421 the end of the temporary relation list */
422 licn = 0;
423 rel_tmp_end = rlen -1;
424 for (r=0; r < rlen; r++) {
425 rel = rp[r];
426 if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
427 len = rel_n_incidences(rel);
428 list = rel_incidence_list(rel);
429 for (c=0; c < len; c++) {
430 if( var_apply_filter(list[c],&vf) ) {
431 licn++;
432 }
433 }
434 } else {
435 rtmp[rel_tmp_end] = rp[r];
436 rel_tmp_end--;
437 }
438 }
439
440 /* sort the inactive out of the unincluded and move to end */
441 len = rlen -1;
442 for(c = rel_tmp_end + 1; len > c; ) {
443 rel = rtmp[c];
444 if (rel_active(rel)) {
445 c++;
446 } else {
447 rtmp[c] = rtmp[len];
448 rtmp[len] = rel;
449 len--;
450 }
451 }
452
453 /* Sort variable list */
454 var_perm_count = 0;
455 var_tmp_end = vlen - 1;
456 for (c = 0; c < vlen; c++) {
457 var = vp[c];
458 if (var_apply_filter(var,&vf)) {
459 vtmp[var_perm_count] = var;
460 var_perm_count++;
461 } else {
462 vtmp[var_tmp_end] = var;
463 var_tmp_end--;
464 }
465 }
466 /* sort the inactive out of the unincluded and move to end */
467 len = vlen -1;
468 for(c = var_tmp_end + 1; len > c; ) {
469 var = vtmp[c];
470 if (var_active(var) && var_active(var)) {
471 c++;
472 } else {
473 vtmp[c] = vtmp[len];
474 vtmp[len] = var;
475 len--;
476 }
477 }
478
479 /* reset solver indicies to make life easier */
480 for (c = 0; c < vlen; c++) {
481 vp[c] = vtmp[c];
482 var_set_sindex(vp[c],c);
483 }
484
485 size = MAX(nrow,ncol);
486 /* Create vectors for fortran calls */
487 icn = ASC_NEW_ARRAY(int32,licn);
488 ip = ASC_NEW_ARRAY(int32,size);
489 lenr = ASC_NEW_ARRAY(int32,size);
490 iperm = ASC_NEW_ARRAY(int32,size);
491 iw1 = ASC_NEW_ARRAY(int32,size);
492 iw2 = ASC_NEW_ARRAY(int32,size);
493 iw3 = ASC_NEW_ARRAY(int32,size);
494 arp = ASC_NEW_ARRAY(int32,size);
495 ipnew = ASC_NEW_ARRAY(int32,size);
496 ib = ASC_NEW_ARRAY(int32,size);
497 row_perm = ASC_NEW_ARRAY(int32,size);
498
499
500 /* Fill incidence vectors and place included relations w/o
501 * incidence in temporary relation list before the unincluded
502 * and inactive relations
503 */
504 col_incidence = 0;
505 row_perm_count = 0;
506 for (r=0; r < rlen; r++) {
507 rel = rp[r];
508 if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
509 len = rel_n_incidences(rel);
510 if (len > 0) {
511 ip[row_perm_count] = col_incidence + 1; /*FORTRAN*/
512 lenr[row_perm_count] = len;
513 row_perm[row_perm_count] = r;
514 row_perm_count++;
515 list = rel_incidence_list(rel);
516 for (c=0; c < len; c++) {
517 if( var_apply_filter(list[c],&vf) ) {
518 col = var_sindex(list[c]);
519 icn[col_incidence] = col + 1; /*FORTRAN*/
520 col_incidence++;
521 }
522 }
523 } else {
524 rtmp[rel_tmp_end] = rel;
525 rel_tmp_end--;
526 }
527 }
528 }
529
530 licn = col_incidence;
531 order = row_perm_count;
532 mc21b(&order,icn,&licn,ip,lenr,iperm,&numnz,iw1,arp,iw2,iw3);
533 if (order == numnz){
534 FPRINTF(stderr,"they are equal\n");
535 } else if(order == numnz + 1){
536 FPRINTF(stderr,"ADD ONE\n");
537 } else {
538 FPRINTF(stderr,"no relationship\n");
539 }
540
541 if (rank == 0 ) {
542 return 1; /* nothing to do, eh? */
543 }
544 row_high = nrow - 1;
545 col_high = ncol -1;
546
547 /* lot of whining about dof */
548 if (rank < nrow) {
549 FPRINTF(stderr,"System is row rank deficient (dependent equations)\n");
550 row_high = rank - 1;
551 /* KHACK: have found that this is executed erroneously */
552 }
553 if (rank < ncol) {
554 if ( nrow != rank) {
555 FPRINTF(stderr,"System is row rank deficient with excess columns.\n");
556 } else {
557 FPRINTF(stderr,"System has degrees of freedom.\n");
558 }
559 }
560 if (ncol == nrow) {
561 FPRINTF(stderr,"System is (%d) square",ncol);
562 if (ncol != rank) {
563 FPRINTF(stderr,"but rank deficient.\n");
564 } else {
565 FPRINTF(stderr,".\n");
566 }
567 }
568
569 dummy_array_ptr = arp;
570 arp = lenr;
571 lenr = dummy_array_ptr;
572
573 /* Fix row pointers for assigned relations. Put unassigned
574 onto temporary relation list (before included w/o incidence)
575 and put dummy dense rows into ipnew */
576 dummy_ptr = licn + 1; /* FORTRAN */
577 for (row = 0; row < order; row++) {
578 row_ptr = iperm[row] - 1; /* FORTRAN */
579 if (row_ptr != -1) { /* indicates unassigned row */
580 ipnew[row_ptr] = ip[row];
581 lenr[row_ptr] = arp[row];
582 } else {
583 ipnew[row_ptr] = dummy_ptr;
584 dummy_ptr += ncol;
585 lenr[row_ptr] = ncol;
586 /* unassigned before fixed*/
587 rtmp[rel_tmp_end] = rp[row_perm[row_ptr]];
588 rel_tmp_end--;
589 }
590 }
591 for (row = order; row < ncol; row++) {
592 ipnew[row] = dummy_ptr;
593 dummy_ptr += ncol;
594 lenr[row] = ncol;
595 }
596
597 mc13emod(&size,icn,&licn,ipnew,lenr,arp,ib,&num,iw1,iw2,iw3);
598
599
600 /* copy the block list. there is at least one if rank >=1 as above */
601 len = num;
602 newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
603 if (newblocks == NULL) {
604 return 2;
605 }
606 /* set up locations of block corners with the row and column
607 orderings which will be set soon */
608 for (c = 0 ; c < len; c++) {
609 newblocks[c].row.low = ib[c]-1;
610 newblocks[c].col.low = newblocks[c].row.low;
611 }
612 for (c = 0 ; c < len -1; c++) {
613 newblocks[c].row.high = newblocks[c+1].row.low -1;
614 newblocks[c].col.high = newblocks[c].row.high;
615 }
616 newblocks[len - 1].row.high = row_high;
617 newblocks[len - 1].col.high = col_high;
618
619 slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
620 d = slv_get_dofdata(sys);
621 d->structural_rank = rank;
622 d->n_rows = nrow;
623 d->n_cols = ncol;
624 /* the next two lines assume only free and fixed incident solvervar
625 * appear in the solvers var list and inequalities are unincluded.
626 */
627 d->n_fixed = vlen - ncol;
628 d->n_unincluded = rlen - nrow;
629 d->reorder.partition = 1; /* yes */
630 d->reorder.basis_selection = 0; /* none yet */
631 d->reorder.block_reordering = 0; /* none */
632
633 for (c = 0; c < ncol; c++) {
634 v = arp[c] -1; /* FORTRAN */
635 vtmp[c] = vp[v];
636 }
637 for (c = 0; c < vlen; c++) {
638 vp[c] = vtmp[c];
639 var_set_sindex(vp[c],c);
640 }
641
642 rel_count = 0;
643 for (c = 0; c < size; c++) {
644 r = arp[c] - 1;
645 if (r < order){ /* weed out fake rows */
646 r = iperm[r] - 1;
647 if (r != -1) { /* unassigned already in rtmp */
648 rtmp[rel_count] = rp[row_perm[r]];
649 rel_count++;
650 }
651 }
652 }
653 for (c = 0; c < rlen; c++) {
654 rp[c] = rtmp[c];
655 rel_set_sindex(rp[c],c);
656 }
657 ascfree(vtmp);
658 ascfree(rtmp);
659 ascfree(icn);
660 ascfree(ip);
661 ascfree(lenr);
662 ascfree(iperm);
663 ascfree(iw1);
664 ascfree(iw2);
665 ascfree(iw3);
666 ascfree(arp);
667 ascfree(ipnew);
668 ascfree(ib);
669 return 0;
670 }
671 #\endif /*STATIC_HARWELL*/
672 #endif /* 0 */
673
674
675 int slv_block_unify(slv_system_t sys)
676 {
677 dof_t *d;
678 const mtx_block_t *mbt;
679 mtx_region_t *newblocks;
680
681 if (sys==NULL) return 1;
682
683 d = slv_get_dofdata(sys);
684 mbt = slv_get_solvers_blocks(sys);
685 assert(d!=NULL && mbt!=NULL);
686 if (d->structural_rank && mbt->nblocks >1) {
687 newblocks = ASC_NEW(mtx_region_t);
688 if (newblocks == NULL) return 2;
689 newblocks->row.low = newblocks->col.low = 0;
690 newblocks->row.high = d->structural_rank - 1;
691 newblocks->col.high = d->n_cols - 1;
692 if ( newblocks->col.high < 0 || newblocks->row.high < 0) {
693 ascfree(newblocks);
694 return 1;
695 }
696 slv_set_solvers_blocks(sys,1,newblocks);
697 }
698 return 0;
699 }
700
701 int slv_set_up_block(slv_system_t sys,int bnum)
702 {
703 struct rel_relation **rp;
704 struct var_variable **vp;
705 mtx_region_t reg;
706 const mtx_block_t *b;
707 int32 c,vlen,rlen;
708 var_filter_t vf;
709 rel_filter_t rf;
710
711 if (sys==NULL) return 1;
712 rlen = slv_get_num_solvers_rels(sys);
713 vlen = slv_get_num_solvers_vars(sys);
714 if (rlen ==0 || vlen == 0) return 1;
715
716 rp = slv_get_solvers_rel_list(sys);
717 vp = slv_get_solvers_var_list(sys);
718 assert(rp!=NULL);
719 assert(vp!=NULL);
720
721 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
722 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
723 vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
724 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
725
726 b = slv_get_solvers_blocks(sys);
727 assert(b!=NULL);
728 if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
729 reg = b->block[bnum];
730 for (c=reg.col.low; c<=reg.col.high; c++) {
731 var_set_in_block(vp[c],1);
732 }
733 for (c=reg.row.low; c<=reg.row.high; c++) {
734 rel_set_in_block(rp[c],1);
735 }
736 return 0;
737 }
738
739 /* returns 0 if ok, OTHERWISE if madness detected.
740 */
741 int slv_spk1_reorder_block(slv_system_t sys,int bnum,int transpose)
742 {
743 struct rel_relation **rp;
744 struct var_variable **vp;
745 mtx_region_t reg;
746 const mtx_block_t *b;
747 mtx_matrix_t mtx;
748 mtx_coord_t coord;
749 int32 c,vlen,rlen;
750 var_filter_t vf;
751 rel_filter_t rf;
752 dof_t *d;
753
754 if (sys==NULL) return 1;
755 rlen = slv_get_num_solvers_rels(sys);
756 vlen = slv_get_num_solvers_vars(sys);
757 if (rlen ==0 || vlen == 0) return 1;
758
759 rp = slv_get_solvers_rel_list(sys);
760 vp = slv_get_solvers_var_list(sys);
761 assert(rp!=NULL);
762 assert(vp!=NULL);
763 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
764 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
765 vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
766 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
767
768 mtx = mtx_create();
769 mtx_set_order(mtx,MAX(rlen,vlen));
770 b = slv_get_solvers_blocks(sys);
771 assert(b!=NULL);
772 if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
773 reg = b->block[bnum];
774 for (c=reg.col.low; c<=reg.col.high; c++) {
775 var_set_in_block(vp[c],1);
776 }
777 for (c=reg.row.low; c<=reg.row.high; c++) {
778 rel_set_in_block(rp[c],1);
779 }
780 if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
781 return 1; /* must be square */
782 /* could also enforce minsize 3x3, but someone my call with a
783 * partitionable region, so don't want to.
784 */
785 }
786
787 if (slv_make_incidence_mtx(sys,mtx,&vf,&rf)) {
788 FPRINTF(stderr,
789 "slv_spk1_reorder_block: failure in creating incidence matrix.\n");
790 mtx_destroy(mtx);
791 return 1;
792 }
793 /* verify that block has no empty columns, though not checking diagonal */
794 for (c = reg.row.low; c <= reg.row.high; c++) {
795 coord.col = mtx_FIRST;
796 coord.row = c;
797 if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
798 mtx_destroy(mtx);
799 FPRINTF(stderr, "slv_spk1_reorder_block: empty row (%d) found.\n",c);
800 return 1;
801 }
802 coord.row = mtx_FIRST;
803 coord.col = c;
804 if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
805 FPRINTF(stderr, "slv_spk1_reorder_block: empty col (%d) found.\n",c);
806 mtx_destroy(mtx);
807 return 1;
808 }
809 }
810 if (transpose) {
811 mtx_reorder(mtx,&reg,mtx_TSPK1);
812 } else {
813 mtx_reorder(mtx,&reg,mtx_SPK1);
814 }
815 if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
816 mtx_destroy(mtx);
817 return 2;
818 }
819 if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
820 mtx_destroy(mtx);
821 return 2;
822 }
823 d = slv_get_dofdata(sys);
824 d->reorder.block_reordering = 1; /* spk1 */
825
826 mtx_destroy(mtx);
827 return 0;
828 }
829
830 /*
831 A function to be called on the leftover diagonal blocks of mr_bisect.
832 This function should be thoroughly investigated, which it has not been.
833 */
834 static int tear_subreorder(slv_system_t server,
835 mtx_matrix_t mtx, mtx_region_t *reg)
836 {
837 assert(mtx != NULL && reg !=NULL);
838 (void)server;
839 if (g_blockmethod==mtx_UNKNOWN) {
840 mtx_reorder(mtx,reg,mtx_SPK1);
841 } else {
842 mtx_reorder(mtx,reg,g_blockmethod);
843 }
844 return 0;
845 }
846
847 int slv_tear_drop_reorder_block(slv_system_t sys, int32 bnum,
848 int32 cutoff, int two,
849 enum mtx_reorder_method blockmethod
850 ){
851 struct rel_relation **rp;
852 struct var_variable **vp;
853 const mtx_block_t *b;
854 dof_t *d;
855 mr_reorder_t *mrsys;
856 mtx_matrix_t mtx;
857 mtx_region_t reg;
858 mtx_coord_t coord;
859 int32 c,vlen,rlen,modcount;
860 var_filter_t vf;
861 rel_filter_t rf;
862
863 if (sys==NULL) return 1;
864 rlen = slv_get_num_solvers_rels(sys);
865 vlen = slv_get_num_solvers_vars(sys);
866 if (rlen ==0 || vlen == 0) return 1;
867
868 rp = slv_get_solvers_rel_list(sys);
869 vp = slv_get_solvers_var_list(sys);
870 assert(rp!=NULL);
871 assert(vp!=NULL);
872 rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
873 rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK| REL_ACTIVE);
874 vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
875 vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);
876
877 mtx = mtx_create();
878 mtx_set_order(mtx,MAX(rlen,vlen));
879 b = slv_get_solvers_blocks(sys);
880 assert(b!=NULL); /* probably shouldn't be an assert here ... */
881 if (bnum <0 || bnum >= b->nblocks || b->block == NULL) {
882 mtx_destroy(mtx);
883 return 1;
884 }
885 reg = b->block[bnum];
886 /* do the flag setting even if we don't do the reorder */
887 for (c=reg.col.low; c<=reg.col.high; c++) {
888 var_set_in_block(vp[c],1);
889 }
890 for (c=reg.row.low; c<=reg.row.high; c++) {
891 rel_set_in_block(rp[c],1);
892 }
893 if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
894 mtx_destroy(mtx);
895 return 1; /* must be square */
896 }
897 if (reg.row.high - reg.row.low < 3) {
898 mtx_destroy(mtx);
899 return 0; /* must be 3x3 or bigger to have any effect */
900 }
901
902 if (slv_make_incidence_mtx(sys,mtx,&vf,&rf)) {
903 FPRINTF(stderr,
904 "slv_tear_drop_reorder_block: failure in creating incidence matrix.\n");
905 mtx_destroy(mtx);
906 return 1;
907 }
908 /* verify that block has no empty columns, though not checking diagonal */
909 for (c = reg.row.low; c <= reg.row.high; c++) {
910 coord.col = mtx_FIRST;
911 coord.row = c;
912 if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
913 mtx_destroy(mtx);
914 FPRINTF(stderr, "slv_tear_drop_reorder_block: empty row (%d) found.\n",c);
915 return 1;
916 }
917 coord.row = mtx_FIRST;
918 coord.col = c;
919 if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
920 FPRINTF(stderr, "slv_tear_drop_reorder_block: empty col (%d) found.\n",c);
921 mtx_destroy(mtx);
922 return 1;
923 }
924 }
925
926 modcount = slv_get_num_models(sys);
927 mrsys = mr_reorder_create(sys,mtx,modcount);
928 if (mrsys==NULL) return 1;
929 mrsys->cutoff = cutoff;
930 g_blockmethod = blockmethod;
931 if (!two) {
932 mr_bisect_partition(mrsys,&reg,1,tear_subreorder);
933 } else {
934 mr_bisect_partition2(mrsys,&reg,1,tear_subreorder);
935 }
936 #if 1
937 FPRINTF(stderr,"Model-based reordering: (tear-style = %d)\n",two);
938 FPRINTF(stderr,"Min. tearable size:\t%d\n",mrsys->cutoff);
939 FPRINTF(stderr,"Tears:\t\t%d\n",mrsys->ntears);
940 FPRINTF(stderr,"Partitionings:\t%d\n",mrsys->nblts);
941 #endif
942 reg = b->block[bnum]; /* bisect likely munged reg */
943 if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
944 mtx_destroy(mtx);
945 mr_reorder_destroy(mrsys);
946 return 2;
947 }
948 if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
949 mtx_destroy(mtx);
950 mr_reorder_destroy(mrsys);
951 return 2;
952 }
953 d = slv_get_dofdata(sys);
954 d->reorder.block_reordering = 2; /* tear_drop_baa */
955
956 mtx_destroy(mtx);
957 mr_reorder_destroy(mrsys);
958 return 0;
959 }
960
961 /*------------------------------------------------------------------------------
962 DEBUG OUTPUT for BLOCK STRUCTURE
963
964 (it *would* belong in the 'mtx' section, but it outputs var and rel names
965 instead of just numbers *)
966 */
967
968 extern int system_block_debug(slv_system_t sys, FILE *fp){
969 int i,j,nr,nc;
970 dof_t *dof;
971 char *relname, *varname;
972 struct var_variable **vlist;
973 struct rel_relation **rlist;
974 mtx_region_t b;
975 dof = slv_get_dofdata(sys);
976 char s[80];
977 char color;
978
979 fprintf(fp,"\n\nSLV_SYSTEM BLOCK INFO\n\n");
980
981 fprintf(fp,"Structural rank: %d\n",dof->structural_rank);
982 fprintf(fp,"Included rels: %d\n",dof->n_rows);
983 fprintf(fp,"Incident, free vars: %d\n",dof->n_cols);
984 fprintf(fp,"Fixed vars: %d\n",dof->n_fixed);
985 fprintf(fp,"Unincluded rels: %d\n",dof->n_unincluded);
986 fprintf(fp,"Number of blocks: %d\n",dof->blocks.nblocks);
987
988 vlist = slv_get_solvers_var_list(sys);
989 rlist = slv_get_solvers_rel_list(sys);
990 color = (fp == stderr || fp==stdout);
991 for(i=0;i<dof->blocks.nblocks;++i){
992 if(color){
993 if(i%2)color_on(fp,ASC_FG_BROWN);
994 else color_on(fp,ASC_FG_BLACK);
995 }
996 b = dof->blocks.block[i];
997 nr = b.row.high - b.row.low + 1;
998 nc = b.col.high - b.col.low + 1;
999 SNPRINTF(s,80,"BLOCK %d (%d x %d)",i,nr,nc);
1000 fprintf(fp,"%-18s",s);
1001 SNPRINTF(s,80,"%-18s","");
1002 for(j=0;j<MAX(nr,nc); ++j){
1003 fprintf(fp,"%s%d",(j?s:""),j);
1004 if(j<nr){
1005 relname = rel_make_name(sys,rlist[b.row.low + j]);
1006 fprintf(fp,"\t%-20s",relname);
1007 ASC_FREE(relname);
1008 }else{
1009 fprintf(fp,"\t%-20s","");
1010 }
1011 if(j<nc){
1012 varname = var_make_name(sys,vlist[b.col.low + j]);
1013 fprintf(fp,"\t%-20s",varname);
1014 ASC_FREE(varname);
1015 }
1016 fprintf(fp,"\n");
1017 }
1018 }
1019 if(color)color_off(fp);
1020 return 0;
1021 }
1022
1023 extern int system_log_block_debug(slv_system_t sys, FILE *fp){
1024 int i,j,nr,nc;
1025 dof_t *dof;
1026 char *relname, *varname;
1027 struct dis_discrete **vlist;
1028 struct logrel_relation **rlist;
1029 mtx_region_t b;
1030 dof = slv_get_log_dofdata(sys);
1031 char s[80];
1032 char color;
1033
1034 fprintf(fp,"\n\nSLV_SYSTEM BLOCK INFO\n\n");
1035
1036 fprintf(fp,"Structural rank: %d\n",dof->structural_rank);
1037 fprintf(fp,"Included rels: %d\n",dof->n_rows);
1038 fprintf(fp,"Incident, free vars: %d\n",dof->n_cols);
1039 fprintf(fp,"Fixed vars: %d\n",dof->n_fixed);
1040 fprintf(fp,"Unincluded rels: %d\n",dof->n_unincluded);
1041 fprintf(fp,"Number of blocks: %d\n",dof->blocks.nblocks);
1042
1043 vlist = slv_get_solvers_dvar_list(sys);
1044 rlist = slv_get_solvers_logrel_list(sys);
1045 color = (fp == stderr || fp==stdout);
1046 for(i=0;i<dof->blocks.nblocks;++i){
1047 if(color){
1048 if(i%2)color_on(fp,ASC_FG_BROWN);
1049 else color_on(fp,ASC_FG_BLACK);
1050 }
1051 b = dof->blocks.block[i];
1052 nr = b.row.high - b.row.low + 1;
1053 nc = b.col.high - b.col.low + 1;
1054 SNPRINTF(s,80,"BLOCK %d (%d x %d)",i,nr,nc);
1055 fprintf(fp,"%-18s",s);
1056 SNPRINTF(s,80,"%-18s","");
1057 for(j=0;j<MAX(nr,nc); ++j){
1058 fprintf(fp,"%s%d",(j?s:""),j);
1059 if(j<nr){
1060 relname = logrel_make_name(sys,rlist[b.row.low + j]);
1061 fprintf(fp,"\t%-20s",relname);
1062 ASC_FREE(relname);
1063 }else{
1064 fprintf(fp,"\t%-20s","");
1065 }
1066 if(j<nc){
1067 varname = dis_make_name(sys,vlist[b.col.low + j]);
1068 fprintf(fp,"\t%-20s",varname);
1069 ASC_FREE(varname);
1070 }
1071 fprintf(fp,"\n");
1072 }
1073 }
1074 if(color)color_off(fp);
1075 return 0;
1076 }
1077
1078 /*------------------------------------------------------------------------------
1079 PARTITIONING for DIFFERENTIAL/ALGEBRAIC SYSTEMS
1080 */
1081
1082 /**
1083 This macro will generate 'system_var_list_debug(sys)' and 'system_rel_list_debug(sys)'
1084 and maybe other useful things if you're lucky.
1085 */
1086 #define LIST_DEBUG(TYPE,FULLTYPE) \
1087 void system_##TYPE##_list_debug(slv_system_t sys){ \
1088 struct FULLTYPE **list; \
1089 int n, i; \
1090 n = slv_get_num_solvers_##TYPE##s(sys); \
1091 list = slv_get_solvers_##TYPE##_list(sys); \
1092 \
1093 CONSOLE_DEBUG("printing " #TYPE " list (%d)",n); \
1094 char *name; \
1095 for(i=0;i<n;++i){ \
1096 name = TYPE##_make_name(sys,list[i]); \
1097 fprintf(stderr,"%d: %s\n",i,name); \
1098 ASC_FREE(name); \
1099 } \
1100 fprintf(stderr,"-----\n\n"); \
1101 }
1102
1103 LIST_DEBUG(var,var_variable)
1104 LIST_DEBUG(rel,rel_relation)
1105
1106 #ifdef CUT_DEBUG
1107 /* if you need this to work with non-GCC, take a page out of the utilities/error.h book */
1108 # define MAYBE_WRITE_LIST(TYPE) system_##TYPE##_list_debug(sys)
1109 # define MAYBE_CONSOLE_DEBUG(MSG,ARGS...) CONSOLE_DEBUG(MSG,##ARGS)
1110 #else
1111 # define MAYBE_WRITE_LIST(TYPE)
1112 # define MAYBE_CONSOLE_DEBUG(MSG,...)
1113 #endif
1114
1115 /* we define system_cut_vars and system_cut_rels with this macro... */
1116 /**
1117 This is a big durtie macro to perform cuts on our solvers_*_lists.
1118 The function will start at position 'begin' and move through all elements
1119 of the list from that point to the end. Any 'good' items matching the filter
1120 are moved to the start of the range traversed. And 'bad' ones that dont
1121 get moved to the end. At the end, the number of items found matching the
1122 filter is returned in 'numgood'. There will be that many 'good' items
1123 in place from position 'begin' onwards.
1124 */
1125 #define SYSTEM_CUT_LIST(TYPE,FULLTYPE) \
1126 int system_cut_##TYPE##s(slv_system_t sys, const int begin, const TYPE##_filter_t *filt, int *numgood){ \
1127 struct FULLTYPE **list, **start, **end, *swap; \
1128 int len, i; \
1129 char *name; \
1130 \
1131 asc_assert(filt); \
1132 \
1133 MAYBE_CONSOLE_DEBUG("CUTTING " #TYPE " LIST (filter=0x%8x,mask=0x%8x)...",filt->matchvalue,filt->matchbits); \
1134 MAYBE_WRITE_LIST(TYPE); \
1135 \
1136 list = slv_get_solvers_##TYPE##_list(sys); \
1137 len = slv_get_num_solvers_##TYPE##s(sys); \
1138 if(len==0){ \
1139 ERROR_REPORTER_HERE(ASC_PROG_ERR,"There are no " #TYPE "s!"); \
1140 return 1; \
1141 } \
1142 \
1143 MAYBE_CONSOLE_DEBUG("SORTING"); \
1144 \
1145 start = list + begin; \
1146 end = list + len; \
1147 while(start < end){ \
1148 if(TYPE##_apply_filter(*start,filt)){ \
1149 start++; continue; \
1150 } \
1151 if(!TYPE##_apply_filter(*(--end),filt))continue; \
1152 swap = *end; \
1153 *end = *start; \
1154 *start = swap; \
1155 start++; \
1156 } \
1157 \
1158 MAYBE_WRITE_LIST(TYPE); \
1159 \
1160 MAYBE_CONSOLE_DEBUG("UPDATING"); \
1161 \
1162 /* update the sindex for each after start */ \
1163 *numgood = 0; \
1164 MAYBE_CONSOLE_DEBUG("numgood = %d",*numgood); \
1165 for(i=begin;i<len;++i){ \
1166 name = TYPE##_make_name(sys,list[i]); \
1167 if(TYPE##_apply_filter(list[i],filt)){ \
1168 MAYBE_CONSOLE_DEBUG("%s: good",name); \
1169 (*numgood)++; \
1170 }else{ \
1171 MAYBE_CONSOLE_DEBUG("'%s' will cut to back (flags=0x%8x)",name,var_flags(list[i])); \
1172 } \
1173 ASC_FREE(name); \
1174 TYPE##_set_sindex(list[i],i); \
1175 } \
1176 MAYBE_CONSOLE_DEBUG("numgood = %d",*numgood); \
1177 \
1178 return 0; \
1179 }
1180
1181 SYSTEM_CUT_LIST(var,var_variable);
1182 SYSTEM_CUT_LIST(rel,rel_relation);
1183

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