1 #include <petscdmda.h> 2 #include <petscts.h> 3 #include <adolc/interfaces.h> 4 #include "contexts.cxx" 5 6 /* 7 REQUIRES configuration of PETSc with option --download-adolc. 8 9 For documentation on ADOL-C, see 10 $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf 11 */ 12 13 /* 14 ADOL-C implementation for Jacobian vector product, using the forward mode of AD. 15 Intended to overload MatMult in matrix-free methods where implicit timestepping 16 has been used. 17 18 For an implicit Jacobian we may use the rule that 19 G = M*xdot + f(x) ==> dG/dx = a*M + df/dx, 20 where a = d(xdot)/dx is a constant. Evaluated at x0 and acting upon a vector x1: 21 (dG/dx)(x0) * x1 = (a*df/d(xdot)(x0) + df/dx)(x0) * x1. 22 23 Input parameters: 24 A_shell - Jacobian matrix of MatShell type 25 X - vector to be multiplied by A_shell 26 27 Output parameters: 28 Y - product of A_shell and X 29 */ 30 PetscErrorCode PetscAdolcIJacobianVectorProduct(Mat A_shell,Vec X,Vec Y) 31 { 32 AdolcMatCtx *mctx; 33 PetscErrorCode ierr; 34 PetscInt m,n,i,j,k = 0,d; 35 const PetscScalar *x0; 36 PetscScalar *action,*x1; 37 Vec localX1; 38 DM da; 39 DMDALocalInfo info; 40 41 PetscFunctionBegin; 42 43 /* Get matrix-free context info */ 44 ierr = MatShellGetContext(A_shell,&mctx);CHKERRQ(ierr); 45 m = mctx->m; 46 n = mctx->n; 47 48 /* Get local input vectors and extract data, x0 and x1*/ 49 ierr = TSGetDM(mctx->ts,&da);CHKERRQ(ierr); 50 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 51 ierr = DMGetLocalVector(da,&localX1);CHKERRQ(ierr); 52 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX1);CHKERRQ(ierr); 53 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX1);CHKERRQ(ierr); 54 55 ierr = VecGetArrayRead(mctx->localX0,&x0);CHKERRQ(ierr); 56 ierr = VecGetArray(localX1,&x1);CHKERRQ(ierr); 57 58 /* dF/dx part */ 59 ierr = PetscMalloc1(m,&action);CHKERRQ(ierr); 60 ierr = PetscLogEventBegin(mctx->event1,0,0,0,0);CHKERRQ(ierr); 61 fos_forward(mctx->tag1,m,n,0,x0,x1,NULL,action); 62 for (j=info.gys; j<info.gys+info.gym; j++) { 63 for (i=info.gxs; i<info.gxs+info.gxm; i++) { 64 for (d=0; d<2; d++) { 65 if ((i >= info.xs) && (i < info.xs+info.xm) && (j >= info.ys) && (j < info.ys+info.ym)) { 66 ierr = VecSetValuesLocal(Y,1,&k,&action[k],INSERT_VALUES);CHKERRQ(ierr); 67 } 68 k++; 69 } 70 } 71 } 72 ierr = PetscLogEventEnd(mctx->event1,0,0,0,0);CHKERRQ(ierr); 73 k = 0; 74 ierr = VecAssemblyBegin(Y);CHKERRQ(ierr); /* Note: Need to assemble between separate calls */ 75 ierr = VecAssemblyEnd(Y);CHKERRQ(ierr); /* to INSERT_VALUES and ADD_VALUES */ 76 77 /* a * dF/d(xdot) part */ 78 ierr = PetscLogEventBegin(mctx->event2,0,0,0,0);CHKERRQ(ierr); 79 fos_forward(mctx->tag2,m,n,0,x0,x1,NULL,action); 80 for (j=info.gys; j<info.gys+info.gym; j++) { 81 for (i=info.gxs; i<info.gxs+info.gxm; i++) { 82 for (d=0; d<2; d++) { 83 if ((i >= info.xs) && (i < info.xs+info.xm) && (j >= info.ys) && (j < info.ys+info.ym)) { 84 action[k] *= mctx->shift; 85 ierr = VecSetValuesLocal(Y,1,&k,&action[k],ADD_VALUES);CHKERRQ(ierr); 86 } 87 k++; 88 } 89 } 90 } 91 ierr = PetscLogEventEnd(mctx->event2,0,0,0,0);CHKERRQ(ierr); 92 ierr = VecAssemblyBegin(Y);CHKERRQ(ierr); 93 ierr = VecAssemblyEnd(Y);CHKERRQ(ierr); 94 ierr = PetscFree(action);CHKERRQ(ierr); 95 96 /* Restore local vector */ 97 ierr = VecRestoreArray(localX1,&x1);CHKERRQ(ierr); 98 ierr = VecRestoreArrayRead(mctx->localX0,&x0);CHKERRQ(ierr); 99 ierr = DMRestoreLocalVector(da,&localX1);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 /* 104 ADOL-C implementation for Jacobian vector product, using the forward mode of AD. 105 Intended to overload MatMult in matrix-free methods where implicit timestepping 106 has been applied to a problem of the form 107 du/dt = F(u). 108 109 Input parameters: 110 A_shell - Jacobian matrix of MatShell type 111 X - vector to be multiplied by A_shell 112 113 Output parameters: 114 Y - product of A_shell and X 115 */ 116 PetscErrorCode PetscAdolcIJacobianVectorProductIDMass(Mat A_shell,Vec X,Vec Y) 117 { 118 AdolcMatCtx *mctx; 119 PetscErrorCode ierr; 120 PetscInt m,n,i,j,k = 0,d; 121 const PetscScalar *x0; 122 PetscScalar *action,*x1; 123 Vec localX1; 124 DM da; 125 DMDALocalInfo info; 126 127 PetscFunctionBegin; 128 129 /* Get matrix-free context info */ 130 ierr = MatShellGetContext(A_shell,&mctx);CHKERRQ(ierr); 131 m = mctx->m; 132 n = mctx->n; 133 134 /* Get local input vectors and extract data, x0 and x1*/ 135 ierr = TSGetDM(mctx->ts,&da);CHKERRQ(ierr); 136 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 137 ierr = DMGetLocalVector(da,&localX1);CHKERRQ(ierr); 138 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX1);CHKERRQ(ierr); 139 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX1);CHKERRQ(ierr); 140 141 ierr = VecGetArrayRead(mctx->localX0,&x0);CHKERRQ(ierr); 142 ierr = VecGetArray(localX1,&x1);CHKERRQ(ierr); 143 144 /* dF/dx part */ 145 ierr = PetscMalloc1(m,&action);CHKERRQ(ierr); 146 ierr = PetscLogEventBegin(mctx->event1,0,0,0,0);CHKERRQ(ierr); 147 fos_forward(mctx->tag1,m,n,0,x0,x1,NULL,action); 148 for (j=info.gys; j<info.gys+info.gym; j++) { 149 for (i=info.gxs; i<info.gxs+info.gxm; i++) { 150 for (d=0; d<2; d++) { 151 if ((i >= info.xs) && (i < info.xs+info.xm) && (j >= info.ys) && (j < info.ys+info.ym)) { 152 ierr = VecSetValuesLocal(Y,1,&k,&action[k],INSERT_VALUES);CHKERRQ(ierr); 153 } 154 k++; 155 } 156 } 157 } 158 ierr = PetscLogEventEnd(mctx->event1,0,0,0,0);CHKERRQ(ierr); 159 k = 0; 160 ierr = VecAssemblyBegin(Y);CHKERRQ(ierr); /* Note: Need to assemble between separate calls */ 161 ierr = VecAssemblyEnd(Y);CHKERRQ(ierr); /* to INSERT_VALUES and ADD_VALUES */ 162 ierr = PetscFree(action);CHKERRQ(ierr); 163 164 /* Restore local vector */ 165 ierr = VecRestoreArray(localX1,&x1);CHKERRQ(ierr); 166 ierr = VecRestoreArrayRead(mctx->localX0,&x0);CHKERRQ(ierr); 167 ierr = DMRestoreLocalVector(da,&localX1);CHKERRQ(ierr); 168 169 /* a * dF/d(xdot) part */ 170 ierr = PetscLogEventBegin(mctx->event2,0,0,0,0);CHKERRQ(ierr); 171 ierr = VecAXPY(Y,mctx->shift,X);CHKERRQ(ierr); 172 ierr = PetscLogEventEnd(mctx->event2,0,0,0,0);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 /* 177 ADOL-C implementation for Jacobian transpose vector product, using the reverse mode of AD. 178 Intended to overload MatMultTranspose in matrix-free methods where implicit timestepping 179 has been used. 180 181 Input parameters: 182 A_shell - Jacobian matrix of MatShell type 183 Y - vector to be multiplied by A_shell transpose 184 185 Output parameters: 186 X - product of A_shell transpose and X 187 */ 188 PetscErrorCode PetscAdolcIJacobianTransposeVectorProduct(Mat A_shell,Vec Y,Vec X) 189 { 190 AdolcMatCtx *mctx; 191 PetscErrorCode ierr; 192 PetscInt m,n,i,j,k = 0,d; 193 const PetscScalar *x; 194 PetscScalar *action,*y; 195 Vec localY; 196 DM da; 197 DMDALocalInfo info; 198 199 PetscFunctionBegin; 200 201 /* Get matrix-free context info */ 202 ierr = MatShellGetContext(A_shell,&mctx);CHKERRQ(ierr); 203 m = mctx->m; 204 n = mctx->n; 205 206 /* Get local input vectors and extract data, x0 and x1*/ 207 ierr = TSGetDM(mctx->ts,&da);CHKERRQ(ierr); 208 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 209 ierr = DMGetLocalVector(da,&localY);CHKERRQ(ierr); 210 ierr = DMGlobalToLocalBegin(da,Y,INSERT_VALUES,localY);CHKERRQ(ierr); 211 ierr = DMGlobalToLocalEnd(da,Y,INSERT_VALUES,localY);CHKERRQ(ierr); 212 ierr = VecGetArrayRead(mctx->localX0,&x);CHKERRQ(ierr); 213 ierr = VecGetArray(localY,&y);CHKERRQ(ierr); 214 215 /* dF/dx part */ 216 ierr = PetscMalloc1(n,&action);CHKERRQ(ierr); 217 ierr = PetscLogEventBegin(mctx->event3,0,0,0,0);CHKERRQ(ierr); 218 if (!mctx->flg) 219 zos_forward(mctx->tag1,m,n,1,x,NULL); 220 fos_reverse(mctx->tag1,m,n,y,action); 221 for (j=info.gys; j<info.gys+info.gym; j++) { 222 for (i=info.gxs; i<info.gxs+info.gxm; i++) { 223 for (d=0; d<2; d++) { 224 if ((i >= info.xs) && (i < info.xs+info.xm) && (j >= info.ys) && (j < info.ys+info.ym)) { 225 ierr = VecSetValuesLocal(X,1,&k,&action[k],INSERT_VALUES);CHKERRQ(ierr); 226 } 227 k++; 228 } 229 } 230 } 231 ierr = PetscLogEventEnd(mctx->event3,0,0,0,0);CHKERRQ(ierr); 232 k = 0; 233 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); /* Note: Need to assemble between separate calls */ 234 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); /* to INSERT_VALUES and ADD_VALUES */ 235 236 /* a * dF/d(xdot) part */ 237 ierr = PetscLogEventBegin(mctx->event4,0,0,0,0);CHKERRQ(ierr); 238 if (!mctx->flg) { 239 zos_forward(mctx->tag2,m,n,1,x,NULL); 240 mctx->flg = PETSC_TRUE; 241 } 242 fos_reverse(mctx->tag2,m,n,y,action); 243 for (j=info.gys; j<info.gys+info.gym; j++) { 244 for (i=info.gxs; i<info.gxs+info.gxm; i++) { 245 for (d=0; d<2; d++) { 246 if ((i >= info.xs) && (i < info.xs+info.xm) && (j >= info.ys) && (j < info.ys+info.ym)) { 247 action[k] *= mctx->shift; 248 ierr = VecSetValuesLocal(X,1,&k,&action[k],ADD_VALUES);CHKERRQ(ierr); 249 } 250 k++; 251 } 252 } 253 } 254 ierr = PetscLogEventEnd(mctx->event4,0,0,0,0);CHKERRQ(ierr); 255 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); 256 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); 257 ierr = PetscFree(action);CHKERRQ(ierr); 258 259 /* Restore local vector */ 260 ierr = VecRestoreArray(localY,&y);CHKERRQ(ierr); 261 ierr = VecRestoreArrayRead(mctx->localX0,&x);CHKERRQ(ierr); 262 ierr = DMRestoreLocalVector(da,&localY);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 /* 267 ADOL-C implementation for Jacobian transpose vector product, using the reverse mode of AD. 268 Intended to overload MatMultTranspose in matrix-free methods where implicit timestepping 269 has been applied to a problem of the form 270 du/dt = F(u). 271 272 Input parameters: 273 A_shell - Jacobian matrix of MatShell type 274 Y - vector to be multiplied by A_shell transpose 275 276 Output parameters: 277 X - product of A_shell transpose and X 278 */ 279 PetscErrorCode PetscAdolcIJacobianTransposeVectorProductIDMass(Mat A_shell,Vec Y,Vec X) 280 { 281 AdolcMatCtx *mctx; 282 PetscErrorCode ierr; 283 PetscInt m,n,i,j,k = 0,d; 284 const PetscScalar *x; 285 PetscScalar *action,*y; 286 Vec localY; 287 DM da; 288 DMDALocalInfo info; 289 290 PetscFunctionBegin; 291 292 /* Get matrix-free context info */ 293 ierr = MatShellGetContext(A_shell,&mctx);CHKERRQ(ierr); 294 m = mctx->m; 295 n = mctx->n; 296 297 /* Get local input vectors and extract data, x0 and x1*/ 298 ierr = TSGetDM(mctx->ts,&da);CHKERRQ(ierr); 299 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 300 ierr = DMGetLocalVector(da,&localY);CHKERRQ(ierr); 301 ierr = DMGlobalToLocalBegin(da,Y,INSERT_VALUES,localY);CHKERRQ(ierr); 302 ierr = DMGlobalToLocalEnd(da,Y,INSERT_VALUES,localY);CHKERRQ(ierr); 303 ierr = VecGetArrayRead(mctx->localX0,&x);CHKERRQ(ierr); 304 ierr = VecGetArray(localY,&y);CHKERRQ(ierr); 305 306 /* dF/dx part */ 307 ierr = PetscMalloc1(n,&action);CHKERRQ(ierr); 308 ierr = PetscLogEventBegin(mctx->event3,0,0,0,0);CHKERRQ(ierr); 309 if (!mctx->flg) zos_forward(mctx->tag1,m,n,1,x,NULL); 310 fos_reverse(mctx->tag1,m,n,y,action); 311 for (j=info.gys; j<info.gys+info.gym; j++) { 312 for (i=info.gxs; i<info.gxs+info.gxm; i++) { 313 for (d=0; d<2; d++) { 314 if ((i >= info.xs) && (i < info.xs+info.xm) && (j >= info.ys) && (j < info.ys+info.ym)) { 315 ierr = VecSetValuesLocal(X,1,&k,&action[k],INSERT_VALUES);CHKERRQ(ierr); 316 } 317 k++; 318 } 319 } 320 } 321 ierr = PetscLogEventEnd(mctx->event3,0,0,0,0);CHKERRQ(ierr); 322 k = 0; 323 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); /* Note: Need to assemble between separate calls */ 324 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); /* to INSERT_VALUES and ADD_VALUES */ 325 ierr = PetscFree(action);CHKERRQ(ierr); 326 327 /* Restore local vector */ 328 ierr = VecRestoreArray(localY,&y);CHKERRQ(ierr); 329 ierr = VecRestoreArrayRead(mctx->localX0,&x);CHKERRQ(ierr); 330 ierr = DMRestoreLocalVector(da,&localY);CHKERRQ(ierr); 331 332 /* a * dF/d(xdot) part */ 333 ierr = PetscLogEventBegin(mctx->event4,0,0,0,0);CHKERRQ(ierr); 334 ierr = VecAXPY(X,mctx->shift,Y);CHKERRQ(ierr); 335 ierr = PetscLogEventEnd(mctx->event4,0,0,0,0);CHKERRQ(ierr); 336 PetscFunctionReturn(0); 337 } 338