/[ascend]/trunk/base/generic/system/block.c
ViewVC logotype

Contents of /trunk/base/generic/system/block.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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