1 #include <petscmat.h>
2 #include <adolc/adolc_sparse.h>
3
4 /*
5 REQUIRES configuration of PETSc with option --download-adolc.
6
7 For documentation on ADOL-C, see
8 $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
9 */
10
11 /*
12 Basic printing for sparsity pattern
13
14 Input parameters:
15 comm - MPI communicator
16 m - number of rows
17
18 Output parameter:
19 sparsity - matrix sparsity pattern, typically computed using an ADOL-C function such as jac_pat
20 */
PrintSparsity(MPI_Comm comm,PetscInt m,unsigned int ** sparsity)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 }
32
33 /*
34 Generate a seed matrix defining the partition of columns of a matrix by a particular coloring,
35 used for matrix compression
36
37 Input parameter:
38 iscoloring - the index set coloring to be used
39
40 Output parameter:
41 S - the resulting seed matrix
42
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 */
GenerateSeedMatrix(ISColoring iscoloring,PetscScalar ** S)48 PetscErrorCode GenerateSeedMatrix(ISColoring iscoloring, PetscScalar **S)
49 {
50 IS *is;
51 PetscInt p, size;
52 const PetscInt *indices;
53
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 }
65
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 */
GenerateSeedMatrixPlusRecovery(ISColoring iscoloring,PetscScalar ** S,PetscScalar * R)78 PetscErrorCode GenerateSeedMatrixPlusRecovery(ISColoring iscoloring, PetscScalar **S, PetscScalar *R)
79 {
80 IS *is;
81 PetscInt p, size, colour, j;
82 const PetscInt *indices;
83
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 }
98
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
102
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)
109
110 Output parameter:
111 R - the recovery matrix to be used for de-compression
112 */
GetRecoveryMatrix(PetscScalar ** S,unsigned int ** sparsity,PetscInt m,PetscInt p,PetscScalar ** R)113 PetscErrorCode GetRecoveryMatrix(PetscScalar **S, unsigned int **sparsity, PetscInt m, PetscInt p, PetscScalar **R)
114 {
115 PetscInt i, j, k, colour;
116
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 }
132
133 /*
134 Recover the values of a sparse matrix from a compressed format and insert these into a matrix
135
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)
143
144 Output parameter:
145 A - Mat to be populated with values from compressed matrix
146 */
RecoverJacobian(Mat A,InsertMode mode,PetscInt m,PetscInt p,PetscScalar ** R,PetscScalar ** C,PetscReal * a)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 }
161
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
165
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)
173
174 Output parameter:
175 A - Mat to be populated with values from compressed matrix
176 */
RecoverJacobianLocal(Mat A,InsertMode mode,PetscInt m,PetscInt p,PetscScalar ** R,PetscScalar ** C,PetscReal * a)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 }
191
192 /*
193 Recover the diagonal of the Jacobian from its compressed matrix format
194
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)
201
202 Output parameter:
203 diag - Vec to be populated with values from compressed matrix
204 */
RecoverDiagonal(Vec diag,InsertMode mode,PetscInt m,PetscScalar * R,PetscScalar ** C,PetscReal * a)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 }
215
216 /*
217 Recover the local portion of the diagonal of the Jacobian from its compressed matrix format
218
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)
225
226 Output parameter:
227 diag - Vec to be populated with values from compressed matrix
228 */
RecoverDiagonalLocal(Vec diag,InsertMode mode,PetscInt m,PetscScalar * R,PetscScalar ** C,PetscReal * a)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 }
239