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