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