Actual source code: ij.c
2: #include <../src/mat/impls/aij/seq/aij.h>
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.
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.
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)
22: Output Parameters:
23: . ia - ia part of IJ representation (row information)
24: . ja - ja part (column indices)
26: Note:
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;
36: PetscFunctionBegin;
37: /* allocate space for row pointers */
38: PetscCall(PetscCalloc1(m + 1, &ia));
39: *iia = ia;
40: PetscCall(PetscMalloc1(m + 1, &work));
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: }
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: }
66: /* allocate space for column pointers */
67: nz = ia[m] + (!shiftin);
68: PetscCall(PetscMalloc1(nz, &ja));
69: *jja = ja;
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(PETSC_SUCCESS);
88: }