1 #ifndef lint 2 static char vcid[] = "$Id: ij.c,v 1.1 1994/03/18 00:27:00 gropp Exp $"; 3 #endif 4 5 #include "aij.h" 6 7 /* 8 SpToSymmetricIJ - 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 This routine is provided for ordering routines that require a 28 symmetric structure. It is used in SpOrder (and derivatives) since 29 those routines call SparsePak routines that expect a symmetric 30 matrix. 31 */ 32 int SpToSymmetricIJ( Matiaij *Matrix, int **iia, int **jja ) 33 { 34 int *work,*ia,*ja,*j,i, nz, n, row, wr; 35 register int col; 36 37 n = Matrix->m; 38 39 /* allocate space for row pointers */ 40 *iia = ia = (int *) MALLOC( (n+1)*sizeof(int) ); CHKPTR(ia); 41 MEMSET(ia,0,(n+1)*sizeof(int)); 42 work = (int *) MALLOC( (n+1)*sizeof(int) ); CHKPTR(work); 43 44 /* determine the number of columns in each row */ 45 ia[0] = 1; 46 for (row = 0; row < n; row++) { 47 nz = Matrix->i[row+1] - Matrix->i[row]; 48 j = Matrix->j + Matrix->i[row] - 1; 49 while (nz--) { 50 col = *j++; 51 if ( col > row ) { 52 break; 53 } 54 if (col != row) ia[row+1]++; 55 ia[col+1]++; 56 } 57 } 58 59 /* shift ia[i] to point to next row */ 60 for ( i=1; i<n+1; i++ ) { 61 row = ia[i-1]; 62 ia[i] += row; 63 work[i-1] = row - 1; 64 } 65 66 /* allocate space for column pointers */ 67 nz = ia[n]; 68 *jja = ja = (int *) MALLOC( nz*sizeof(int) ); CHKPTR(ja); 69 70 /* loop over lower triangular part putting into ja */ 71 for (row = 0; row < n; row++) { 72 nz = Matrix->i[row+1] - Matrix->i[row]; 73 j = Matrix->j + Matrix->i[row] - 1; 74 while (nz--) { 75 col = *j++; 76 if ( col > row ) { 77 break; 78 } 79 if (col != row) {wr = work[col]; work[col] = wr + 1; 80 ja[wr] = row + 1; } 81 wr = work[row]; work[row] = wr + 1; 82 ja[wr] = col + 1; 83 } 84 } 85 FREE(work); 86 return 0; 87 } 88 89