1 #ifndef lint 2 static char vcid[] = "$Id: ij.c,v 1.7 1995/06/08 03:09:10 bsmith Exp bsmith $"; 3 #endif 4 5 #include "aij.h" 6 7 /* 8 MatToSymmetricIJ_AIJ - Convert a sparse AIJ matrix to IJ format 9 (ignore the "A" part) 10 Allocates the space needed. Uses only the lower triangular 11 part of the matrix. 12 13 Description: 14 Take the data in the row-oriented sparse storage and build the 15 IJ data for the Matrix. Return 0 on success, row + 1 on failure 16 at that row. Produces the ij for a symmetric matrix by only using 17 the lower triangular part of the matrix. 18 19 Input Parameters: 20 . Matrix - matrix to convert 21 22 Output Parameters: 23 . ia - ia part of IJ representation (row information) 24 . ja - ja part (column indices) 25 26 Notes: 27 $ Both ia and ja maybe freed with PETSCFREE(); 28 $ This routine is provided for ordering routines that require a 29 $ symmetric structure. It is used in SpOrder (and derivatives) since 30 $ those routines call SparsePak routines that expect a symmetric 31 $ matrix. 32 */ 33 int MatToSymmetricIJ_AIJ( Mat_AIJ *Matrix, int **iia, int **jja ) 34 { 35 int *work,*ia,*ja,*j,i, nz, n, row, wr; 36 register int col; 37 38 n = Matrix->m; 39 40 /* allocate space for row pointers */ 41 *iia = ia = (int *) PETSCMALLOC( (n+1)*sizeof(int) ); CHKPTRQ(ia); 42 PETSCMEMSET(ia,0,(n+1)*sizeof(int)); 43 work = (int *) PETSCMALLOC( (n+1)*sizeof(int) ); CHKPTRQ(work); 44 45 /* determine the number of columns in each row */ 46 ia[0] = 1; 47 for (row = 0; row < n; row++) { 48 nz = Matrix->i[row+1] - Matrix->i[row]; 49 j = Matrix->j + Matrix->i[row] - 1; 50 while (nz--) { 51 col = *j++ - 1; 52 if ( col > row ) { 53 break; 54 } 55 if (col != row) ia[row+1]++; 56 ia[col+1]++; 57 } 58 } 59 60 /* shift ia[i] to point to next row */ 61 for ( i=1; i<n+1; i++ ) { 62 row = ia[i-1]; 63 ia[i] += row; 64 work[i-1] = row - 1; 65 } 66 67 /* allocate space for column pointers */ 68 nz = ia[n]; 69 *jja = ja = (int *) PETSCMALLOC( nz*sizeof(int) ); CHKPTRQ(ja); 70 71 /* loop over lower triangular part putting into ja */ 72 for (row = 0; row < n; row++) { 73 nz = Matrix->i[row+1] - Matrix->i[row]; 74 j = Matrix->j + Matrix->i[row] - 1; 75 while (nz--) { 76 col = *j++ - 1; 77 if ( col > row ) { 78 break; 79 } 80 if (col != row) {wr = work[col]; work[col] = wr + 1; 81 ja[wr] = row + 1; } 82 wr = work[row]; work[row] = wr + 1; 83 ja[wr] = col + 1; 84 } 85 86 } 87 PETSCFREE(work); 88 return 0; 89 } 90 91