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