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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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