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