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