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