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