1 2 #include <../src/mat/impls/aij/seq/aij.h> 3 4 /* 5 MatToSymmetricIJ_SeqAIJ - Convert a (generally nonsymmetric) sparse AIJ matrix 6 to IJ format (ignore the "A" part) Allocates the space needed. Uses only 7 the lower triangular part of the matrix. 8 9 Description: 10 Take the data in the row-oriented sparse storage and build the 11 IJ data for the Matrix. Return 0 on success,row + 1 on failure 12 at that row. Produces the ij for a symmetric matrix by using 13 the lower triangular part of the matrix if lower_triangular is PETSC_TRUE; 14 it uses the upper triangular otherwise. 15 16 Input Parameters: 17 . Matrix - matrix to convert 18 . lower_triangular - symmetrize the lower triangular part 19 . shiftin - the shift for the original matrix (0 or 1) 20 . shiftout - the shift required for the ordering routine (0 or 1) 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 may be freed with PetscFree(); 28 This routine is provided for ordering routines that require a 29 symmetric structure. It is required since those routines call 30 SparsePak routines that expect a symmetric matrix. 31 */ 32 PetscErrorCode MatToSymmetricIJ_SeqAIJ(PetscInt m,PetscInt *ai,PetscInt *aj,PetscBool lower_triangular,PetscInt shiftin,PetscInt shiftout,PetscInt **iia,PetscInt **jja) 33 { 34 PetscInt *work,*ia,*ja,*j,i,nz,row,col; 35 36 PetscFunctionBegin; 37 /* allocate space for row pointers */ 38 PetscCall(PetscCalloc1(m+1,&ia)); 39 *iia = ia; 40 PetscCall(PetscMalloc1(m+1,&work)); 41 42 /* determine the number of columns in each row */ 43 ia[0] = shiftout; 44 for (row = 0; row < m; row++) { 45 nz = ai[row+1] - ai[row]; 46 j = aj + ai[row] + shiftin; 47 while (nz--) { 48 col = *j++ + shiftin; 49 if (lower_triangular) { 50 if (col > row) break; 51 } else { 52 if (col < row) break; 53 } 54 if (col != row) ia[row+1]++; 55 ia[col+1]++; 56 } 57 } 58 59 /* shiftin ia[i] to point to next row */ 60 for (i=1; i<m+1; i++) { 61 row = ia[i-1]; 62 ia[i] += row; 63 work[i-1] = row - shiftout; 64 } 65 66 /* allocate space for column pointers */ 67 nz = ia[m] + (!shiftin); 68 PetscCall(PetscMalloc1(nz,&ja)); 69 *jja = ja; 70 71 /* loop over lower triangular part putting into ja */ 72 for (row = 0; row < m; row++) { 73 nz = ai[row+1] - ai[row]; 74 j = aj + ai[row] + shiftin; 75 while (nz--) { 76 col = *j++ + shiftin; 77 if (lower_triangular) { 78 if (col > row) break; 79 } else { 80 if (col < row) break; 81 } 82 if (col != row) ja[work[col]++] = row + shiftout; 83 ja[work[row]++] = col + shiftout; 84 } 85 } 86 PetscCall(PetscFree(work)); 87 PetscFunctionReturn(0); 88 } 89