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