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