Actual source code: sparse.cxx
1: #include <petscmat.h>
2: #include <adolc/adolc_sparse.h>
4: /*
5: REQUIRES configuration of PETSc with option --download-adolc.
7: For documentation on ADOL-C, see
8: $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
9: */
11: /*
12: Basic printing for sparsity pattern
14: Input parameters:
15: comm - MPI communicator
16: m - number of rows
18: Output parameter:
19: sparsity - matrix sparsity pattern, typically computed using an ADOL-C function such as jac_pat
20: */
21: PetscErrorCode PrintSparsity(MPI_Comm comm, PetscInt m, unsigned int **sparsity)
22: {
23: PetscFunctionBegin;
24: PetscCall(PetscPrintf(comm, "Sparsity pattern:\n"));
25: for (PetscInt i = 0; i < m; i++) {
26: PetscCall(PetscPrintf(comm, "\n %2d: ", i));
27: for (PetscInt j = 1; j <= (PetscInt)sparsity[i][0]; j++) PetscCall(PetscPrintf(comm, " %2d ", sparsity[i][j]));
28: }
29: PetscCall(PetscPrintf(comm, "\n\n"));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: /*
34: Generate a seed matrix defining the partition of columns of a matrix by a particular coloring,
35: used for matrix compression
37: Input parameter:
38: iscoloring - the index set coloring to be used
40: Output parameter:
41: S - the resulting seed matrix
43: Notes:
44: Before calling GenerateSeedMatrix, Seed should be allocated as a logically 2d array with number of
45: rows equal to the matrix to be compressed and number of columns equal to the number of colors used
46: in iscoloring.
47: */
48: PetscErrorCode GenerateSeedMatrix(ISColoring iscoloring, PetscScalar **S)
49: {
50: IS *is;
51: PetscInt p, size;
52: const PetscInt *indices;
54: PetscFunctionBegin;
55: PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &p, &is));
56: for (PetscInt colour = 0; colour < p; colour++) {
57: PetscCall(ISGetLocalSize(is[colour], &size));
58: PetscCall(ISGetIndices(is[colour], &indices));
59: for (PetscInt j = 0; j < size; j++) S[indices[j]][colour] = 1.;
60: PetscCall(ISRestoreIndices(is[colour], &indices));
61: }
62: PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: /*
67: Establish a look-up vector whose entries contain the colour used for that diagonal entry. Clearly
68: we require the number of dependent and independent variables to be equal in this case.
69: Input parameters:
70: S - the seed matrix defining the coloring
71: sparsity - the sparsity pattern of the matrix to be recovered, typically computed using an ADOL-C
72: function, such as jac_pat or hess_pat
73: m - the number of rows of Seed (and the matrix to be recovered)
74: p - the number of colors used (also the number of columns in Seed)
75: Output parameter:
76: R - the recovery vector to be used for de-compression
77: */
78: PetscErrorCode GenerateSeedMatrixPlusRecovery(ISColoring iscoloring, PetscScalar **S, PetscScalar *R)
79: {
80: IS *is;
81: PetscInt p, size, colour, j;
82: const PetscInt *indices;
84: PetscFunctionBegin;
85: PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &p, &is));
86: for (colour = 0; colour < p; colour++) {
87: PetscCall(ISGetLocalSize(is[colour], &size));
88: PetscCall(ISGetIndices(is[colour], &indices));
89: for (j = 0; j < size; j++) {
90: S[indices[j]][colour] = 1.;
91: R[indices[j]] = colour;
92: }
93: PetscCall(ISRestoreIndices(is[colour], &indices));
94: }
95: PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: /*
100: Establish a look-up matrix whose entries contain the column coordinates of the corresponding entry
101: in a matrix which has been compressed using the coloring defined by some seed matrix
103: Input parameters:
104: S - the seed matrix defining the coloring
105: sparsity - the sparsity pattern of the matrix to be recovered, typically computed using an ADOL-C
106: function, such as jac_pat or hess_pat
107: m - the number of rows of Seed (and the matrix to be recovered)
108: p - the number of colors used (also the number of columns in Seed)
110: Output parameter:
111: R - the recovery matrix to be used for de-compression
112: */
113: PetscErrorCode GetRecoveryMatrix(PetscScalar **S, unsigned int **sparsity, PetscInt m, PetscInt p, PetscScalar **R)
114: {
115: PetscInt i, j, k, colour;
117: PetscFunctionBegin;
118: for (i = 0; i < m; i++) {
119: for (colour = 0; colour < p; colour++) {
120: R[i][colour] = -1.;
121: for (k = 1; k <= (PetscInt)sparsity[i][0]; k++) {
122: j = (PetscInt)sparsity[i][k];
123: if (S[j][colour] == 1.) {
124: R[i][colour] = j;
125: break;
126: }
127: }
128: }
129: }
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /*
134: Recover the values of a sparse matrix from a compressed format and insert these into a matrix
136: Input parameters:
137: mode - use INSERT_VALUES or ADD_VALUES, as required
138: m - number of rows of matrix.
139: p - number of colors used in the compression of J (also the number of columns of R)
140: R - recovery matrix to use in the decompression procedure
141: C - compressed matrix to recover values from
142: a - shift value for implicit problems (select NULL or unity for explicit problems)
144: Output parameter:
145: A - Mat to be populated with values from compressed matrix
146: */
147: PetscErrorCode RecoverJacobian(Mat A, InsertMode mode, PetscInt m, PetscInt p, PetscScalar **R, PetscScalar **C, PetscReal *a)
148: {
149: PetscFunctionBegin;
150: for (PetscInt i = 0; i < m; i++) {
151: for (PetscInt colour = 0; colour < p; colour++) {
152: PetscInt j = (PetscInt)R[i][colour];
153: if (j != -1) {
154: if (a) C[i][colour] *= *a;
155: PetscCall(MatSetValues(A, 1, &i, 1, &j, &C[i][colour], mode));
156: }
157: }
158: }
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: /*
163: Recover the values of the local portion of a sparse matrix from a compressed format and insert
164: these into the local portion of a matrix
166: Input parameters:
167: mode - use INSERT_VALUES or ADD_VALUES, as required
168: m - number of rows of matrix.
169: p - number of colors used in the compression of J (also the number of columns of R)
170: R - recovery matrix to use in the decompression procedure
171: C - compressed matrix to recover values from
172: a - shift value for implicit problems (select NULL or unity for explicit problems)
174: Output parameter:
175: A - Mat to be populated with values from compressed matrix
176: */
177: PetscErrorCode RecoverJacobianLocal(Mat A, InsertMode mode, PetscInt m, PetscInt p, PetscScalar **R, PetscScalar **C, PetscReal *a)
178: {
179: PetscFunctionBegin;
180: for (PetscInt i = 0; i < m; i++) {
181: for (PetscInt colour = 0; colour < p; colour++) {
182: PetscInt j = (PetscInt)R[i][colour];
183: if (j != -1) {
184: if (a) C[i][colour] *= *a;
185: PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &C[i][colour], mode));
186: }
187: }
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: /*
193: Recover the diagonal of the Jacobian from its compressed matrix format
195: Input parameters:
196: mode - use INSERT_VALUES or ADD_VALUES, as required
197: m - number of rows of matrix.
198: R - recovery vector to use in the decompression procedure
199: C - compressed matrix to recover values from
200: a - shift value for implicit problems (select NULL or unity for explicit problems)
202: Output parameter:
203: diag - Vec to be populated with values from compressed matrix
204: */
205: PetscErrorCode RecoverDiagonal(Vec diag, InsertMode mode, PetscInt m, PetscScalar *R, PetscScalar **C, PetscReal *a)
206: {
207: PetscFunctionBegin;
208: for (PetscInt i = 0; i < m; i++) {
209: PetscInt colour = (PetscInt)R[i];
210: if (a) C[i][colour] *= *a;
211: PetscCall(VecSetValues(diag, 1, &i, &C[i][colour], mode));
212: }
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
216: /*
217: Recover the local portion of the diagonal of the Jacobian from its compressed matrix format
219: Input parameters:
220: mode - use INSERT_VALUES or ADD_VALUES, as required
221: m - number of rows of matrix.
222: R - recovery vector to use in the decompression procedure
223: C - compressed matrix to recover values from
224: a - shift value for implicit problems (select NULL or unity for explicit problems)
226: Output parameter:
227: diag - Vec to be populated with values from compressed matrix
228: */
229: PetscErrorCode RecoverDiagonalLocal(Vec diag, InsertMode mode, PetscInt m, PetscScalar *R, PetscScalar **C, PetscReal *a)
230: {
231: PetscFunctionBegin;
232: for (PetscInt i = 0; i < m; i++) {
233: PetscInt colour = (PetscInt)R[i];
234: if (a) C[i][colour] *= *a;
235: PetscCall(VecSetValuesLocal(diag, 1, &i, &C[i][colour], mode));
236: }
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }