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: /* allocate space for row pointers */
37: PetscCalloc1(m + 1, &ia);
38: *iia = ia;
39: 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: 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: PetscFree(work);
86: return 0;
87: }