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

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

Parent Directory Parent Directory | Revision Log Revision Log


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