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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 52600 byte(s)
Setting up web subdirectory in repository
1 /* mtx: Ascend Sparse Matrix Package
2 * by Karl Michael Westerberg
3 * Created: 2/6/90
4 * Version: $Revision: 1.12 $
5 * Version control file: $RCSfile: mtx_perms.c,v $
6 * Date last modified: $Date: 1998/05/06 17:28:51 $
7 * Last modified by: $Author: ballan $
8 *
9 * This file is part of the SLV solver.
10 *
11 * Copyright (C) 1990 Karl Michael Westerberg
12 * Copyright (C) 1993 Joseph Zaher
13 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
14 * Copyright (C) 1995 Kirk Andre Abbott, Benjamin Andrew Allan
15 * Copyright (C) 1996 Benjamin Andrew Allan
16 *
17 * The SLV solver is free software; you can redistribute
18 * it and/or modify it under the terms of the GNU General Public License as
19 * published by the Free Software Foundation; either version 2 of the
20 * License, or (at your option) any later version.
21 *
22 * The SLV solver is distributed in hope that it will be
23 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25 * General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License along with
28 * the program; if not, write to the Free Software Foundation, Inc., 675
29 * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
30 * COPYING is found in ../compiler.
31 */
32 #include "utilities/ascConfig.h"
33 #include "utilities/ascMalloc.h"
34 #include "utilities/mem.h"
35 #include "solver/mtx.h"
36 /* grab our private parts */
37 #define __MTX_C_SEEN__
38 #include "solver/mtx_use_only.h"
39
40
41 static mtx_block_perm_t alloc_block_perm()
42 /**
43 *** Allocates a block_perm header and returns a pointer to it.
44 **/
45 {
46 mtx_block_perm_t bp;
47 bp = (mtx_block_perm_t)asccalloc(1,sizeof(struct mtx_block_perm_structure));
48 bp->integrity = OK;
49 bp->data =
50 (struct structural_data_t *)ascmalloc(sizeof(struct structural_data_t));
51 return(bp);
52 }
53 #define free_block_perm(bp) ascfree(bp)
54
55 /**
56 This is a strange little function which, at the expense of much array logic,
57 sets an org_to_cur permutation array based on the data in a cur_to_org.
58 Method:
59 for i = start..len-1
60 o2c[c2o[i]] = i;
61 If the majority of an org_to_cur is being copied, it is probably much
62 better to copy the whole permutation from the relevant original org_to_cur
63 rather than piecewise as is done here.
64 WARNING: ABSOLUTELY NOTHING is being done to insure sanity in the function.
65 If you can't draw the picture of what this function does by hand and
66 compute an example result, don't use it. That said, this saves lots of
67 otherwise very annoying time spent in re-reordering.
68 **/
69 static void fix_org2cur_from_cur2org(int32 *c2o, int32 start,
70 int32 *o2c, int32 len) {
71 int32 i;
72 if (ISNULL(c2o) || ISNULL(o2c) || start <0 || len < 1) {
73 FPRINTF(g_mtxerr,"ERROR (mtx-internal) fix_org2cur_from_cur2org");
74 FPRINTF(g_mtxerr," got bogus args");
75 return;
76 }
77 for (i=start; i < len ; i++) {
78 o2c[c2o[i]] = i;
79 }
80 return;
81 }
82
83 extern mtx_block_perm_t mtx_create_block_perm(mtx_matrix_t mtx)
84 {
85 mtx_block_perm_t bp = NULL;
86 int32 bnum;
87
88 if (ISSLAVE(mtx)) {
89 return (mtx_create_block_perm(mtx->master));
90 }
91 if (!mtx_check_matrix(mtx)) return NULL;
92 if (mtx->capacity < 2) {
93 FPRINTF(g_mtxerr,"ERROR: mtx_create_block_perm: mtx given is too\n");
94 FPRINTF(g_mtxerr," small to permute.\n");
95 FPRINTF(g_mtxerr," Returning NULL.\n");
96 return NULL;
97 }
98 if (mtx->data->symbolic_rank < 0) {
99 FPRINTF(g_mtxerr,"ERROR: mtx_create_block_perm: mtx given is not\n");
100 FPRINTF(g_mtxerr," output assigned.\n");
101 FPRINTF(g_mtxerr," Returning NULL.\n");
102 return NULL;
103 }
104
105 /* create same size matrix information */
106 bp = alloc_block_perm();
107 if (ISNULL(bp)) {
108 FPRINTF(g_mtxerr,"ERROR: mtx_create_block_perm: insufficient memory\n");
109 FPRINTF(g_mtxerr," allocating bperm (1).\n");
110 FPRINTF(g_mtxerr," Returning NULL.\n");
111 return NULL;
112 }
113 bp->order = mtx->order;
114 bp->capacity = mtx->capacity;
115 bp->mtx = mtx;
116
117 /* copy permutation information */
118 bp->perm.row.org_to_cur = mtx_alloc_perm(bp->capacity);
119 bp->perm.row.cur_to_org = mtx_alloc_perm(bp->capacity);
120 bp->perm.col.org_to_cur = mtx_alloc_perm(bp->capacity);
121 bp->perm.col.cur_to_org = mtx_alloc_perm(bp->capacity);
122 if (ISNULL(bp->perm.row.org_to_cur) ||
123 ISNULL(bp->perm.row.cur_to_org) ||
124 ISNULL(bp->perm.col.org_to_cur) ||
125 ISNULL(bp->perm.col.cur_to_org) ) {
126 FPRINTF(g_mtxerr,"ERROR: mtx_create_block_perm: insufficient memory\n");
127 FPRINTF(g_mtxerr," allocating bperm (2).\n");
128 FPRINTF(g_mtxerr," Returning NULL.\n");
129 mtx_free_perm(bp->perm.row.org_to_cur);
130 mtx_free_perm(bp->perm.row.cur_to_org);
131 mtx_free_perm(bp->perm.col.org_to_cur);
132 mtx_free_perm(bp->perm.col.cur_to_org);
133 free_block_perm(bp);
134 return NULL;
135 }
136
137 mtx_copy_perm(bp->perm.row.org_to_cur,
138 mtx->perm.row.org_to_cur,bp->capacity);
139 mtx_copy_perm(bp->perm.row.cur_to_org,
140 mtx->perm.row.cur_to_org,bp->capacity);
141 bp->perm.row.parity = mtx_row_parity(mtx);
142
143 mtx_copy_perm(bp->perm.col.org_to_cur,
144 mtx->perm.col.org_to_cur,bp->capacity);
145 mtx_copy_perm(bp->perm.col.cur_to_org,
146 mtx->perm.col.cur_to_org,bp->capacity);
147 bp->perm.col.parity = mtx_col_parity(mtx);
148
149 bp->data = (struct structural_data_t *)
150 ascmalloc(sizeof(struct structural_data_t));
151 bp->data->symbolic_rank = mtx->data->symbolic_rank;
152 bp->data->nblocks = mtx->data->nblocks;
153 bp->data->block = mtx->data->nblocks > 0 ? (mtx_region_t *)
154 ascmalloc( mtx->data->nblocks*sizeof(mtx_region_t) ) : NULL;
155 for( bnum=0; bnum < mtx->data->nblocks; bnum++ ) {
156 bp->data->block[bnum] = mtx->data->block[bnum];
157 }
158 return bp;
159 }
160
161 int mtx_update_block_perm(mtx_matrix_t mtx, int32 bnum,
162 mtx_block_perm_t bp)
163 {
164 if (ISSLAVE(mtx)) {
165 return (mtx_update_block_perm(mtx->master,bnum,bp));
166 }
167 return 3;
168 }
169
170 int mtx_restore_block_perm(mtx_matrix_t mtx, int32 bnum,
171 mtx_block_perm_t bp)
172 {
173 int32 len;
174 if (mtx_check_matrix(mtx) != 1) {
175 FPRINTF(g_mtxerr,"ERROR: (mtx) restore_block_perm: Called with bad\n");
176 FPRINTF(g_mtxerr," mtx.\n");
177 return 1;
178 }
179 if (ISSLAVE(mtx)) {
180 return(mtx_restore_block_perm(mtx,bnum,bp));
181 }
182 if (ISNULL(bp) || bp->integrity != OK) {
183 FPRINTF(g_mtxerr,"ERROR: (mtx) restore_block_perm: Called with bad\n");
184 FPRINTF(g_mtxerr," bperm.\n");
185 return 2;
186 }
187 if (bp->mtx != mtx) {
188 FPRINTF(g_mtxerr,"ERROR: (mtx) restore_block_perm: Mismatched mtx\n");
189 FPRINTF(g_mtxerr," and bperm given.\n");
190 return 3;
191 }
192 if (bp->order != mtx->order) {
193 FPRINTF(g_mtxerr,"ERROR: (mtx) restore_block_perm: Mismatched order mtx\n");
194 FPRINTF(g_mtxerr," and bperm given.\n");
195 return 3;
196 }
197 if (bnum != mtx_ALL_BLOCKS && bp->data->nblocks != 0) {
198 /* check block sanity */
199 if (mtx->data->symbolic_rank < 0) {
200 FPRINTF(g_mtxerr,"ERROR: mtx_restore_block_perm: mtx given is not\n");
201 FPRINTF(g_mtxerr," output assigned.\n");
202 FPRINTF(g_mtxerr," Returning NULL.\n");
203 return 3;
204 }
205 if ( bp->data->nblocks != mtx->data->nblocks ||
206 bp->data->block[bnum].row.low != mtx->data->block[bnum].row.low ||
207 bp->data->block[bnum].col.low != mtx->data->block[bnum].col.low ||
208 bp->data->block[bnum].row.high != mtx->data->block[bnum].row.high ||
209 bp->data->block[bnum].col.high != mtx->data->block[bnum].col.high ) {
210 FPRINTF(g_mtxerr,
211 "ERROR: (mtx) restore_block_perm: Mismatched block mtx\n");
212 FPRINTF(g_mtxerr,
213 " and bperm given.\n");
214 return 3;
215 }
216 /* now copy block bnumth info accordingly, resetting parity */
217 /* remember that copy_perm is for some weird reason cp(targ,src,len) */
218
219 len = bp->data->block[bnum].row.high - bp->data->block[bnum].row.low +1;
220
221 /* block given is in cur coordinates, so we can copy relevant stretches
222 of the cur_to_org permutations sequentially. pointer arithmetic. */
223 mtx_copy_perm((mtx->perm.row.cur_to_org) + bp->data->block[bnum].row.low,
224 (bp->perm.row.cur_to_org) + bp->data->block[bnum].row.low,
225 len);
226 fix_org2cur_from_cur2org(mtx->perm.row.cur_to_org,
227 bp->data->block[bnum].row.low,
228 mtx->perm.row.org_to_cur, len);
229 mtx->perm.row.parity = bp->perm.row.parity;
230
231 len = bp->data->block[bnum].col.high - bp->data->block[bnum].col.low +1;
232 mtx_copy_perm((mtx->perm.col.cur_to_org) + bp->data->block[bnum].col.low,
233 (bp->perm.col.cur_to_org) + bp->data->block[bnum].col.low,
234 len);
235 fix_org2cur_from_cur2org(mtx->perm.col.cur_to_org,
236 bp->data->block[bnum].col.low,
237 mtx->perm.col.org_to_cur, len);
238 mtx->perm.col.parity = bp->perm.col.parity;
239
240 } else {
241 /* copy all perm info */
242 mtx_copy_perm(mtx->perm.row.cur_to_org,
243 bp->perm.row.cur_to_org,bp->capacity);
244 mtx_copy_perm(mtx->perm.row.org_to_cur,
245 bp->perm.row.org_to_cur,bp->capacity);
246 mtx->perm.row.parity = bp->perm.row.parity;
247
248 mtx_copy_perm(mtx->perm.col.org_to_cur,
249 bp->perm.col.org_to_cur,bp->capacity);
250 mtx_copy_perm(mtx->perm.col.cur_to_org,
251 bp->perm.col.cur_to_org,bp->capacity);
252 mtx->perm.col.parity = bp->perm.col.parity;
253
254 /* copy all oa/block info if any. */
255 mtx->data->symbolic_rank = bp->data->symbolic_rank;
256 if ( mtx->data->nblocks != bp->data->nblocks) {
257 if (NOTNULL(mtx->data->block)) {
258 ascfree(mtx->data->block);
259 }
260 mtx->data->nblocks = bp->data->nblocks;
261 mtx->data->block = bp->data->nblocks > 0 ? (mtx_region_t *)
262 ascmalloc( mtx->data->nblocks*sizeof(mtx_region_t) ) : NULL;
263 }
264 for( bnum=0; bnum < bp->data->nblocks; bnum++ ) {
265 mtx->data->block[bnum] = bp->data->block[bnum];
266 }
267 }
268 #if MTX_DEBUG
269 if (super_check_matrix(mtx)) {
270 FPRINTF(g_mtxerr,"ERROR: mtx_restore_block_perm: mtx restored now\n");
271 FPRINTF(g_mtxerr," fails sanity.\n");
272 return 3;
273 }
274 #endif
275 return 0;
276 }
277
278 extern int mtx_destroy_block_perm(mtx_block_perm_t bp)
279 {
280 if (bp->integrity != OK) {
281 if (bp->integrity == DESTROYED) {
282 FPRINTF(g_mtxerr,
283 "ERROR: mtx_destroy_block_perm: Called with recently\n");
284 FPRINTF(g_mtxerr," destroyed bperm.\n");
285 } else {
286 FPRINTF(g_mtxerr,
287 "ERROR: mtx_destroy_block_perm: Called with apparently\n");
288 FPRINTF(g_mtxerr," corrupted bperm.\n");
289 }
290 return 1;
291 }
292
293 mtx_free_perm(bp->perm.row.org_to_cur);
294 mtx_free_perm(bp->perm.row.cur_to_org);
295 mtx_free_perm(bp->perm.col.org_to_cur);
296 mtx_free_perm(bp->perm.col.cur_to_org);
297
298 if( NOTNULL(bp->data->block) )
299 ascfree(bp->data->block);
300 ascfree(bp->data);
301
302 bp->integrity=DESTROYED;
303 free_block_perm(bp);
304 return 0;
305 }
306
307
308 extern size_t mtx_block_perm_size(mtx_block_perm_t bp)
309 {
310 size_t size=0;
311 if (ISNULL(bp)) return size;
312 size += sizeof(struct mtx_block_perm_structure);
313 /* data */
314 size += sizeof(struct structural_data_t);
315 /* block data */
316 if (bp->data->nblocks >0)
317 size += sizeof(mtx_region_t)*(size_t)(bp->data->nblocks-1);
318 /* permutations */
319 size += (size_t)4*(sizeof(int32)*(size_t)bp->capacity);
320 return size;
321 }
322
323 /* DO NOT CALL this on the perm of a matrix that is a slave */
324 static void swap( struct permutation_t *perm, int32 cur1,int32 cur2)
325 {
326 register org1,org2;
327
328 if( cur1 == cur2 ) return;
329 org1 = perm->cur_to_org[cur1];
330 org2 = perm->cur_to_org[cur2];
331 perm->cur_to_org[cur1] = org2;
332 perm->cur_to_org[cur2] = org1;
333 perm->org_to_cur[org1] = cur2;
334 perm->org_to_cur[org2] = cur1;
335 perm->parity = !(perm->parity);
336 }
337
338 void mtx_swap_rows( mtx_matrix_t mtx, int32 cur1,int32 cur2)
339 {
340 #if MTX_DEBUG
341 if(!mtx_check_matrix(mtx)) return;
342 #endif
343 if (ISSLAVE(mtx)) {
344 mtx = mtx->master;
345 #if MTX_DEBUG
346 if(!mtx_check_matrix(mtx)) return;
347 #endif
348 }
349 swap(&(mtx->perm.row),cur1,cur2);
350 }
351
352 void mtx_swap_cols( mtx_matrix_t mtx, int32 cur1,int32 cur2)
353 {
354 #if MTX_DEBUG
355 if(!mtx_check_matrix(mtx)) return;
356 #endif
357 if (ISSLAVE(mtx)) {
358 mtx = mtx->master;
359 #if MTX_DEBUG
360 if(!mtx_check_matrix(mtx)) return;
361 #endif
362 }
363 swap(&(mtx->perm.col),cur1,cur2);
364 }
365
366 void mtx_drag( mtx_matrix_t mtx, int32 begin, int32 end)
367 {
368 register int32 *rcto, *ccto, *rotc, *cotc;
369 int32 holdrow,holdcol;
370
371 #if MTX_DEBUG
372 if(!mtx_check_matrix(mtx)) return;
373 #endif
374 if (ISSLAVE(mtx)) {
375 mtx = mtx->master;
376 #if MTX_DEBUG
377 if(!mtx_check_matrix(mtx)) return;
378 #endif
379 }
380 if( begin == end ) return;
381 if( abs(begin-end) % 2 ) {
382 mtx->perm.row.parity = !(mtx->perm.row.parity);
383 mtx->perm.col.parity = !(mtx->perm.col.parity);
384 }
385 rcto = mtx->perm.row.cur_to_org;
386 ccto = mtx->perm.col.cur_to_org;
387 rotc = mtx->perm.row.org_to_cur;
388 cotc = mtx->perm.col.org_to_cur;
389 holdrow = rcto[begin];
390 holdcol = ccto[begin];
391 while( begin < end ) {
392 rotc[rcto[begin]=rcto[begin+1]] = begin;
393 cotc[ccto[begin]=ccto[begin+1]] = begin;
394 begin++;
395 }
396 while( begin > end ) {
397 rotc[rcto[begin]=rcto[begin-1]] = begin;
398 cotc[ccto[begin]=ccto[begin-1]] = begin;
399 begin--;
400 }
401 rotc[rcto[end]=holdrow] = end;
402 cotc[ccto[end]=holdcol] = end;
403 }
404
405 void mtx_reverse_diagonal( mtx_matrix_t mtx, int32 begin, int32 end)
406 {
407 register int32 center,j;
408 register struct permutation_t *rperm, *cperm;
409 #if MTX_DEBUG
410 if(!mtx_check_matrix(mtx)) return;
411 #endif
412 if (ISSLAVE(mtx)) {
413 mtx = mtx->master;
414 #if MTX_DEBUG
415 if(!mtx_check_matrix(mtx)) return;
416 #endif
417 }
418 if (begin > end) return;
419 if (end >= mtx->order) return;
420 center=(end-begin+1)/2;
421 rperm = &(mtx->perm.row);
422 cperm = &(mtx->perm.col);
423 for (j=0; j<center; j++) {
424 swap(cperm,begin,end);
425 swap(rperm,begin,end);
426 begin++;
427 end--;
428 }
429 }
430
431 int32 mtx_row_to_org( mtx_matrix_t mtx, int32 row)
432 {
433 #if MTX_DEBUG
434 if(!mtx_check_matrix(mtx)) return -1;
435 #endif
436 if (ISSLAVE(mtx)) {
437 mtx = mtx->master;
438 #if MTX_DEBUG
439 if(!mtx_check_matrix(mtx)) return -1;
440 #endif
441 }
442 return( mtx->perm.row.cur_to_org[row] );
443 }
444
445 int32 mtx_col_to_org(mtx,col)
446 mtx_matrix_t mtx;
447 int32 col;
448 {
449 #if MTX_DEBUG
450 if(!mtx_check_matrix(mtx)) return -1;
451 #endif
452 if (ISSLAVE(mtx)) {
453 mtx = mtx->master;
454 #if MTX_DEBUG
455 if(!mtx_check_matrix(mtx)) return -1;
456 #endif
457 }
458 return( mtx->perm.col.cur_to_org[col] );
459 }
460
461 int32 mtx_org_to_row( mtx_matrix_t mtx, int32 org)
462 {
463 #if MTX_DEBUG
464 if(!mtx_check_matrix(mtx)) return -1;
465 #endif
466 if (ISSLAVE(mtx)) {
467 mtx = mtx->master;
468 #if MTX_DEBUG
469 if(!mtx_check_matrix(mtx)) return -1;
470 #endif
471 }
472 return( mtx->perm.row.org_to_cur[org] );
473 }
474
475 int32 mtx_org_to_col( mtx_matrix_t mtx, int32 org)
476 {
477 #if MTX_DEBUG
478 if(!mtx_check_matrix(mtx)) return -1;
479 #endif
480 if (ISSLAVE(mtx)) {
481 mtx = mtx->master;
482 #if MTX_DEBUG
483 if(!mtx_check_matrix(mtx)) return -1;
484 #endif
485 }
486 return( mtx->perm.col.org_to_cur[org] );
487 }
488
489 boolean mtx_row_parity(mtx_matrix_t mtx)
490 {
491 #if MTX_DEBUG
492 if(!mtx_check_matrix(mtx)) return -1;
493 if (ISSLAVE(mtx)) {
494 if(!mtx_check_matrix(mtx->master)) return -1;
495 }
496 #endif
497 if (ISSLAVE(mtx)) {
498 return( mtx->master->perm.row.parity );
499 } else {
500 return( mtx->perm.row.parity );
501 }
502 }
503
504 boolean mtx_col_parity( mtx_matrix_t mtx)
505 {
506 #if MTX_DEBUG
507 if(!mtx_check_matrix(mtx)) return -1;
508 if (ISSLAVE(mtx)) {
509 if(!mtx_check_matrix(mtx->master)) return -1;
510 }
511 #endif
512 if (ISSLAVE(mtx)) {
513 return( mtx->master->perm.col.parity );
514 } else {
515 return( mtx->perm.col.parity );
516 }
517 }
518
519 /***********************************************************************\
520 structural analysis stuff
521 \***********************************************************************/
522 struct assign_row_vars {
523 int32 rv_indicator; /* rows visited indicator */
524 int32 *row_visited;
525 mtx_range_t unassigned_cols;
526 mtx_range_t assigned_cols;
527 };
528
529 /* MUST NOT BE CALLED ON slave matrix */
530 static int32 assign_row( mtx_matrix_t mtx, int32 row,
531 struct assign_row_vars *vars)
532 /**
533 *** Attempts to assign one of the columns in the range
534 *** vars->unassigned_columns to the given row (exchange that column with
535 *** column "row" in order to place a non-zero on the diagonal). Returns
536 *** the column that is assigned, or < 0 if not possible.
537 ***
538 *** Assignment is done row-wise by starting with the top row, finding the
539 *** nearest column in the row to the right of the diagonal to assign to it,
540 *** and advancing to the next row. If no column exists to the right of the
541 *** diagonal, then it precedes along a steward's path to find amongst the
542 *** previously assigned columns one which will be more suitably assigned to
543 *** it, thereby forcing a row which was visited earlier to be re-assigned to
544 *** some new column. Termination occurs when there are no more unvisited
545 *** rows. If there still remains non-empty yet unassigned rows, they will
546 *** be positioned just after the last assigned row. Same goes for any
547 *** remaining non-empty unassigned columns. The condensed grouping of
548 *** assigned rows and columns will be of dimension equal to the symbolic
549 *** rank of the matrix and are used to form an initial block.
550 **/
551 {
552 struct element_t Rewind, *elt;
553 int32 *tocur;
554
555 tocur = mtx->perm.col.org_to_cur;
556 Rewind.next.col = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
557 if( (elt=mtx_next_col(&Rewind,&(vars->unassigned_cols),tocur)) ) {
558 /* Cheap assignment */
559 register int32 col;
560 col = tocur[elt->col];
561 swap(&(mtx->perm.col),col,row);
562 return(col);
563 }
564
565 vars->row_visited[row] = vars->rv_indicator;
566 for( elt = &Rewind ; (elt=mtx_next_col(elt,&(vars->assigned_cols),tocur));) {
567 int32 col;
568 col = tocur[elt->col];
569 if( vars->row_visited[col] != vars->rv_indicator ) {
570 register int32 assigned_col;
571 if( (assigned_col = assign_row(mtx,col,vars)) >= ZERO ) {
572 swap(&(mtx->perm.col),assigned_col,row);
573 return( assigned_col );
574 }
575 }
576 }
577 return( -1 );
578 }
579
580 /*
581 * This function is very similar to mtx_output_assign. It however
582 * takes a region in which to do this output_assignement.
583 * It returns the number of rows that could be assigned.
584 * Remember that rv stands for row_visited, for newcomers to this file.
585 */
586 int mtx_output_assign_region(mtx_matrix_t mtx,
587 mtx_region_t *region,
588 int *orphaned_rows)
589 {
590 struct assign_row_vars vars;
591 struct element_t Rewind;
592 int32 nrows, *tocur;
593 mtx_range_t single_col;
594 mtx_coord_t nz;
595 int notassigned = 0;
596
597 #if MTX_DEBUG
598 if( !mtx_check_matrix(mtx) ) return 0;
599 #endif
600
601 if (ISSLAVE(mtx)) {
602 mtx = mtx->master;
603 #if MTX_DEBUG
604 if( !mtx_check_matrix(mtx) ) return 0;
605 #endif
606 }
607
608 vars.row_visited = mtx->order > ZERO ?
609 (int32 *)asccalloc(mtx->order, sizeof(int32)) : NULL;
610 vars.rv_indicator = 1;
611
612 if (region!=mtx_ENTIRE_MATRIX) {
613 vars.unassigned_cols.high = region->col.high;
614 vars.assigned_cols.low = region->col.low;
615 tocur = mtx->perm.col.org_to_cur;
616 nrows = region->row.high - region->row.low + 1;
617 }
618 else{
619 vars.unassigned_cols.high = mtx->order - 1;
620 vars.assigned_cols.low = ZERO;
621 tocur = mtx->perm.col.org_to_cur;
622 nrows = mtx->order;
623 }
624
625
626 for( nz.row = region->row.low ; nz.row < nrows ; ) {
627 single_col.low = single_col.high = nz.row;
628 vars.unassigned_cols.low = nz.row;
629 vars.assigned_cols.high = nz.row - 1;
630
631 /*
632 * Find next nnz in row and attempt to assign it using the
633 * recursive function assign_row. If successful, proceed
634 * to the next row. Otherwise, move the row down to end of the
635 * matrix. -- we are structurally singular.
636 */
637 Rewind.next.col = mtx->hdr.row[mtx->perm.row.cur_to_org[nz.row]];
638 if (mtx_next_col(&Rewind, &single_col, tocur) ||
639 assign_row(mtx, nz.row, &vars) >= ZERO) {
640 ++nz.row; /* Go to next row */
641 }
642 else {
643 swap(&(mtx->perm.row), nz.row, --nrows); /* Move row to the end */
644 notassigned++;
645 }
646 ++vars.rv_indicator; /* Effectively clear vars.row_visited */
647 }
648
649 if (NOTNULL(vars.row_visited)) ascfree(vars.row_visited);
650 *orphaned_rows = notassigned; /* return count of unassigned rows */
651 return nrows;
652 }
653
654 void mtx_output_assign( mtx_matrix_t mtx, int32 hirow, int32 hicol)
655 {
656 struct assign_row_vars vars;
657 struct element_t Rewind;
658 int32 nrows, ndxhigh, *tocur;
659 mtx_range_t single_col;
660 mtx_coord_t nz;
661
662 #if MTX_DEBUG
663 if( !mtx_check_matrix(mtx) ) return; /*ben*/
664 #endif
665 if (ISSLAVE(mtx)) {
666 mtx = mtx->master;
667 #if MTX_DEBUG
668 if( !mtx_check_matrix(mtx) ) return;
669 #endif
670 }
671 vars.row_visited = mtx->order > ZERO ?
672 (int32 *)asccalloc( mtx->order,sizeof(int32) ) : NULL;
673 vars.rv_indicator = 1;
674 vars.unassigned_cols.high = mtx->order-1;
675 vars.assigned_cols.low = ZERO;
676 tocur = mtx->perm.col.org_to_cur;
677
678 nrows = mtx->order;
679 for( nz.row = ZERO ; nz.row < nrows ; ) {
680 single_col.low = single_col.high = nz.row;
681 vars.unassigned_cols.low = nz.row;
682 vars.assigned_cols.high = nz.row-1;
683 Rewind.next.col = mtx->hdr.row[mtx->perm.row.cur_to_org[nz.row]];
684 if( mtx_next_col(&Rewind,&single_col,tocur) ||
685 assign_row(mtx,nz.row,&vars) >= ZERO ) {
686 ++nz.row; /* Go to next row */
687 } else {
688 swap(&(mtx->perm.row),nz.row,--nrows); /* Move row to the end */
689 }
690 ++vars.rv_indicator; /* Effectively clear vars.row_visited */
691 }
692
693 mtx->data->symbolic_rank = nrows;
694 if( NOTNULL(mtx->data->block) )
695 ascfree(mtx->data->block);
696
697 /**
698 *** Place all assigned rows and columns into one
699 *** initial block.
700 **/
701 if( mtx->data->symbolic_rank > 0 ) {
702 mtx->data->nblocks = 1;
703 mtx->data->block = (mtx_region_t *)ascmalloc( sizeof(mtx_region_t) );
704 mtx->data->block[0].row.low = mtx->data->block[0].col.low = 0;
705 mtx->data->block[0].row.high = mtx->data->block[0].col.high = nrows - 1;
706 } else {
707 mtx->data->nblocks = 0;
708 mtx->data->block = NULL;
709 }
710
711 /**
712 *** Condense all unassigned rows on the outside of the block.
713 **/
714 ndxhigh = mtx->order-1;
715 for (nz.row = nrows; nz.row <= ndxhigh; ) {
716 if( NOTNULL(mtx->hdr.row[mtx->perm.row.cur_to_org[nz.row]]) ) ++nz.row;
717 else swap(&(mtx->perm.row),nz.row,ndxhigh--);
718 }
719 /**
720 *** Condense any incidenceless but wanted orgrows on
721 *** the outside of the assigned and unassigned block if sensible.
722 **/
723 if (nz.row < hirow) {
724 ndxhigh = mtx->order-1;
725 for (; nz.row <= ndxhigh; ) {
726 if( mtx->perm.row.cur_to_org[nz.row] < hirow ) ++nz.row;
727 else swap(&(mtx->perm.row),nz.row,ndxhigh--);
728 }
729 }
730 /**
731 *** Condense all unassigned columns on
732 *** the outside of the block.
733 **/
734 ndxhigh = mtx->order-1;
735 for (nz.col = nrows; nz.col <= ndxhigh; ) {
736 if( NOTNULL(mtx->hdr.col[mtx->perm.col.cur_to_org[nz.col]]) ) ++nz.col;
737 else swap(&(mtx->perm.col),nz.col,ndxhigh--);
738 }
739 /**
740 *** Condense any incidenceless but wanted orgcols on
741 *** the outside of the assigned and unassigned block if sensible.
742 **/
743 if (nz.col < hicol) {
744 ndxhigh = mtx->order-1;
745 for (; nz.col <= ndxhigh; ) {
746 if( mtx->perm.col.cur_to_org[nz.col] < hicol ) ++nz.col;
747 else swap(&(mtx->perm.col),nz.col,ndxhigh--);
748 }
749 }
750
751 if( NOTNULL(vars.row_visited) )
752 ascfree( vars.row_visited );
753 }
754
755 boolean mtx_output_assigned( mtx_matrix_t mtx)
756 {
757 if(!mtx_check_matrix(mtx)) return 0;
758 if (ISSLAVE(mtx)) {
759 mtx = mtx->master;
760 #if MTX_DEBUG
761 if( !mtx_check_matrix(mtx) ) return 0;
762 #endif
763 }
764 return (mtx->data->symbolic_rank > -1);
765 }
766
767 int32 mtx_symbolic_rank( mtx_matrix_t mtx)
768 {
769 #if MTX_DEBUG
770 if( !mtx_check_matrix(mtx) ) return -2;
771 #endif
772 if (ISSLAVE(mtx)) {
773 mtx = mtx->master;
774 #if MTX_DEBUG
775 if( !mtx_check_matrix(mtx) ) return -2;
776 #endif
777 }
778
779 if( !mtx_output_assigned(mtx) ) { /* matrix will be checked */
780 #if MTX_DEBUG
781 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_symbolic_rank\n");
782 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
783 #endif
784 }
785 return (mtx->data->symbolic_rank);
786 }
787
788 void mtx_set_symbolic_rank( mtx_matrix_t mtx, int32 rank)
789 {
790 mtx->data->symbolic_rank = rank;
791 }
792
793 boolean mtx_make_col_independent( mtx_matrix_t mtx, int32 col, mtx_range_t *rng)
794 {
795 struct assign_row_vars vars;
796 mtx_coord_t nz;
797 int32 ncol;
798
799 if (mtx!=NULL && ISSLAVE(mtx)) {
800 mtx = mtx->master;
801 #if MTX_DEBUG
802 if( !mtx_check_matrix(mtx) ) return -2;
803 #endif
804 }
805 if( !mtx_output_assigned(mtx) ) { /* matrix will be checked */
806 #if MTX_DEBUG
807 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_make_col_independent\n");
808 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
809 #endif
810 if(!mtx_check_matrix(mtx)) return -2;
811 }
812 if( col >= mtx->data->symbolic_rank ) return TRUE; /* col already ind */
813 if( rng->high < rng->low ) return FALSE; /* nobody to choose from */
814 if( rng->low < mtx->data->symbolic_rank ) return FALSE; /* bad choices */
815
816 vars.row_visited = mtx->order > ZERO ?
817 (int32 *)asccalloc( mtx->order,sizeof(int32) ) : NULL;
818 vars.rv_indicator = 1;
819 vars.unassigned_cols.low = rng->low;
820 vars.unassigned_cols.high = rng->high;
821 vars.assigned_cols.low = ZERO;
822 vars.assigned_cols.high = mtx->data->symbolic_rank-2;
823
824 /**
825 *** Move (col,col) to the bottom right corner
826 **/
827 nz.row = nz.col = col;
828 swap(&(mtx->perm.row),nz.row,mtx->data->symbolic_rank-1);
829 swap(&(mtx->perm.col),nz.col,mtx->data->symbolic_rank-1);
830
831 ncol= assign_row(mtx,mtx->data->symbolic_rank-1,&vars);
832 swap(&(mtx->perm.row),nz.row,mtx->data->symbolic_rank-1);
833 swap(&(mtx->perm.col),nz.col,mtx->data->symbolic_rank-1);
834 if(ncol < ZERO ) {
835 /* put it back the way it was */
836 return FALSE;
837 }
838
839 if( NOTNULL(mtx->data->block) && mtx->data->nblocks > 1 ) {
840 ascfree(mtx->data->block);
841 mtx->data->nblocks = 1;
842 mtx->data->block = (mtx_region_t *)ascmalloc( sizeof(mtx_region_t) );
843 mtx->data->block[0].row.low = ZERO;
844 mtx->data->block[0].col.low = ZERO;
845 mtx->data->block[0].row.high = mtx->data->symbolic_rank-1;
846 mtx->data->block[0].col.high = mtx->data->symbolic_rank-1;
847 }
848 if( NOTNULL(vars.row_visited) )
849 ascfree( vars.row_visited );
850 return TRUE;
851 }
852
853 struct analyze_row_vars {
854 int32 done;
855 int32 vstack;
856 int32 *lowlink;
857 int32 *blocksize;
858 int32 nblocks;
859 };
860
861 static
862 int32 analyze_row(mtx_matrix_t mtx, int32 row, struct analyze_row_vars *vars,
863 mtx_range_t *rng)
864 /*
865 * Places the row on the stack (row must be initially unvisited) and
866 * then determines lowlink[row]. If, in the end, lowlink[row] = row,
867 * then this row and all rows "above it on the stack" (i.e.
868 * vstack <= ndx <= row) belong to the same block, which is then moved
869 * from the stack to done (nblocks and blocksize updated). The new row
870 * number is returned (this may be different from the given row number
871 * since rows tend to move around).
872 */
873 {
874 mtx_range_t cols;
875 struct element_t Rewind, *elt;
876 int32 *tocur;
877
878 /* Push row onto the stack */
879 --(vars->vstack);
880 swap(&(mtx->perm.row),row,vars->vstack);
881 swap(&(mtx->perm.col),row,vars->vstack);
882 row = vars->vstack;
883 vars->lowlink[row] = row; /* Starting value of lowlink */
884
885 /*
886 * Scan the row for non-zeros. If unvisited indices are found, analyze
887 * them first (recursive call to this function) and in any case, if they
888 * are still on the stack, update lowlink[row] according to the new
889 * information. When the entire row is scanned, lowlink[row] should be
890 * completely computed.
891 */
892
893 /*
894 * THESE WERE: KAA_MODIFICATION
895 * cols.low = 0; and
896 * cols.high = mtx->data->symbolic_rank-1;
897 */
898 cols.low = rng->low;
899 cols.high = rng->high;
900 tocur = mtx->perm.col.org_to_cur;
901 Rewind.next.col = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
902 for( elt = &Rewind ; (elt=mtx_next_col(elt,&cols,tocur)) ; ) {
903 int32 col = tocur[elt->col];
904 if( vars->done <= col && col < vars->vstack ) {
905 /* Column unvisited, compute lowlink[col] */
906 #if SWAPS_PRESERVE_ORDER
907 col = analyze_row(mtx,col,vars,rng);
908 #else
909 if( col != analyze_row(mtx,col,vars,rng) ) {
910 elt = &Rewind;
911 continue;
912 }
913 #endif
914 }
915
916 if( col >= vars->vstack ) {
917 /* Still on stack: update lowlink[row] */
918 if( vars->lowlink[col] > vars->lowlink[row] )
919 vars->lowlink[row] = vars->lowlink[col];
920 } else {
921 /**
922 *** col is done (already in a block) and therefore
923 *** won't contribute anything to lowlink[row]
924 **/
925 }
926 }
927
928 if( vars->lowlink[row] == row ) {
929 /* Rows currently in vstack..row form the next block */
930 ++(vars->nblocks);
931 vars->blocksize[vars->done] = row - vars->vstack + 1; /* Block size */
932 /* Move rows/columns from the stack to done */
933 for( ; vars->vstack <= row ; ++(vars->vstack),++(vars->done) ) {
934 swap(&(mtx->perm.row),vars->vstack,vars->done);
935 swap(&(mtx->perm.col),vars->vstack,vars->done);
936 }
937 return( vars->done - 1 ); /* That's where the row ends up */
938 }
939 return(row);
940 }
941
942 static
943 int32 analyze_col(mtx_matrix_t mtx, int32 col, struct analyze_row_vars *vars,
944 mtx_range_t *rng)
945 /* (swap row for col here)
946 * Places the row on the stack (row must be initially unvisited) and
947 * then determines lowlink[row]. If, in the end, lowlink[row] = row,
948 * then this row and all rows "above it on the stack" (i.e.
949 * vstack <= ndx <= row) belong to the same block, which is then moved
950 * from the stack to done (nblocks and blocksize updated). The new row
951 * number is returned (this may be different from the given row number
952 * since rows tend to move around).
953 */
954 {
955 mtx_range_t rows;
956 struct element_t Rewind, *elt;
957 int32 *tocur;
958
959 /* Push row onto the stack */
960 --(vars->vstack);
961 swap(&(mtx->perm.col),col,vars->vstack);
962 swap(&(mtx->perm.row),col,vars->vstack);
963 col = vars->vstack;
964 vars->lowlink[col] = col; /* Starting value of lowlink */
965
966 /*
967 * Scan the row for non-zeros. If unvisited indices are found, analyze
968 * them first (recursive call to this function) and in any case, if they
969 * are still on the stack, update lowlink[row] according to the new
970 * information. When the entire row is scanned, lowlink[row] should be
971 * completely computed.
972 */
973
974 /* for full version:
975 * cols.low = 0; and
976 * cols.high = mtx->data->symbolic_rank-1;
977 */
978 rows.low = rng->low;
979 rows.high = rng->high;
980 tocur = mtx->perm.row.org_to_cur;
981 Rewind.next.row = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
982 for( elt = &Rewind ; (elt=mtx_next_row(elt,&rows,tocur)) ; ) {
983 int32 row = tocur[elt->row];
984 if( vars->done <= row && row < vars->vstack ) {
985 /* Column unvisited, compute lowlink[col] */
986 #if SWAPS_PRESERVE_ORDER
987 row = analyze_col(mtx,row,vars,rng);
988 #else
989 if( row != analyze_col(mtx,row,vars,rng) ) {
990 elt = &Rewind;
991 continue;
992 }
993 #endif
994 }
995
996 if( row >= vars->vstack ) {
997 /* Still on stack: update lowlink[row] */
998 if( vars->lowlink[row] > vars->lowlink[col] ) {
999 vars->lowlink[col] = vars->lowlink[row];
1000 }
1001 } else {
1002 /*
1003 * col is done (already in a block) and therefore
1004 * won't contribute anything to lowlink[row]
1005 */
1006 }
1007 }
1008
1009 if( vars->lowlink[col] == col ) {
1010 /* Rows currently in vstack..row form the next block */
1011 ++(vars->nblocks);
1012 vars->blocksize[vars->done] = col - vars->vstack + 1; /* Block size */
1013 /* Move rows/columns from the stack to done */
1014 for( ; vars->vstack <= col ; ++(vars->vstack),++(vars->done) ) {
1015 swap(&(mtx->perm.col),vars->vstack,vars->done);
1016 swap(&(mtx->perm.row),vars->vstack,vars->done);
1017 }
1018 return( vars->done - 1 ); /* That's where the row ends up */
1019 }
1020 return(col);
1021 }
1022
1023
1024 /* a function to check the diagonal for holes. if symbolic_rank is
1025 * set, and rng->high < rank, returns immediately.
1026 * Worst Cost: O(nnz) where nnz = incidence count sum over rows in rng.
1027 * Does not pass calls up to master.
1028 * Returns number of holes detected in diagonal in rng given.
1029 * If noisy != 0, writes the hole locations to g_mtxerr.
1030 * On detectably bad input, returns -1.
1031 */
1032 int32 mtx_full_diagonal(mtx_matrix_t mtx, mtx_range_t *rng, int noisy)
1033 {
1034 mtx_coord_t nz;
1035 mtx_range_t diag;
1036 int32 holecount=0;
1037
1038 if (mtx==NULL || rng == NULL) {
1039 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_full_diagonal sent NULL.\n");
1040 return -1;
1041 }
1042 if (mtx_output_assigned(mtx) && rng->high < mtx->data->symbolic_rank) {
1043 return 0;
1044 }
1045 if (!mtx_check_matrix(mtx)) {
1046 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_full_diagonal sent bad mtx.\n");
1047 return -1;
1048 }
1049 if (rng->low <0 || rng->high >= mtx->order) {
1050 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_full_diagonal sent bad rng.\n");
1051 return -1;
1052 }
1053 if (noisy) {
1054 for (nz.row=rng->low; nz.row <= rng->high; nz.row++) {
1055 diag.high = diag.low = nz.row;
1056 nz.col = mtx_FIRST;
1057 mtx_next_in_row(mtx,&nz,&diag); /* find diagonal, if there */
1058 if (nz.col==mtx_LAST) {
1059 holecount++;
1060 FPRINTF(g_mtxerr,"mtx: mtx_full_diagonal found hole at %d.\n",nz.row);
1061 }
1062 }
1063 } else {
1064 for (nz.row=rng->low; nz.row <= rng->high; nz.row++) {
1065 diag.high = diag.low = nz.row;
1066 nz.col = mtx_FIRST;
1067 mtx_next_in_row(mtx,&nz,&diag); /* find diagonal, if there */
1068 if (nz.col==mtx_LAST) holecount++;
1069 }
1070 }
1071 return holecount;
1072 }
1073
1074 /* so long as the mesh is sane, traversing it along one axis transposes
1075 * it in both directions. This must be called just before swapping headers
1076 * or the semantics are not right.
1077 */
1078 static
1079 void transpose_mesh(mtx_matrix_t mtx)
1080 {
1081 struct element_t *elttmp; /* swap var for nextelt pointers in elements */
1082 struct element_t *next, *elt;
1083 int32 ndxtmp; /* swap var for orgrow/orgcol indices */
1084 int32 nrows, ndx;
1085
1086 if (mtx == NULL) {
1087 return;
1088 }
1089 nrows = mtx->order;
1090 for( ndx = ZERO ; ndx < nrows ; ndx++) {
1091 next = mtx->hdr.row[ndx];
1092 /* traverse the row, turning it into a column behind us as we go. */
1093 while (next != NULL) {
1094 elt = next; /* the use of elt is redundant, but much more readable. */
1095 /* swap org row/col indices */
1096 ndxtmp = elt->col;
1097 elt->col = elt->row;
1098 elt->row = ndxtmp;
1099 /* swap next elt pointers */
1100 elttmp = elt->next.col;
1101 elt->next.col = elt->next.row;
1102 elt->next.row = elttmp;
1103 next = elttmp;
1104 }
1105 }
1106 }
1107
1108 int32 mtx_transpose(mtx_matrix_t mtx)
1109 {
1110 struct element_t **hdrtmp; /* swap var for headers */
1111 int32 *permtmp; /* swap var for permutation vectors. */
1112 int32 ndxtmp,bnum; /* swap var for indices */
1113 mtx_matrix_t master;
1114 int32 slave;
1115
1116
1117 #if MTX_DEBUG
1118 if( !mtx_check_matrix(mtx) ) {
1119 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_transpose\n");
1120 FPRINTF(g_mtxerr," Matrix given is in error.\n");
1121 return mtx_NONE; /*ben*/
1122 }
1123 #endif
1124 if (ISSLAVE(mtx)) {
1125 mtx = mtx->master;
1126 #if MTX_DEBUG
1127 if( !mtx_check_matrix(mtx) ) {
1128 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_transpose\n");
1129 FPRINTF(g_mtxerr," Matrix given is in error.\n");
1130 return mtx_NONE;
1131 }
1132 #endif
1133 }
1134 master = mtx;
1135 if (master->perm.transpose == mtx_NONE) {
1136 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_transpose\n");
1137 FPRINTF(g_mtxerr," Matrix given is 0 order.\n");
1138 return mtx_NONE;
1139 }
1140 master->perm.transpose = !(master->perm.transpose);
1141 /* swap perms on master */
1142 permtmp = mtx->perm.col.org_to_cur; /* do o2v */
1143 mtx->perm.col.org_to_cur = mtx->perm.row.org_to_cur;
1144 mtx->perm.row.org_to_cur = permtmp;
1145 permtmp = mtx->perm.col.cur_to_org; /* do v2o */
1146 mtx->perm.col.cur_to_org = mtx->perm.row.cur_to_org;
1147 mtx->perm.row.cur_to_org = permtmp;
1148
1149 /* transpose mesh on master */
1150 transpose_mesh(mtx);
1151
1152 /* swap headers on master */
1153 hdrtmp = mtx->hdr.col;
1154 mtx->hdr.col = mtx->hdr.row;
1155 mtx->hdr.row = hdrtmp;
1156
1157 /* transpose blocks on master */
1158 if (mtx->data->block != NULL) {
1159 for (bnum = 0; bnum < mtx->data->nblocks; bnum++) {
1160 /* upper left */
1161 ndxtmp = mtx->data->block[bnum].row.low;
1162 mtx->data->block[bnum].row.low = mtx->data->block[bnum].col.low;
1163 mtx->data->block[bnum].col.low = ndxtmp;
1164 /* lower right */
1165 ndxtmp = mtx->data->block[bnum].row.high;
1166 mtx->data->block[bnum].row.high = mtx->data->block[bnum].col.high;
1167 mtx->data->block[bnum].col.high = ndxtmp;
1168 }
1169 }
1170
1171 /* swap hdrs, update perm references, and do mesh on slaves */
1172 for (slave = 0; slave < master->nslaves; slave++) {
1173 mtx = master->slaves[slave];
1174 mtx->perm.row.org_to_cur = master->perm.row.org_to_cur;
1175 mtx->perm.row.cur_to_org = master->perm.row.cur_to_org;
1176 mtx->perm.col.org_to_cur = master->perm.col.org_to_cur;
1177 mtx->perm.col.cur_to_org = master->perm.col.cur_to_org;
1178 transpose_mesh(mtx);
1179 hdrtmp = mtx->hdr.col;
1180 mtx->hdr.col = mtx->hdr.row;
1181 mtx->hdr.row = hdrtmp;
1182 }
1183
1184 return 0;
1185 }
1186
1187 int32 mtx_isa_transpose(mtx_matrix_t mtx)
1188 {
1189 #if MTX_DEBUG
1190 if( !mtx_check_matrix(mtx) ) {
1191 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_isa_transpose\n");
1192 FPRINTF(g_mtxerr," Matrix given is in error.\n");
1193 return mtx_NONE; /*ben*/
1194 }
1195 #endif
1196 if (ISSLAVE(mtx)) {
1197 mtx = mtx->master;
1198 #if MTX_DEBUG
1199 if( !mtx_check_matrix(mtx) ) {
1200 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_isa_transpose\n");
1201 FPRINTF(g_mtxerr," Matrix given's master is in error.\n");
1202 return mtx_NONE;
1203 }
1204 #endif
1205 }
1206 return mtx->perm.transpose;
1207 }
1208
1209 /*
1210 * This function uses the same algorithm as its sister function
1211 * mtx_partition (at the moment). It is more general in that it
1212 * explicitly takes a region, and will return a block_list
1213 * structure rather than manipulating the one associated with
1214 * the matrix. It now works (we hope) on an arbitrary square
1215 * region.
1216 */
1217 mtx_block_t *mtx_block_partition( mtx_matrix_t mtx, mtx_region_t *reg)
1218 {
1219 struct analyze_row_vars vars;
1220 int32 blocknum, start, size;
1221 mtx_block_t *blocklist;
1222
1223 if (mtx!=NULL && ISSLAVE(mtx)) {
1224 mtx = mtx->master;
1225 }
1226
1227 if (!mtx_output_assigned(mtx)) {
1228 int mc = mtx_check_matrix(mtx);
1229 if (!mc) return NULL;
1230 if (mtx_full_diagonal(mtx,&(reg->row),1)) {
1231 FPRINTF(g_mtxerr,"WARNING: (mtx) mtx_block_partition\n");
1232 FPRINTF(g_mtxerr," Assignment bad. partitioning may be bad..\n");
1233 }
1234 }
1235
1236 vars.vstack = reg->row.high + 1; /* was symbolic_rank */
1237 if (vars.vstack > 0) {
1238 vars.lowlink = vars.blocksize =
1239 (int32 *)ascmalloc(vars.vstack*sizeof(int32));
1240 } else {
1241 vars.lowlink = vars.blocksize = NULL;
1242 }
1243 vars.done = 0;
1244 vars.nblocks = 0;
1245
1246 /*
1247 * The core of the partition routine. Analyze row
1248 * now takes a range, which represents a square region.
1249 */
1250 while( vars.done < vars.vstack ) {
1251 analyze_row(mtx,vars.done,&vars,&(reg->row));
1252 }
1253
1254 /*
1255 * Prepare the block list structure and copy the data
1256 * from into it from vars.
1257 */
1258 if (vars.nblocks > 0) {
1259 blocklist = (mtx_block_t *)ascmalloc(sizeof(mtx_block_t));
1260 blocklist->nblocks = vars.nblocks;
1261 blocklist->block = (mtx_region_t *)
1262 ascmalloc(vars.nblocks * sizeof(mtx_region_t));
1263 } else {
1264 return NULL;
1265 }
1266
1267 for(start = blocknum = 0; blocknum < vars.nblocks; blocknum++) {
1268 size = vars.blocksize[start];
1269 blocklist->block[blocknum].row.low =
1270 blocklist->block[blocknum].col.low = start;
1271 blocklist->block[blocknum].row.high =
1272 blocklist->block[blocknum].col.high = start + size - 1;
1273 start += size;
1274 }
1275 if(NOTNULL(vars.lowlink))
1276 ascfree((POINTER)vars.lowlink);
1277
1278 return blocklist;
1279 }
1280
1281
1282 void mtx_partition( mtx_matrix_t mtx)
1283 /**
1284 *** Matrix must be previously output assigned. Then the assigned
1285 *** region (symbolic basis) portion of the matrix is partitioned
1286 *** leaving unassigned rows and columns on the outside.
1287 ***
1288 *** vstack -----------------------+
1289 *** done -------+ |
1290 *** V V
1291 *** done | unvisited | stack (visited)
1292 ***
1293 *** 0 <= "done" < done <= "unvisited" < vstack <= "stack" < order
1294 ***
1295 *** For all ndx on the stack, lowlink[ndx] "points to" the deepest part of
1296 *** the stack (highest index) which can be reached from row ndx by the
1297 *** associated directed graph (with edges r-->c, where (r,c) is non-zero).
1298 ***
1299 *** For blocks which are already completed, blocksize[start] = block size,
1300 *** where start is the starting row/column index of the block. In every
1301 *** case, 0 <= start < done, so lowlink and blocksize can be the same
1302 *** array.
1303 **/
1304 {
1305 struct analyze_row_vars vars;
1306 int32 blocknum, start, size;
1307 mtx_range_t rng;
1308
1309 if (mtx!=NULL && ISSLAVE(mtx)) {
1310 mtx = mtx->master;
1311 }
1312 if( !mtx_output_assigned(mtx) ) {
1313 int mc=mtx_check_matrix(mtx);
1314 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_partition\n");
1315 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
1316 if(!mc) return;
1317 }
1318
1319 vars.vstack = mtx->data->symbolic_rank;
1320 vars.lowlink = vars.blocksize = vars.vstack > 0 ?
1321 (int32 *)ascmalloc( vars.vstack*sizeof(int32) ) : NULL;
1322 vars.done = 0;
1323 vars.nblocks = 0;
1324
1325 rng.low = 0; /* KAA_MODIFICATION */
1326 rng.high = mtx->data->symbolic_rank-1;
1327 while( vars.done < vars.vstack )
1328 analyze_row(mtx,vars.done,&vars,&rng);
1329
1330 if( NOTNULL(mtx->data->block) )
1331 ascfree(mtx->data->block);
1332 mtx->data->nblocks = vars.nblocks;
1333 mtx->data->block = vars.nblocks > 0 ?
1334 (mtx_region_t *)ascmalloc( vars.nblocks*sizeof(mtx_region_t) ) : NULL;
1335 for( start=blocknum=0 ; blocknum < vars.nblocks ; blocknum++ ) {
1336 size = vars.blocksize[start];
1337 mtx->data->block[blocknum].row.low =
1338 mtx->data->block[blocknum].col.low = start;
1339 mtx->data->block[blocknum].row.high =
1340 mtx->data->block[blocknum].col.high = start + size - 1;
1341 start += size;
1342 }
1343 if( NOTNULL(vars.lowlink) )
1344 ascfree( (POINTER)vars.lowlink );
1345 }
1346
1347 void mtx_ut_partition( mtx_matrix_t mtx)
1348 /*
1349 * This is the upper triangular version of the blt mtx_partition function.
1350 * Swap 'row' and 'column' or r and c in all related comments.
1351 * Matrix must be previously output assigned. Then the assigned
1352 * region (symbolic basis) portion of the matrix is partitioned
1353 * leaving unassigned rows and columns on the outside.
1354 *
1355 * vstack -----------------------+
1356 * done -------+ |
1357 * V V
1358 * done | unvisited | stack (visited)
1359 *
1360 * 0 <= "done" < done <= "unvisited" < vstack <= "stack" < order
1361 *
1362 * For all ndx on the stack, lowlink[ndx] "points to" the deepest part of
1363 * the stack (highest index) which can be reached from row ndx by the
1364 * associated directed graph (with edges r-->c, where (r,c) is non-zero).
1365 *
1366 * For blocks which are already completed, blocksize[start] = block size,
1367 * where start is the starting row/column index of the block. In every
1368 * case, 0 <= start < done, so lowlink and blocksize can be the same
1369 * array.
1370 */
1371 {
1372 struct analyze_row_vars vars;
1373 int32 blocknum, start, size;
1374 mtx_range_t rng;
1375
1376 if (mtx!=NULL && ISSLAVE(mtx)) {
1377 mtx = mtx->master;
1378 }
1379 if( !mtx_output_assigned(mtx) ) {
1380 int mc=mtx_check_matrix(mtx);
1381 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_ut_partition\n");
1382 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
1383 if(!mc) return;
1384 }
1385
1386 vars.vstack = mtx->data->symbolic_rank;
1387 vars.lowlink = vars.blocksize = vars.vstack > 0 ?
1388 (int32 *)ascmalloc( vars.vstack*sizeof(int32) ) : NULL;
1389 vars.done = 0;
1390 vars.nblocks = 0;
1391
1392 rng.low = 0;
1393 rng.high = mtx->data->symbolic_rank-1;
1394 while( vars.done < vars.vstack ) {
1395 analyze_col(mtx,vars.done,&vars,&rng);
1396 }
1397
1398 if( NOTNULL(mtx->data->block) ) {
1399 ascfree(mtx->data->block);
1400 }
1401 mtx->data->nblocks = vars.nblocks;
1402 mtx->data->block = vars.nblocks > 0 ?
1403 (mtx_region_t *)ascmalloc( vars.nblocks*sizeof(mtx_region_t) ) : NULL;
1404 for( start=blocknum=0 ; blocknum < vars.nblocks ; blocknum++ ) {
1405 size = vars.blocksize[start];
1406 mtx->data->block[blocknum].row.low =
1407 mtx->data->block[blocknum].col.low = start;
1408 mtx->data->block[blocknum].row.high =
1409 mtx->data->block[blocknum].col.high = start + size - 1;
1410 start += size;
1411 }
1412 if( NOTNULL(vars.lowlink) ) {
1413 ascfree( (POINTER)vars.lowlink );
1414 }
1415 }
1416
1417 #ifdef DELETE_THIS_DEAD_FUNCTION
1418 /* a comparator meeting the qsort spec exactly for any pair i1,i2 */
1419 static int mtx_cmpint(int *i1, int *i2)
1420 {
1421 if (*i1 < *i2) return -1;
1422 return (*i1 > *i2);
1423 }
1424 #endif /* DELETE_THIS_DEAD_FUNCTION */
1425
1426 /* a faster comparator which assumes i1 =/= i2 ever */
1427 static int mtx_cmpintfast(int *i1, int *i2)
1428 {
1429 if (*i1 < *i2) return -1;
1430 return 1;
1431 }
1432
1433 /* this could be faster, but as it's a 1 time deal, who cares?
1434 we're at the mercy of the system qsort. */
1435 /* ah, don't you hate prototypes? now it looks like this
1436 function matters... */
1437 void mtx_org_permute(mtx_matrix_t mtx, mtx_region_t *reg)
1438 {
1439 if (ISNULL(mtx)) {
1440 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_org_permute\n");
1441 FPRINTF(g_mtxerr," NULL matrix given.\n");
1442 return;
1443 }
1444 if (ISSLAVE(mtx)) {
1445 mtx = mtx->master;
1446 }
1447 if (ISNULL(reg)) {
1448 mtx_reset_perm(mtx);
1449 } else {
1450 int32 i,j,len;
1451 int32 *sort,*top;
1452 len=reg->col.high-reg->col.low+1;
1453 if (len>1) {
1454 /* create sort space */
1455 top=sort=(int *)ascmalloc(sizeof(int)*len);
1456 if (ISNULL(sort)) {
1457 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_org_permute\n");
1458 FPRINTF(g_mtxerr," Insufficient memory. Not permuted.\n");
1459 return;
1460 }
1461 /* copy current org ordering into array */
1462 for (i = reg->col.low; i <= reg->col.high; i++) {
1463 *top = mtx->perm.col.cur_to_org[i];
1464 top++;
1465 }
1466 /* sort org columns */
1467 qsort((void *)sort,(size_t)len,sizeof(int),
1468 (int (*)(const void *,const void *))mtx_cmpintfast);
1469 /* permute columns */
1470 top = sort;
1471 for (i= reg->col.low; i < reg->col.high; i++ ) {
1472 if (mtx->perm.col.cur_to_org[i] != *top) {
1473 j = mtx->perm.col.org_to_cur[*top];
1474 swap(&(mtx->perm.col),i,j);
1475 }
1476 top++;
1477 }
1478 ascfree(sort);
1479 }
1480 /* ditto for rows */
1481 len=reg->row.high-reg->row.low+1;
1482 if (len>1) {
1483 top=sort=(int *)ascmalloc(sizeof(int)*len);
1484 if (ISNULL(sort)) {
1485 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_org_permute\n");
1486 FPRINTF(g_mtxerr," Insufficient memory. Not permuted.\n");
1487 return;
1488 }
1489 for (i = reg->row.low; i <= reg->row.high; i++) {
1490 *top = mtx->perm.row.cur_to_org[i];
1491 top++;
1492 }
1493 qsort((void *)sort,(size_t)len,sizeof(int),
1494 (int (*)(const void *,const void *))mtx_cmpintfast);
1495 top = sort;
1496 for (i= reg->row.low; i < reg->row.high; i++ ) {
1497 if (mtx->perm.row.cur_to_org[i] != *top) {
1498 j = mtx->perm.row.org_to_cur[*top];
1499 swap(&(mtx->perm.row),i,j);
1500 }
1501 top++;
1502 }
1503 ascfree(sort);
1504 }
1505 }
1506 }
1507
1508 /***********************************************************************\
1509 end structural analysis stuff
1510 \***********************************************************************/
1511 boolean mtx_check_blocks( mtx_matrix_t mtx)
1512 {
1513 mtx_range_t invalid;
1514 int32 blocknum;
1515
1516 if (ISSLAVE(mtx)) {
1517 mtx = mtx->master;
1518 }
1519 if( !mtx_output_assigned(mtx) ) { /* matrix will be checked */
1520 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_check_blocks\n");
1521 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
1522 return( FALSE );
1523 }
1524
1525 invalid.high = mtx->data->symbolic_rank-1;
1526 /* Last block need not be checked */
1527 for( blocknum = 0 ; blocknum < mtx->data->nblocks-1 ; ++blocknum ) {
1528 mtx_coord_t nz;
1529 invalid.low = mtx->data->block[blocknum].row.high + 1;
1530 for( nz.row=mtx->data->block[blocknum].row.low ;
1531 nz.row<=mtx->data->block[blocknum].row.high ;
1532 ++nz.row ) {
1533 struct element_t Rewind;
1534 Rewind.next.col = mtx->hdr.row[mtx->perm.row.cur_to_org[nz.row]];
1535 if( mtx_next_col(&Rewind,&invalid,mtx->perm.col.org_to_cur) ) {
1536 mtx->data->nblocks = -1;
1537 if( NOTNULL(mtx->data->block) )
1538 ascfree(mtx->data->block);
1539 return(FALSE);
1540 }
1541 }
1542 }
1543 return(TRUE);
1544 }
1545
1546 int32 mtx_number_of_blocks( mtx_matrix_t mtx)
1547 {
1548 if( !mtx_output_assigned(mtx) ) { /* matrix will be checked */
1549 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_number_of_blocks\n");
1550 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
1551 }
1552
1553 #if MTX_DEBUG
1554 if( !mtx_check_blocks(mtx) ) {
1555 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_number_of_blocks\n");
1556 FPRINTF(g_mtxerr," Invalid partition.\n");
1557 return mtx_NONE;
1558 }
1559 #endif
1560
1561 return (mtx->data->nblocks);
1562 }
1563
1564 int32 mtx_block( mtx_matrix_t mtx, int32 block_number, mtx_region_t *block)
1565 {
1566 mtx_region_t *reg;
1567 int error=0;
1568
1569 if( !mtx_output_assigned(mtx) ) { /* matrix will be checked */
1570 #if MTX_DEBUG
1571 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_block\n");
1572 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
1573 #endif
1574 error=1;
1575 }
1576
1577 #if MTX_DEBUG
1578 if( !mtx_check_blocks(mtx) ) {
1579 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_block\n");
1580 FPRINTF(g_mtxerr," Invalid partition.\n");
1581 }
1582 #endif
1583
1584 if( (block_number > mtx->data->nblocks-1) || (block_number<0) ) {
1585 #if MTX_DEBUG
1586 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_block\n");
1587 FPRINTF(g_mtxerr," Block number doesn't exist.\n");
1588 #endif
1589 error=1;
1590 }
1591 if (error) {
1592 block->row.low=0;
1593 block->col.low=0;
1594 block->row.high=mtx->order-1;
1595 block->col.high=mtx->order-1;
1596 return mtx_NONE;
1597 }
1598
1599 reg = &(mtx->data->block[block_number]);
1600 mem_copy_cast(reg,block,sizeof(mtx_region_t));
1601 return 0;
1602 }
1603
1604 int32 mtx_block_containing_row(mtx_matrix_t mtx, int32 row, mtx_region_t *block)
1605 {
1606 int32 blow,bhigh;
1607 mtx_region_t *reg;
1608
1609 if( !mtx_output_assigned(mtx) ) { /* matrix will be checked */
1610 #if MTX_DEBUG
1611 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_block_containing_row\n");
1612 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
1613 #endif
1614 return 0;
1615 }
1616
1617 #if MTX_DEBUG
1618 if( !mtx_check_blocks(mtx) ) {
1619 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_block_containing_row\n");
1620 FPRINTF(g_mtxerr," Invalid partition.\n");
1621 }
1622 #endif
1623
1624 blow = 0;
1625 bhigh = mtx->data->nblocks-1;
1626 while( blow <= bhigh ) {
1627 int32 block_number = (blow+bhigh)/2;
1628 if( row > mtx->data->block[block_number].row.high ) {
1629 blow = block_number+1;
1630 } else if( row < mtx->data->block[block_number].row.low ) {
1631 bhigh = block_number-1;
1632 } else {
1633 reg = &(mtx->data->block[block_number]);
1634 mem_copy_cast(reg,block,sizeof(mtx_region_t));
1635 return(block_number);
1636 }
1637 }
1638 return(mtx_NONE);
1639 }
1640
1641 int32 mtx_block_containing_col(mtx_matrix_t mtx, int32 col, mtx_region_t *block)
1642 {
1643 int32 blow,bhigh;
1644 mtx_region_t *reg;
1645
1646 if( !mtx_output_assigned(mtx) ) { /* matrix will be checked */
1647 #if MTX_DEBUG
1648 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_block_containing_col\n");
1649 FPRINTF(g_mtxerr," Matrix not output assigned.\n");
1650 #endif
1651 return 0;
1652 }
1653
1654 #if MTX_DEBUG
1655 if( !mtx_check_blocks(mtx) ) {
1656 FPRINTF(g_mtxerr,"ERROR: (mtx) mtx_block_containing_col\n");
1657 FPRINTF(g_mtxerr," Invalid partition.\n");
1658 }
1659 #endif
1660
1661 blow = 0;
1662 bhigh = mtx->data->nblocks-1;
1663 while( blow <= bhigh ) {
1664 int32 block_number = (blow+bhigh)/2;
1665 if( col > mtx->data->block[block_number].col.high ) {
1666 blow = block_number+1;
1667 } else if( col < mtx->data->block[block_number].col.low ) {
1668 bhigh = block_number-1;
1669 } else {
1670 reg = &(mtx->data->block[block_number]);
1671 mem_copy_cast(reg,block,sizeof(mtx_region_t));
1672 return(block_number);
1673 }
1674 }
1675 return(mtx_NONE);
1676 }
1677
1678 #undef __MTX_C_SEEN__

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