1 #ifndef lint 2 static char vcid[] = "$Id: ij.c,v 1.13 1995/11/01 23:17:58 bsmith Exp bsmith $"; 3 #endif 4 5 #include "aij.h" 6 7 /* 8 MatToSymmetricIJ_SeqAIJ - Convert a sparse AIJ matrix to IJ format 9 (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 21 Output Parameters: 22 . ia - ia part of IJ representation (row information) 23 . ja - ja part (column indices) 24 25 Notes: 26 $ Both ia and ja may be freed with PetscFree(); 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 MatToSymmetricIJ_SeqAIJ( Mat_SeqAIJ *A, int **iia, int **jja ) 33 { 34 int *work,*ia,*ja,*j,i, nz, n = A->m, row, wr, col, shift = A->indexshift; 35 36 /* allocate space for row pointers */ 37 *iia = ia = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(ia); 38 PetscMemzero(ia,(n+1)*sizeof(int)); 39 work = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(work); 40 41 /* determine the number of columns in each row */ 42 ia[0] = 1; 43 for (row = 0; row < n; row++) { 44 nz = A->i[row+1] - A->i[row]; 45 j = A->j + A->i[row] + shift; 46 while (nz--) { 47 col = *j++ + shift; 48 if ( col > row ) { break;} 49 if (col != row) ia[row+1]++; 50 ia[col+1]++; 51 } 52 } 53 54 /* shift ia[i] to point to next row */ 55 for ( i=1; i<n+1; i++ ) { 56 row = ia[i-1]; 57 ia[i] += row; 58 work[i-1] = row - 1; 59 } 60 61 /* allocate space for column pointers */ 62 nz = ia[n] + (!shift); 63 *jja = ja = (int *) PetscMalloc( nz*sizeof(int) ); CHKPTRQ(ja); 64 65 /* loop over lower triangular part putting into ja */ 66 for (row = 0; row < n; row++) { 67 nz = A->i[row+1] - A->i[row]; 68 j = A->j + A->i[row] + shift; 69 while (nz--) { 70 col = *j++ + shift; 71 if ( col > row ) { break;} 72 if (col != row) {wr = work[col]; work[col] = wr + 1; ja[wr] = row + 1; } 73 wr = work[row]; work[row] = wr + 1; 74 ja[wr] = col + 1; 75 } 76 } 77 PetscFree(work); 78 return 0; 79 } 80 81 82 83