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 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 . flg - flag indicating matrix structure 480 481 */ 482 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 483 { 484 AppCtx *user = (AppCtx *)ptr; 485 486 PetscFunctionBegin; 487 /* Evaluate the Hessian entries*/ 488 PetscCall(QuadraticH(user, X, H)); 489 PetscFunctionReturn(PETSC_SUCCESS); 490 } 491 492 /* ------------------------------------------------------------------- */ 493 /* 494 QuadraticH - Evaluates Hessian matrix. 495 496 Input Parameters: 497 . user - user-defined context, as set by TaoSetHessian() 498 . X - input vector 499 500 Output Parameter: 501 . H - Hessian matrix 502 */ 503 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) 504 { 505 PetscInt i, j, k; 506 PetscInt mx = user->mx, my = user->my; 507 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym; 508 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy; 509 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 510 PetscReal hl, hr, ht, hb, hc, htl, hbr; 511 PetscReal **x, v[7]; 512 MatStencil col[7], row; 513 Vec localX; 514 PetscBool assembled; 515 516 PetscFunctionBegin; 517 /* Get local mesh boundaries */ 518 PetscCall(DMGetLocalVector(user->dm, &localX)); 519 520 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 521 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 522 523 /* Scatter ghost points to local vector */ 524 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX)); 525 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX)); 526 527 /* Get pointers to vector data */ 528 PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x)); 529 530 /* Initialize matrix entries to zero */ 531 PetscCall(MatAssembled(Hessian, &assembled)); 532 if (assembled) PetscCall(MatZeroEntries(Hessian)); 533 534 /* Set various matrix options */ 535 PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE)); 536 537 /* Compute Hessian over the locally owned part of the mesh */ 538 for (j = ys; j < ys + ym; j++) { 539 for (i = xs; i < xs + xm; i++) { 540 xc = x[j][i]; 541 xlt = xrb = xl = xr = xb = xt = xc; 542 543 /* Left side */ 544 if (i == 0) { 545 xl = user->left[j - ys + 1]; 546 xlt = user->left[j - ys + 2]; 547 } else { 548 xl = x[j][i - 1]; 549 } 550 551 if (j == 0) { 552 xb = user->bottom[i - xs + 1]; 553 xrb = user->bottom[i - xs + 2]; 554 } else { 555 xb = x[j - 1][i]; 556 } 557 558 if (i + 1 == mx) { 559 xr = user->right[j - ys + 1]; 560 xrb = user->right[j - ys]; 561 } else { 562 xr = x[j][i + 1]; 563 } 564 565 if (j + 1 == my) { 566 xt = user->top[i - xs + 1]; 567 xlt = user->top[i - xs]; 568 } else { 569 xt = x[j + 1][i]; 570 } 571 572 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; 573 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; 574 575 d1 = (xc - xl) / hx; 576 d2 = (xc - xr) / hx; 577 d3 = (xc - xt) / hy; 578 d4 = (xc - xb) / hy; 579 d5 = (xrb - xr) / hy; 580 d6 = (xrb - xb) / hx; 581 d7 = (xlt - xl) / hy; 582 d8 = (xlt - xt) / hx; 583 584 f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7); 585 f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4); 586 f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8); 587 f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2); 588 f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5); 589 f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6); 590 591 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 592 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 593 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 594 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 595 596 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 597 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 598 599 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); 600 601 hl /= 2.0; 602 hr /= 2.0; 603 ht /= 2.0; 604 hb /= 2.0; 605 hbr /= 2.0; 606 htl /= 2.0; 607 hc /= 2.0; 608 609 row.j = j; 610 row.i = i; 611 k = 0; 612 if (j > 0) { 613 v[k] = hb; 614 col[k].j = j - 1; 615 col[k].i = i; 616 k++; 617 } 618 619 if (j > 0 && i < mx - 1) { 620 v[k] = hbr; 621 col[k].j = j - 1; 622 col[k].i = i + 1; 623 k++; 624 } 625 626 if (i > 0) { 627 v[k] = hl; 628 col[k].j = j; 629 col[k].i = i - 1; 630 k++; 631 } 632 633 v[k] = hc; 634 col[k].j = j; 635 col[k].i = i; 636 k++; 637 638 if (i < mx - 1) { 639 v[k] = hr; 640 col[k].j = j; 641 col[k].i = i + 1; 642 k++; 643 } 644 645 if (i > 0 && j < my - 1) { 646 v[k] = htl; 647 col[k].j = j + 1; 648 col[k].i = i - 1; 649 k++; 650 } 651 652 if (j < my - 1) { 653 v[k] = ht; 654 col[k].j = j + 1; 655 col[k].i = i; 656 k++; 657 } 658 659 /* 660 Set matrix values using local numbering, which was defined 661 earlier, in the main routine. 662 */ 663 PetscCall(MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES)); 664 } 665 } 666 667 PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x)); 668 PetscCall(DMRestoreLocalVector(user->dm, &localX)); 669 670 PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY)); 671 PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY)); 672 673 PetscCall(PetscLogFlops(199.0 * xm * ym)); 674 PetscFunctionReturn(PETSC_SUCCESS); 675 } 676 677 /* ------------------------------------------------------------------- */ 678 /* 679 MSA_BoundaryConditions - Calculates the boundary conditions for 680 the region. 681 682 Input Parameter: 683 . user - user-defined application context 684 685 Output Parameter: 686 . user - user-defined application context 687 */ 688 static PetscErrorCode MSA_BoundaryConditions(AppCtx *user) 689 { 690 PetscInt i, j, k, limit = 0, maxits = 5; 691 PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym; 692 PetscInt mx = user->mx, my = user->my; 693 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 694 PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10; 695 PetscReal fnorm, det, hx, hy, xt = 0, yt = 0; 696 PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 697 PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5; 698 PetscReal *boundary; 699 PetscBool flg; 700 701 PetscFunctionBegin; 702 /* Get local mesh boundaries */ 703 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 704 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 705 706 bsize = xm + 2; 707 lsize = ym + 2; 708 rsize = ym + 2; 709 tsize = xm + 2; 710 711 PetscCall(PetscMalloc1(bsize, &user->bottom)); 712 PetscCall(PetscMalloc1(tsize, &user->top)); 713 PetscCall(PetscMalloc1(lsize, &user->left)); 714 PetscCall(PetscMalloc1(rsize, &user->right)); 715 716 hx = (r - l) / (mx + 1); 717 hy = (t - b) / (my + 1); 718 719 for (j = 0; j < 4; j++) { 720 if (j == 0) { 721 yt = b; 722 xt = l + hx * xs; 723 limit = bsize; 724 boundary = user->bottom; 725 } else if (j == 1) { 726 yt = t; 727 xt = l + hx * xs; 728 limit = tsize; 729 boundary = user->top; 730 } else if (j == 2) { 731 yt = b + hy * ys; 732 xt = l; 733 limit = lsize; 734 boundary = user->left; 735 } else { /* if (j==3) */ 736 yt = b + hy * ys; 737 xt = r; 738 limit = rsize; 739 boundary = user->right; 740 } 741 742 for (i = 0; i < limit; i++) { 743 u1 = xt; 744 u2 = -yt; 745 for (k = 0; k < maxits; k++) { 746 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 747 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 748 fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2); 749 if (fnorm <= tol) break; 750 njac11 = one + u2 * u2 - u1 * u1; 751 njac12 = two * u1 * u2; 752 njac21 = -two * u1 * u2; 753 njac22 = -one - u1 * u1 + u2 * u2; 754 det = njac11 * njac22 - njac21 * njac12; 755 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 756 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 757 } 758 759 boundary[i] = u1 * u1 - u2 * u2; 760 if (j == 0 || j == 1) { 761 xt = xt + hx; 762 } else { /* if (j==2 || j==3) */ 763 yt = yt + hy; 764 } 765 } 766 } 767 768 /* Scale the boundary if desired */ 769 if (1 == 1) { 770 PetscReal scl = 1.0; 771 772 PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg)); 773 if (flg) { 774 for (i = 0; i < bsize; i++) user->bottom[i] *= scl; 775 } 776 777 PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg)); 778 if (flg) { 779 for (i = 0; i < tsize; i++) user->top[i] *= scl; 780 } 781 782 PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg)); 783 if (flg) { 784 for (i = 0; i < rsize; i++) user->right[i] *= scl; 785 } 786 787 PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg)); 788 if (flg) { 789 for (i = 0; i < lsize; i++) user->left[i] *= scl; 790 } 791 } 792 PetscFunctionReturn(PETSC_SUCCESS); 793 } 794 795 /* ------------------------------------------------------------------- */ 796 /* 797 MSA_InitialPoint - Calculates the initial guess in one of three ways. 798 799 Input Parameters: 800 . user - user-defined application context 801 . X - vector for initial guess 802 803 Output Parameters: 804 . X - newly computed initial guess 805 */ 806 static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X) 807 { 808 PetscInt start2 = -1, i, j; 809 PetscReal start1 = 0; 810 PetscBool flg1, flg2; 811 812 PetscFunctionBegin; 813 PetscCall(PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1)); 814 PetscCall(PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2)); 815 816 if (flg1) { /* The zero vector is reasonable */ 817 PetscCall(VecSet(X, start1)); 818 } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */ 819 PetscRandom rctx; 820 PetscReal np5 = -0.5; 821 822 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 823 PetscCall(VecSetRandom(X, rctx)); 824 PetscCall(PetscRandomDestroy(&rctx)); 825 PetscCall(VecShift(X, np5)); 826 } else { /* Take an average of the boundary conditions */ 827 PetscInt xs, xm, ys, ym; 828 PetscInt mx = user->mx, my = user->my; 829 PetscReal **x; 830 831 /* Get local mesh boundaries */ 832 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 833 834 /* Get pointers to vector data */ 835 PetscCall(DMDAVecGetArray(user->dm, X, (void **)&x)); 836 837 /* Perform local computations */ 838 for (j = ys; j < ys + ym; j++) { 839 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; 840 } 841 PetscCall(DMDAVecRestoreArray(user->dm, X, (void **)&x)); 842 PetscCall(PetscLogFlops(9.0 * xm * ym)); 843 } 844 PetscFunctionReturn(PETSC_SUCCESS); 845 } 846 847 /*-----------------------------------------------------------------------*/ 848 PetscErrorCode My_Monitor(Tao tao, void *ctx) 849 { 850 Vec X; 851 852 PetscFunctionBegin; 853 PetscCall(TaoGetSolution(tao, &X)); 854 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD)); 855 PetscFunctionReturn(PETSC_SUCCESS); 856 } 857 858 /*TEST 859 860 build: 861 requires: !complex 862 863 test: 864 args: -tao_monitor_short -tao_type lmvm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3 865 requires: !single 866 867 test: 868 suffix: 2 869 nsize: 2 870 args: -tao_monitor_short -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4 871 filter: grep -v "nls ksp" 872 requires: !single 873 874 test: 875 suffix: 2_snes 876 nsize: 2 877 args: -tao_monitor_short -tao_type snes -ksp_converged_maxits -ksp_max_it 15 -snes_atol 1.e-4 878 filter: grep -v "nls ksp" 879 requires: !single 880 881 test: 882 suffix: 3 883 nsize: 3 884 args: -tao_monitor_short -tao_type cg -tao_cg_type fr -da_grid_x 10 -da_grid_y 10 -tao_gatol 1.e-3 885 requires: !single 886 887 test: 888 suffix: 3_snes 889 nsize: 3 890 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 891 requires: !single 892 893 test: 894 suffix: 4_snes_ngmres 895 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} 896 requires: !single 897 898 test: 899 suffix: 5 900 nsize: 2 901 args: -tao_monitor_short -tao_type bmrm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3 902 requires: !single 903 904 TEST*/ 905