1 /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */ 2 3 /* 4 Include "petsctao.h" so we can use TAO solvers. 5 petscdmda.h for distributed array 6 */ 7 #include <petsctao.h> 8 #include <petscdmda.h> 9 10 static char help[] = "This example demonstrates use of the TAO package to \n\ 11 solve an unconstrained minimization problem. This example is based on a \n\ 12 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\ 13 boundary values along the edges of the domain, the objective is to find the\n\ 14 surface with the minimal area that satisfies the boundary conditions.\n\ 15 The command line options are:\n\ 16 -da_grid_x <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ 17 -da_grid_y <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ 18 -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\ 19 for an average of the boundary conditions\n\n"; 20 21 /* 22 User-defined application context - contains data needed by the 23 application-provided call-back routines, FormFunction(), 24 FormFunctionGradient(), and FormHessian(). 25 */ 26 typedef struct { 27 PetscInt mx, my; /* discretization in x, y directions */ 28 PetscReal *bottom, *top, *left, *right; /* boundary values */ 29 DM dm; /* distributed array data structure */ 30 Mat H; /* Hessian */ 31 } AppCtx; 32 33 /* -------- User-defined Routines --------- */ 34 35 static PetscErrorCode MSA_BoundaryConditions(AppCtx *); 36 static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec); 37 static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat); 38 static PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *); 39 static PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *); 40 static PetscErrorCode FormGradient(Tao, Vec, Vec, void *); 41 static PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *); 42 static PetscErrorCode My_Monitor(Tao, void *); 43 44 int main(int argc, char **argv) 45 { 46 Vec x; /* solution, gradient vectors */ 47 PetscBool viewmat; /* flags */ 48 PetscBool fddefault, fdcoloring; /* flags */ 49 Tao tao; /* TAO solver context */ 50 AppCtx user; /* user-defined work context */ 51 ISColoring iscoloring; 52 MatFDColoring matfdcoloring; 53 54 /* Initialize TAO */ 55 PetscFunctionBeginUser; 56 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 57 58 /* Create distributed array (DM) to manage parallel grid and vectors */ 59 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &user.dm)); 60 PetscCall(DMSetFromOptions(user.dm)); 61 PetscCall(DMSetUp(user.dm)); 62 PetscCall(DMDASetUniformCoordinates(user.dm, -0.5, 0.5, -0.5, 0.5, PETSC_DECIDE, PETSC_DECIDE)); 63 PetscCall(DMDAGetInfo(user.dm, PETSC_IGNORE, &user.mx, &user.my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 64 65 PetscCall(PetscPrintf(MPI_COMM_WORLD, "\n---- Minimum Surface Area Problem -----\n")); 66 PetscCall(PetscPrintf(MPI_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my)); 67 68 /* Create TAO solver and set desired solution method.*/ 69 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 70 PetscCall(TaoSetType(tao, TAOCG)); 71 72 /* 73 Extract global vector from DA for the vector of variables -- PETSC routine 74 Compute the initial solution -- application specific, see below 75 Set this vector for use by TAO -- TAO routine 76 */ 77 PetscCall(DMCreateGlobalVector(user.dm, &x)); 78 PetscCall(MSA_BoundaryConditions(&user)); 79 PetscCall(MSA_InitialPoint(&user, x)); 80 PetscCall(TaoSetSolution(tao, x)); 81 82 /* 83 Initialize the Application context for use in function evaluations -- application specific, see below. 84 Set routines for function and gradient evaluation 85 */ 86 PetscCall(TaoSetObjective(tao, FormFunction, (void *)&user)); 87 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user)); 88 89 /* 90 Given the command line arguments, calculate the hessian with either the user- 91 provided function FormHessian, or the default finite-difference driven Hessian 92 functions 93 */ 94 PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fddefault", &fddefault)); 95 PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fdcoloring", &fdcoloring)); 96 97 /* 98 Create a matrix data structure to store the Hessian and set 99 the Hessian evaluation routine. 100 Set the matrix nonzero structure to be used for Hessian evaluations 101 */ 102 PetscCall(DMCreateMatrix(user.dm, &user.H)); 103 PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE)); 104 105 if (fdcoloring) { 106 PetscCall(DMCreateColoring(user.dm, IS_COLORING_GLOBAL, &iscoloring)); 107 PetscCall(MatFDColoringCreate(user.H, iscoloring, &matfdcoloring)); 108 PetscCall(ISColoringDestroy(&iscoloring)); 109 PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode (*)(void))FormGradient, (void *)&user)); 110 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 111 PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessianColor, (void *)matfdcoloring)); 112 } else if (fddefault) { 113 PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessian, (void *)NULL)); 114 } else { 115 PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user)); 116 } 117 118 /* 119 If my_monitor option is in command line, then use the user-provided 120 monitoring function 121 */ 122 PetscCall(PetscOptionsHasName(NULL, NULL, "-my_monitor", &viewmat)); 123 if (viewmat) PetscCall(TaoMonitorSet(tao, My_Monitor, NULL, NULL)); 124 125 /* Check for any tao command line options */ 126 PetscCall(TaoSetFromOptions(tao)); 127 128 /* SOLVE THE APPLICATION */ 129 PetscCall(TaoSolve(tao)); 130 131 PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD)); 132 133 /* Free TAO data structures */ 134 PetscCall(TaoDestroy(&tao)); 135 136 /* Free PETSc data structures */ 137 PetscCall(VecDestroy(&x)); 138 PetscCall(MatDestroy(&user.H)); 139 if (fdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring)); 140 PetscCall(PetscFree(user.bottom)); 141 PetscCall(PetscFree(user.top)); 142 PetscCall(PetscFree(user.left)); 143 PetscCall(PetscFree(user.right)); 144 PetscCall(DMDestroy(&user.dm)); 145 PetscCall(PetscFinalize()); 146 return 0; 147 } 148 149 /* -------------------------------------------------------------------- */ 150 /* 151 FormFunction - Evaluates the objective function. 152 153 Input Parameters: 154 . tao - the Tao context 155 . X - input vector 156 . userCtx - optional user-defined context, as set by TaoSetObjective() 157 158 Output Parameters: 159 . fcn - the newly evaluated function 160 */ 161 PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *fcn, void *userCtx) 162 { 163 AppCtx *user = (AppCtx *)userCtx; 164 PetscInt i, j; 165 PetscInt mx = user->mx, my = user->my; 166 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym; 167 PetscReal ft = 0.0; 168 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), area = 0.5 * hx * hy; 169 PetscReal rhx = mx + 1, rhy = my + 1; 170 PetscReal f2, f4, d1, d2, d3, d4, xc, xl, xr, xt, xb; 171 PetscReal **x; 172 Vec localX; 173 174 PetscFunctionBegin; 175 /* Get local mesh boundaries */ 176 PetscCall(DMGetLocalVector(user->dm, &localX)); 177 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 178 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 179 180 /* Scatter ghost points to local vector */ 181 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX)); 182 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX)); 183 184 /* Get pointers to vector data */ 185 PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x)); 186 187 /* Compute function and gradient over the locally owned part of the mesh */ 188 for (j = ys; j < ys + ym; j++) { 189 for (i = xs; i < xs + xm; i++) { 190 xc = x[j][i]; 191 192 if (i == 0) { /* left side */ 193 xl = user->left[j - ys + 1]; 194 } else { 195 xl = x[j][i - 1]; 196 } 197 198 if (j == 0) { /* bottom side */ 199 xb = user->bottom[i - xs + 1]; 200 } else { 201 xb = x[j - 1][i]; 202 } 203 204 if (i + 1 == gxs + gxm) { /* right side */ 205 xr = user->right[j - ys + 1]; 206 } else { 207 xr = x[j][i + 1]; 208 } 209 210 if (j + 1 == gys + gym) { /* top side */ 211 xt = user->top[i - xs + 1]; 212 } else { 213 xt = x[j + 1][i]; 214 } 215 216 d1 = (xc - xl); 217 d2 = (xc - xr); 218 d3 = (xc - xt); 219 d4 = (xc - xb); 220 221 d1 *= rhx; 222 d2 *= rhx; 223 d3 *= rhy; 224 d4 *= rhy; 225 226 f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 227 f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 228 229 ft = ft + (f2 + f4); 230 } 231 } 232 233 /* Compute triangular areas along the border of the domain. */ 234 if (xs == 0) { /* left side */ 235 for (j = ys; j < ys + ym; j++) { 236 d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy; 237 d2 = (user->left[j - ys + 1] - x[j][0]) * rhx; 238 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 239 } 240 } 241 if (ys == 0) { /* bottom side */ 242 for (i = xs; i < xs + xm; i++) { 243 d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx; 244 d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy; 245 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 246 } 247 } 248 if (xs + xm == mx) { /* right side */ 249 for (j = ys; j < ys + ym; j++) { 250 d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx; 251 d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy; 252 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 253 } 254 } 255 if (ys + ym == my) { /* top side */ 256 for (i = xs; i < xs + xm; i++) { 257 d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy; 258 d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx; 259 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 260 } 261 } 262 if (ys == 0 && xs == 0) { 263 d1 = (user->left[0] - user->left[1]) / hy; 264 d2 = (user->bottom[0] - user->bottom[1]) * rhx; 265 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2); 266 } 267 if (ys + ym == my && xs + xm == mx) { 268 d1 = (user->right[ym + 1] - user->right[ym]) * rhy; 269 d2 = (user->top[xm + 1] - user->top[xm]) * rhx; 270 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2); 271 } 272 273 ft = ft * area; 274 PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD)); 275 276 /* Restore vectors */ 277 PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x)); 278 PetscCall(DMRestoreLocalVector(user->dm, &localX)); 279 PetscFunctionReturn(PETSC_SUCCESS); 280 } 281 282 /* -------------------------------------------------------------------- */ 283 /* 284 FormFunctionGradient - Evaluates the function and corresponding gradient. 285 286 Input Parameters: 287 . tao - the Tao context 288 . X - input vector 289 . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient() 290 291 Output Parameters: 292 . fcn - the newly evaluated function 293 . G - vector containing the newly evaluated gradient 294 */ 295 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx) 296 { 297 AppCtx *user = (AppCtx *)userCtx; 298 PetscInt i, j; 299 PetscInt mx = user->mx, my = user->my; 300 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym; 301 PetscReal ft = 0.0; 302 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy; 303 PetscReal rhx = mx + 1, rhy = my + 1; 304 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 305 PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc; 306 PetscReal **g, **x; 307 Vec localX; 308 309 PetscFunctionBegin; 310 /* Get local mesh boundaries */ 311 PetscCall(DMGetLocalVector(user->dm, &localX)); 312 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 313 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 314 315 /* Scatter ghost points to local vector */ 316 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX)); 317 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX)); 318 319 /* Get pointers to vector data */ 320 PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x)); 321 PetscCall(DMDAVecGetArray(user->dm, G, (void **)&g)); 322 323 /* Compute function and gradient over the locally owned part of the mesh */ 324 for (j = ys; j < ys + ym; j++) { 325 for (i = xs; i < xs + xm; i++) { 326 xc = x[j][i]; 327 xlt = xrb = xl = xr = xb = xt = xc; 328 329 if (i == 0) { /* left side */ 330 xl = user->left[j - ys + 1]; 331 xlt = user->left[j - ys + 2]; 332 } else { 333 xl = x[j][i - 1]; 334 } 335 336 if (j == 0) { /* bottom side */ 337 xb = user->bottom[i - xs + 1]; 338 xrb = user->bottom[i - xs + 2]; 339 } else { 340 xb = x[j - 1][i]; 341 } 342 343 if (i + 1 == gxs + gxm) { /* right side */ 344 xr = user->right[j - ys + 1]; 345 xrb = user->right[j - ys]; 346 } else { 347 xr = x[j][i + 1]; 348 } 349 350 if (j + 1 == gys + gym) { /* top side */ 351 xt = user->top[i - xs + 1]; 352 xlt = user->top[i - xs]; 353 } else { 354 xt = x[j + 1][i]; 355 } 356 357 if (i > gxs && j + 1 < gys + gym) xlt = x[j + 1][i - 1]; 358 if (j > gys && i + 1 < gxs + gxm) xrb = x[j - 1][i + 1]; 359 360 d1 = (xc - xl); 361 d2 = (xc - xr); 362 d3 = (xc - xt); 363 d4 = (xc - xb); 364 d5 = (xr - xrb); 365 d6 = (xrb - xb); 366 d7 = (xlt - xl); 367 d8 = (xt - xlt); 368 369 df1dxc = d1 * hydhx; 370 df2dxc = (d1 * hydhx + d4 * hxdhy); 371 df3dxc = d3 * hxdhy; 372 df4dxc = (d2 * hydhx + d3 * hxdhy); 373 df5dxc = d2 * hydhx; 374 df6dxc = d4 * hxdhy; 375 376 d1 *= rhx; 377 d2 *= rhx; 378 d3 *= rhy; 379 d4 *= rhy; 380 d5 *= rhy; 381 d6 *= rhx; 382 d7 *= rhy; 383 d8 *= rhx; 384 385 f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7); 386 f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 387 f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8); 388 f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 389 f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5); 390 f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6); 391 392 ft = ft + (f2 + f4); 393 394 df1dxc /= f1; 395 df2dxc /= f2; 396 df3dxc /= f3; 397 df4dxc /= f4; 398 df5dxc /= f5; 399 df6dxc /= f6; 400 401 g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5; 402 } 403 } 404 405 /* Compute triangular areas along the border of the domain. */ 406 if (xs == 0) { /* left side */ 407 for (j = ys; j < ys + ym; j++) { 408 d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy; 409 d2 = (user->left[j - ys + 1] - x[j][0]) * rhx; 410 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 411 } 412 } 413 if (ys == 0) { /* bottom side */ 414 for (i = xs; i < xs + xm; i++) { 415 d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx; 416 d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy; 417 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 418 } 419 } 420 421 if (xs + xm == mx) { /* right side */ 422 for (j = ys; j < ys + ym; j++) { 423 d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx; 424 d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy; 425 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 426 } 427 } 428 if (ys + ym == my) { /* top side */ 429 for (i = xs; i < xs + xm; i++) { 430 d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy; 431 d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx; 432 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 433 } 434 } 435 436 if (ys == 0 && xs == 0) { 437 d1 = (user->left[0] - user->left[1]) / hy; 438 d2 = (user->bottom[0] - user->bottom[1]) * rhx; 439 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2); 440 } 441 if (ys + ym == my && xs + xm == mx) { 442 d1 = (user->right[ym + 1] - user->right[ym]) * rhy; 443 d2 = (user->top[xm + 1] - user->top[xm]) * rhx; 444 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2); 445 } 446 447 ft = ft * area; 448 PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD)); 449 450 /* Restore vectors */ 451 PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x)); 452 PetscCall(DMDAVecRestoreArray(user->dm, G, (void **)&g)); 453 PetscCall(DMRestoreLocalVector(user->dm, &localX)); 454 PetscCall(PetscLogFlops(67.0 * xm * ym)); 455 PetscFunctionReturn(PETSC_SUCCESS); 456 } 457 458 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *userCtx) 459 { 460 PetscReal fcn; 461 462 PetscFunctionBegin; 463 PetscCall(FormFunctionGradient(tao, X, &fcn, G, userCtx)); 464 PetscFunctionReturn(PETSC_SUCCESS); 465 } 466 467 /* ------------------------------------------------------------------- */ 468 /* 469 FormHessian - Evaluates Hessian matrix. 470 471 Input Parameters: 472 . tao - the Tao context 473 . x - input vector 474 . ptr - optional user-defined context, as set by TaoSetHessian() 475 476 Output Parameters: 477 . H - Hessian matrix 478 . Hpre - optionally different preconditioning matrix 479 480 */ 481 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 482 { 483 AppCtx *user = (AppCtx *)ptr; 484 485 PetscFunctionBegin; 486 /* Evaluate the Hessian entries*/ 487 PetscCall(QuadraticH(user, X, H)); 488 PetscFunctionReturn(PETSC_SUCCESS); 489 } 490 491 /* ------------------------------------------------------------------- */ 492 /* 493 QuadraticH - Evaluates Hessian matrix. 494 495 Input Parameters: 496 . user - user-defined context, as set by TaoSetHessian() 497 . X - input vector 498 499 Output Parameter: 500 . H - Hessian matrix 501 */ 502 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) 503 { 504 PetscInt i, j, k; 505 PetscInt mx = user->mx, my = user->my; 506 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym; 507 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy; 508 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 509 PetscReal hl, hr, ht, hb, hc, htl, hbr; 510 PetscReal **x, v[7]; 511 MatStencil col[7], row; 512 Vec localX; 513 PetscBool assembled; 514 515 PetscFunctionBegin; 516 /* Get local mesh boundaries */ 517 PetscCall(DMGetLocalVector(user->dm, &localX)); 518 519 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 520 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 521 522 /* Scatter ghost points to local vector */ 523 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX)); 524 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX)); 525 526 /* Get pointers to vector data */ 527 PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x)); 528 529 /* Initialize matrix entries to zero */ 530 PetscCall(MatAssembled(Hessian, &assembled)); 531 if (assembled) PetscCall(MatZeroEntries(Hessian)); 532 533 /* Set various matrix options */ 534 PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE)); 535 536 /* Compute Hessian over the locally owned part of the mesh */ 537 for (j = ys; j < ys + ym; j++) { 538 for (i = xs; i < xs + xm; i++) { 539 xc = x[j][i]; 540 xlt = xrb = xl = xr = xb = xt = xc; 541 542 /* Left side */ 543 if (i == 0) { 544 xl = user->left[j - ys + 1]; 545 xlt = user->left[j - ys + 2]; 546 } else { 547 xl = x[j][i - 1]; 548 } 549 550 if (j == 0) { 551 xb = user->bottom[i - xs + 1]; 552 xrb = user->bottom[i - xs + 2]; 553 } else { 554 xb = x[j - 1][i]; 555 } 556 557 if (i + 1 == mx) { 558 xr = user->right[j - ys + 1]; 559 xrb = user->right[j - ys]; 560 } else { 561 xr = x[j][i + 1]; 562 } 563 564 if (j + 1 == my) { 565 xt = user->top[i - xs + 1]; 566 xlt = user->top[i - xs]; 567 } else { 568 xt = x[j + 1][i]; 569 } 570 571 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; 572 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; 573 574 d1 = (xc - xl) / hx; 575 d2 = (xc - xr) / hx; 576 d3 = (xc - xt) / hy; 577 d4 = (xc - xb) / hy; 578 d5 = (xrb - xr) / hy; 579 d6 = (xrb - xb) / hx; 580 d7 = (xlt - xl) / hy; 581 d8 = (xlt - xt) / hx; 582 583 f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7); 584 f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 585 f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8); 586 f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 587 f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5); 588 f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6); 589 590 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 591 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 592 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 593 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 594 595 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 596 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 597 598 hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2 * d2 * d3) / (f4 * f4 * f4); 599 600 hl /= 2.0; 601 hr /= 2.0; 602 ht /= 2.0; 603 hb /= 2.0; 604 hbr /= 2.0; 605 htl /= 2.0; 606 hc /= 2.0; 607 608 row.j = j; 609 row.i = i; 610 k = 0; 611 if (j > 0) { 612 v[k] = hb; 613 col[k].j = j - 1; 614 col[k].i = i; 615 k++; 616 } 617 618 if (j > 0 && i < mx - 1) { 619 v[k] = hbr; 620 col[k].j = j - 1; 621 col[k].i = i + 1; 622 k++; 623 } 624 625 if (i > 0) { 626 v[k] = hl; 627 col[k].j = j; 628 col[k].i = i - 1; 629 k++; 630 } 631 632 v[k] = hc; 633 col[k].j = j; 634 col[k].i = i; 635 k++; 636 637 if (i < mx - 1) { 638 v[k] = hr; 639 col[k].j = j; 640 col[k].i = i + 1; 641 k++; 642 } 643 644 if (i > 0 && j < my - 1) { 645 v[k] = htl; 646 col[k].j = j + 1; 647 col[k].i = i - 1; 648 k++; 649 } 650 651 if (j < my - 1) { 652 v[k] = ht; 653 col[k].j = j + 1; 654 col[k].i = i; 655 k++; 656 } 657 658 /* 659 Set matrix values using local numbering, which was defined 660 earlier, in the main routine. 661 */ 662 PetscCall(MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES)); 663 } 664 } 665 666 PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x)); 667 PetscCall(DMRestoreLocalVector(user->dm, &localX)); 668 669 PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY)); 670 PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY)); 671 672 PetscCall(PetscLogFlops(199.0 * xm * ym)); 673 PetscFunctionReturn(PETSC_SUCCESS); 674 } 675 676 /* ------------------------------------------------------------------- */ 677 /* 678 MSA_BoundaryConditions - Calculates the boundary conditions for 679 the region. 680 681 Input Parameter: 682 . user - user-defined application context 683 684 Output Parameter: 685 . user - user-defined application context 686 */ 687 static PetscErrorCode MSA_BoundaryConditions(AppCtx *user) 688 { 689 PetscInt i, j, k, limit = 0, maxits = 5; 690 PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym; 691 PetscInt mx = user->mx, my = user->my; 692 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 693 PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10; 694 PetscReal fnorm, det, hx, hy, xt = 0, yt = 0; 695 PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 696 PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5; 697 PetscReal *boundary; 698 PetscBool flg; 699 700 PetscFunctionBegin; 701 /* Get local mesh boundaries */ 702 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 703 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 704 705 bsize = xm + 2; 706 lsize = ym + 2; 707 rsize = ym + 2; 708 tsize = xm + 2; 709 710 PetscCall(PetscMalloc1(bsize, &user->bottom)); 711 PetscCall(PetscMalloc1(tsize, &user->top)); 712 PetscCall(PetscMalloc1(lsize, &user->left)); 713 PetscCall(PetscMalloc1(rsize, &user->right)); 714 715 hx = (r - l) / (mx + 1); 716 hy = (t - b) / (my + 1); 717 718 for (j = 0; j < 4; j++) { 719 if (j == 0) { 720 yt = b; 721 xt = l + hx * xs; 722 limit = bsize; 723 boundary = user->bottom; 724 } else if (j == 1) { 725 yt = t; 726 xt = l + hx * xs; 727 limit = tsize; 728 boundary = user->top; 729 } else if (j == 2) { 730 yt = b + hy * ys; 731 xt = l; 732 limit = lsize; 733 boundary = user->left; 734 } else { /* if (j==3) */ 735 yt = b + hy * ys; 736 xt = r; 737 limit = rsize; 738 boundary = user->right; 739 } 740 741 for (i = 0; i < limit; i++) { 742 u1 = xt; 743 u2 = -yt; 744 for (k = 0; k < maxits; k++) { 745 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 746 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 747 fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2); 748 if (fnorm <= tol) break; 749 njac11 = one + u2 * u2 - u1 * u1; 750 njac12 = two * u1 * u2; 751 njac21 = -two * u1 * u2; 752 njac22 = -one - u1 * u1 + u2 * u2; 753 det = njac11 * njac22 - njac21 * njac12; 754 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 755 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 756 } 757 758 boundary[i] = u1 * u1 - u2 * u2; 759 if (j == 0 || j == 1) { 760 xt = xt + hx; 761 } else { /* if (j==2 || j==3) */ 762 yt = yt + hy; 763 } 764 } 765 } 766 767 /* Scale the boundary if desired */ 768 if (1 == 1) { 769 PetscReal scl = 1.0; 770 771 PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg)); 772 if (flg) { 773 for (i = 0; i < bsize; i++) user->bottom[i] *= scl; 774 } 775 776 PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg)); 777 if (flg) { 778 for (i = 0; i < tsize; i++) user->top[i] *= scl; 779 } 780 781 PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg)); 782 if (flg) { 783 for (i = 0; i < rsize; i++) user->right[i] *= scl; 784 } 785 786 PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg)); 787 if (flg) { 788 for (i = 0; i < lsize; i++) user->left[i] *= scl; 789 } 790 } 791 PetscFunctionReturn(PETSC_SUCCESS); 792 } 793 794 /* ------------------------------------------------------------------- */ 795 /* 796 MSA_InitialPoint - Calculates the initial guess in one of three ways. 797 798 Input Parameters: 799 . user - user-defined application context 800 . X - vector for initial guess 801 802 Output Parameters: 803 . X - newly computed initial guess 804 */ 805 static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X) 806 { 807 PetscInt start2 = -1, i, j; 808 PetscReal start1 = 0; 809 PetscBool flg1, flg2; 810 811 PetscFunctionBegin; 812 PetscCall(PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1)); 813 PetscCall(PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2)); 814 815 if (flg1) { /* The zero vector is reasonable */ 816 PetscCall(VecSet(X, start1)); 817 } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */ 818 PetscRandom rctx; 819 PetscReal np5 = -0.5; 820 821 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 822 PetscCall(VecSetRandom(X, rctx)); 823 PetscCall(PetscRandomDestroy(&rctx)); 824 PetscCall(VecShift(X, np5)); 825 } else { /* Take an average of the boundary conditions */ 826 PetscInt xs, xm, ys, ym; 827 PetscInt mx = user->mx, my = user->my; 828 PetscReal **x; 829 830 /* Get local mesh boundaries */ 831 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 832 833 /* Get pointers to vector data */ 834 PetscCall(DMDAVecGetArray(user->dm, X, (void **)&x)); 835 836 /* Perform local computations */ 837 for (j = ys; j < ys + ym; j++) { 838 for (i = xs; i < xs + xm; i++) x[j][i] = (((j + 1) * user->bottom[i - xs + 1] + (my - j + 1) * user->top[i - xs + 1]) / (my + 2) + ((i + 1) * user->left[j - ys + 1] + (mx - i + 1) * user->right[j - ys + 1]) / (mx + 2)) / 2.0; 839 } 840 PetscCall(DMDAVecRestoreArray(user->dm, X, (void **)&x)); 841 PetscCall(PetscLogFlops(9.0 * xm * ym)); 842 } 843 PetscFunctionReturn(PETSC_SUCCESS); 844 } 845 846 /*-----------------------------------------------------------------------*/ 847 PetscErrorCode My_Monitor(Tao tao, void *ctx) 848 { 849 Vec X; 850 851 PetscFunctionBegin; 852 PetscCall(TaoGetSolution(tao, &X)); 853 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD)); 854 PetscFunctionReturn(PETSC_SUCCESS); 855 } 856 857 /*TEST 858 859 build: 860 requires: !complex 861 862 test: 863 args: -tao_monitor_short -tao_type lmvm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3 864 requires: !single 865 866 test: 867 suffix: 2 868 nsize: 2 869 args: -tao_monitor_short -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4 870 filter: grep -v "nls ksp" 871 requires: !single 872 873 test: 874 suffix: 2_snes 875 nsize: 2 876 args: -tao_monitor_short -tao_type snes -ksp_converged_maxits -ksp_max_it 15 -snes_atol 1.e-4 877 filter: grep -v "nls ksp" 878 requires: !single 879 880 test: 881 suffix: 3 882 nsize: 3 883 args: -tao_monitor_short -tao_type cg -tao_cg_type fr -da_grid_x 10 -da_grid_y 10 -tao_gatol 1.e-3 884 requires: !single 885 886 test: 887 suffix: 3_snes 888 nsize: 3 889 args: -tao_monitor_short -tao_type snes -snes_type ncg -snes_ncg_type fr -da_grid_x 10 -da_grid_y 10 -snes_atol 1.e-4 890 requires: !single 891 892 test: 893 suffix: 4_snes_ngmres 894 args: -tao_type snes -snes_type ngmres -npc_snes_type ncg -da_grid_x 10 -da_grid_y 10 -snes_atol 1.e-4 -npc_snes_ncg_type fr -snes_converged_reason -snes_monitor ::ascii_info_detail -snes_ngmres_monitor -snes_ngmres_select_type {{linesearch difference}separate output} 895 requires: !single 896 897 test: 898 suffix: 5 899 nsize: 2 900 args: -tao_monitor_short -tao_type bmrm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3 901 requires: !single 902 903 TEST*/ 904