1 #ifndef lint 2 static char vcid[] = "$Id: ij.c,v 1.10 1995/09/21 20:10:16 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) 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 $ Both ia and ja maybe freed with PETSCFREE(); 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 MatToSymmetricIJ_SeqAIJ( Mat_SeqAIJ *A, int **iia, int **jja ) 34 { 35 int *work,*ia,*ja,*j,i, nz, n = A->m, row, wr, col, shift = A->indexshift; 36 37 /* allocate space for row pointers */ 38 *iia = ia = (int *) PETSCMALLOC( (n+1)*sizeof(int) ); CHKPTRQ(ia); 39 PetscZero(ia,(n+1)*sizeof(int)); 40 work = (int *) PETSCMALLOC( (n+1)*sizeof(int) ); CHKPTRQ(work); 41 42 /* determine the number of columns in each row */ 43 ia[0] = 1; 44 for (row = 0; row < n; row++) { 45 nz = A->i[row+1] - A->i[row]; 46 j = A->j + A->i[row] + shift; 47 while (nz--) { 48 col = *j++ + shift; 49 if ( col > row ) { break;} 50 if (col != row) ia[row+1]++; 51 ia[col+1]++; 52 } 53 } 54 55 /* shift ia[i] to point to next row */ 56 for ( i=1; i<n+1; i++ ) { 57 row = ia[i-1]; 58 ia[i] += row; 59 work[i-1] = row - 1; 60 } 61 62 /* allocate space for column pointers */ 63 nz = ia[n] + (!shift); 64 *jja = ja = (int *) PETSCMALLOC( nz*sizeof(int) ); CHKPTRQ(ja); 65 66 /* loop over lower triangular part putting into ja */ 67 for (row = 0; row < n; row++) { 68 nz = A->i[row+1] - A->i[row]; 69 j = A->j + A->i[row] + shift; 70 while (nz--) { 71 col = *j++ + shift; 72 if ( col > row ) { break;} 73 if (col != row) {wr = work[col]; work[col] = wr + 1; ja[wr] = row + 1; } 74 wr = work[row]; work[row] = wr + 1; 75 ja[wr] = col + 1; 76 } 77 } 78 PETSCFREE(work); 79 return 0; 80 } 81 82