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 */
MatToSymmetricIJ_SeqAIJ(PetscInt m,PetscInt * ai,PetscInt * aj,PetscBool lower_triangular,PetscInt shiftin,PetscInt shiftout,PetscInt ** iia,PetscInt ** jja)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