1 #ifndef lint 2 static char vcid[] = "$Id: $"; 3 #endif 4 5 6 #include "aij.h" 7 8 /* 9 SpToSymmetricIJ - Convert a sparse AIJ matrix to IJ format 10 (ignore the "A" part) 11 Allocates the space needed. Uses only the lower triangular 12 part of the matrix. 13 14 Description: 15 Take the data in the row-oriented sparse storage and build the 16 IJ data for the Matrix. Return 0 on success, row + 1 on failure 17 at that row. Produces the ij for a symmetric matrix by only using 18 the lower triangular part of the matrix. 19 20 Input Parameters: 21 . Matrix - matrix to convert 22 23 Output Parameters: 24 . ia - ia part of IJ representation (row information) 25 . ja - ja part (column indices) 26 27 Notes: 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 SpToSymmetricIJ( Matiaij *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 *) MALLOC( (n+1)*sizeof(int) ); CHKPTR(ia); 42 MEMSET(ia,0,(n+1)*sizeof(int)); 43 work = (int *) MALLOC( (n+1)*sizeof(int) ); CHKPTR(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 *) MALLOC( nz*sizeof(int) ); CHKPTR(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 FREE(work); 88 return 0; 89 } 90 91