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