1 /* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */ 2 3 /* Include "petsctao.h" so we can use TAO solvers. */ 4 #include <petsctao.h> 5 6 static char help[] = "This example demonstrates use of the TAO package to \n\ 7 solve an unconstrained minimization problem. This example is based on a \n\ 8 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\ 9 boundary values along the edges of the domain, the objective is to find the\n\ 10 surface with the minimal area that satisfies the boundary conditions.\n\ 11 The command line options are:\n\ 12 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ 13 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ 14 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n"; 15 16 /* 17 User-defined application context - contains data needed by the 18 application-provided call-back routines, FormFunctionGradient() 19 and FormHessian(). 20 */ 21 typedef struct { 22 PetscInt mx, my; /* discretization in x, y directions */ 23 PetscReal *bottom, *top, *left, *right; /* boundary values */ 24 Mat H; 25 } AppCtx; 26 27 /* -------- User-defined Routines --------- */ 28 29 static PetscErrorCode MSA_BoundaryConditions(AppCtx *); 30 static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec); 31 static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat); 32 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *); 33 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *); 34 35 int main(int argc, char **argv) 36 { 37 PetscInt N; /* Size of vector */ 38 PetscMPIInt size; /* Number of processors */ 39 Vec x; /* solution, gradient vectors */ 40 KSP ksp; /* PETSc Krylov subspace method */ 41 PetscBool flg; /* A return value when checking for user options */ 42 Tao tao; /* Tao solver context */ 43 AppCtx user; /* user-defined work context */ 44 45 /* Initialize TAO,PETSc */ 46 PetscFunctionBeginUser; 47 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 48 49 PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size)); 50 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors"); 51 52 /* Specify default dimension of the problem */ 53 user.mx = 4; 54 user.my = 4; 55 56 /* Check for any command line arguments that override defaults */ 57 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg)); 58 PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg)); 59 60 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Minimum Surface Area Problem -----\n")); 61 PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my)); 62 63 /* Calculate any derived values from parameters */ 64 N = user.mx * user.my; 65 66 /* Create TAO solver and set desired solution method */ 67 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao)); 68 PetscCall(TaoSetType(tao, TAOLMVM)); 69 70 /* Initialize minsurf application data structure for use in the function evaluations */ 71 PetscCall(MSA_BoundaryConditions(&user)); /* Application specific routine */ 72 73 /* 74 Create a vector to hold the variables. Compute an initial solution. 75 Set this vector, which will be used by TAO. 76 */ 77 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x)); 78 PetscCall(MSA_InitialPoint(&user, x)); /* Application specific routine */ 79 PetscCall(TaoSetSolution(tao, x)); /* A TAO routine */ 80 81 /* Provide TAO routines for function, gradient, and Hessian evaluation */ 82 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user)); 83 84 /* Create a matrix data structure to store the Hessian. This structure will be used by TAO */ 85 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 7, NULL, &(user.H))); 86 PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE)); 87 PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user)); 88 89 /* Check for any TAO command line options */ 90 PetscCall(TaoSetFromOptions(tao)); 91 92 /* Limit the number of iterations in the KSP linear solver */ 93 PetscCall(TaoGetKSP(tao, &ksp)); 94 if (ksp) PetscCall(KSPSetTolerances(ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, user.mx * user.my)); 95 96 /* SOLVE THE APPLICATION */ 97 PetscCall(TaoSolve(tao)); 98 99 PetscCall(TaoDestroy(&tao)); 100 PetscCall(VecDestroy(&x)); 101 PetscCall(MatDestroy(&user.H)); 102 PetscCall(PetscFree(user.bottom)); 103 PetscCall(PetscFree(user.top)); 104 PetscCall(PetscFree(user.left)); 105 PetscCall(PetscFree(user.right)); 106 107 PetscCall(PetscFinalize()); 108 return 0; 109 } 110 111 /* -------------------------------------------------------------------- */ 112 113 /* FormFunctionGradient - Evaluates function and corresponding gradient. 114 115 Input Parameters: 116 . tao - the Tao context 117 . X - input vector 118 . userCtx - optional user-defined context, as set by TaoSetFunctionGradient() 119 120 Output Parameters: 121 . fcn - the newly evaluated function 122 . G - vector containing the newly evaluated gradient 123 */ 124 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx) 125 { 126 AppCtx *user = (AppCtx *)userCtx; 127 PetscInt i, j, row; 128 PetscInt mx = user->mx, my = user->my; 129 PetscReal rhx = mx + 1, rhy = my + 1; 130 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy, ft = 0; 131 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 132 PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc; 133 PetscReal zero = 0.0; 134 PetscScalar *g; 135 const PetscScalar *x; 136 137 PetscFunctionBeginUser; 138 PetscCall(VecSet(G, zero)); 139 140 PetscCall(VecGetArrayRead(X, &x)); 141 PetscCall(VecGetArray(G, &g)); 142 143 /* Compute function over the locally owned part of the mesh */ 144 for (j = 0; j < my; j++) { 145 for (i = 0; i < mx; i++) { 146 row = (j)*mx + (i); 147 xc = x[row]; 148 xlt = xrb = xl = xr = xb = xt = xc; 149 if (i == 0) { /* left side */ 150 xl = user->left[j + 1]; 151 xlt = user->left[j + 2]; 152 } else { 153 xl = x[row - 1]; 154 } 155 156 if (j == 0) { /* bottom side */ 157 xb = user->bottom[i + 1]; 158 xrb = user->bottom[i + 2]; 159 } else { 160 xb = x[row - mx]; 161 } 162 163 if (i + 1 == mx) { /* right side */ 164 xr = user->right[j + 1]; 165 xrb = user->right[j]; 166 } else { 167 xr = x[row + 1]; 168 } 169 170 if (j + 1 == 0 + my) { /* top side */ 171 xt = user->top[i + 1]; 172 xlt = user->top[i]; 173 } else { 174 xt = x[row + mx]; 175 } 176 177 if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx]; 178 if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx]; 179 180 d1 = (xc - xl); 181 d2 = (xc - xr); 182 d3 = (xc - xt); 183 d4 = (xc - xb); 184 d5 = (xr - xrb); 185 d6 = (xrb - xb); 186 d7 = (xlt - xl); 187 d8 = (xt - xlt); 188 189 df1dxc = d1 * hydhx; 190 df2dxc = (d1 * hydhx + d4 * hxdhy); 191 df3dxc = d3 * hxdhy; 192 df4dxc = (d2 * hydhx + d3 * hxdhy); 193 df5dxc = d2 * hydhx; 194 df6dxc = d4 * hxdhy; 195 196 d1 *= rhx; 197 d2 *= rhx; 198 d3 *= rhy; 199 d4 *= rhy; 200 d5 *= rhy; 201 d6 *= rhx; 202 d7 *= rhy; 203 d8 *= rhx; 204 205 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 206 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 207 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 208 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 209 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 210 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 211 212 ft = ft + (f2 + f4); 213 214 df1dxc /= f1; 215 df2dxc /= f2; 216 df3dxc /= f3; 217 df4dxc /= f4; 218 df5dxc /= f5; 219 df6dxc /= f6; 220 221 g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0; 222 } 223 } 224 225 for (j = 0; j < my; j++) { /* left side */ 226 d3 = (user->left[j + 1] - user->left[j + 2]) * rhy; 227 d2 = (user->left[j + 1] - x[j * mx]) * rhx; 228 ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 229 } 230 231 for (i = 0; i < mx; i++) { /* bottom */ 232 d2 = (user->bottom[i + 1] - user->bottom[i + 2]) * rhx; 233 d3 = (user->bottom[i + 1] - x[i]) * rhy; 234 ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 235 } 236 237 for (j = 0; j < my; j++) { /* right side */ 238 d1 = (x[(j + 1) * mx - 1] - user->right[j + 1]) * rhx; 239 d4 = (user->right[j] - user->right[j + 1]) * rhy; 240 ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 241 } 242 243 for (i = 0; i < mx; i++) { /* top side */ 244 d1 = (x[(my - 1) * mx + i] - user->top[i + 1]) * rhy; 245 d4 = (user->top[i + 1] - user->top[i]) * rhx; 246 ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 247 } 248 249 /* Bottom left corner */ 250 d1 = (user->left[0] - user->left[1]) * rhy; 251 d2 = (user->bottom[0] - user->bottom[1]) * rhx; 252 ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2); 253 254 /* Top right corner */ 255 d1 = (user->right[my + 1] - user->right[my]) * rhy; 256 d2 = (user->top[mx + 1] - user->top[mx]) * rhx; 257 ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2); 258 259 (*fcn) = ft * area; 260 261 /* Restore vectors */ 262 PetscCall(VecRestoreArrayRead(X, &x)); 263 PetscCall(VecRestoreArray(G, &g)); 264 PetscCall(PetscLogFlops(67.0 * mx * my)); 265 PetscFunctionReturn(PETSC_SUCCESS); 266 } 267 268 /* ------------------------------------------------------------------- */ 269 /* 270 FormHessian - Evaluates the Hessian matrix. 271 272 Input Parameters: 273 . tao - the Tao context 274 . x - input vector 275 . ptr - optional user-defined context, as set by TaoSetHessian() 276 277 Output Parameters: 278 . H - Hessian matrix 279 . Hpre - optionally different preconditioning matrix 280 . flg - flag indicating matrix structure 281 282 */ 283 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 284 { 285 AppCtx *user = (AppCtx *)ptr; 286 287 PetscFunctionBeginUser; 288 /* Evaluate the Hessian entries*/ 289 PetscCall(QuadraticH(user, X, H)); 290 PetscFunctionReturn(PETSC_SUCCESS); 291 } 292 293 /* ------------------------------------------------------------------- */ 294 /* 295 QuadraticH - Evaluates the Hessian matrix. 296 297 Input Parameters: 298 . user - user-defined context, as set by TaoSetHessian() 299 . X - input vector 300 301 Output Parameter: 302 . H - Hessian matrix 303 */ 304 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) 305 { 306 PetscInt i, j, k, row; 307 PetscInt mx = user->mx, my = user->my; 308 PetscInt col[7]; 309 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy; 310 PetscReal rhx = mx + 1, rhy = my + 1; 311 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 312 PetscReal hl, hr, ht, hb, hc, htl, hbr; 313 const PetscScalar *x; 314 PetscReal v[7]; 315 316 PetscFunctionBeginUser; 317 /* Get pointers to vector data */ 318 PetscCall(VecGetArrayRead(X, &x)); 319 320 /* Initialize matrix entries to zero */ 321 PetscCall(MatZeroEntries(Hessian)); 322 323 /* Set various matrix options */ 324 PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE)); 325 326 /* Compute Hessian over the locally owned part of the mesh */ 327 for (i = 0; i < mx; i++) { 328 for (j = 0; j < my; j++) { 329 row = (j)*mx + (i); 330 331 xc = x[row]; 332 xlt = xrb = xl = xr = xb = xt = xc; 333 334 /* Left side */ 335 if (i == 0) { 336 xl = user->left[j + 1]; 337 xlt = user->left[j + 2]; 338 } else { 339 xl = x[row - 1]; 340 } 341 342 if (j == 0) { 343 xb = user->bottom[i + 1]; 344 xrb = user->bottom[i + 2]; 345 } else { 346 xb = x[row - mx]; 347 } 348 349 if (i + 1 == mx) { 350 xr = user->right[j + 1]; 351 xrb = user->right[j]; 352 } else { 353 xr = x[row + 1]; 354 } 355 356 if (j + 1 == my) { 357 xt = user->top[i + 1]; 358 xlt = user->top[i]; 359 } else { 360 xt = x[row + mx]; 361 } 362 363 if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx]; 364 if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx]; 365 366 d1 = (xc - xl) * rhx; 367 d2 = (xc - xr) * rhx; 368 d3 = (xc - xt) * rhy; 369 d4 = (xc - xb) * rhy; 370 d5 = (xrb - xr) * rhy; 371 d6 = (xrb - xb) * rhx; 372 d7 = (xlt - xl) * rhy; 373 d8 = (xlt - xt) * rhx; 374 375 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 376 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 377 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 378 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 379 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 380 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 381 382 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 383 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 384 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 385 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 386 387 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 388 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 389 390 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); 391 392 hl *= 0.5; 393 hr *= 0.5; 394 ht *= 0.5; 395 hb *= 0.5; 396 hbr *= 0.5; 397 htl *= 0.5; 398 hc *= 0.5; 399 400 k = 0; 401 if (j > 0) { 402 v[k] = hb; 403 col[k] = row - mx; 404 k++; 405 } 406 407 if (j > 0 && i < mx - 1) { 408 v[k] = hbr; 409 col[k] = row - mx + 1; 410 k++; 411 } 412 413 if (i > 0) { 414 v[k] = hl; 415 col[k] = row - 1; 416 k++; 417 } 418 419 v[k] = hc; 420 col[k] = row; 421 k++; 422 423 if (i < mx - 1) { 424 v[k] = hr; 425 col[k] = row + 1; 426 k++; 427 } 428 429 if (i > 0 && j < my - 1) { 430 v[k] = htl; 431 col[k] = row + mx - 1; 432 k++; 433 } 434 435 if (j < my - 1) { 436 v[k] = ht; 437 col[k] = row + mx; 438 k++; 439 } 440 441 /* 442 Set matrix values using local numbering, which was defined 443 earlier, in the main routine. 444 */ 445 PetscCall(MatSetValues(Hessian, 1, &row, k, col, v, INSERT_VALUES)); 446 } 447 } 448 449 /* Restore vectors */ 450 PetscCall(VecRestoreArrayRead(X, &x)); 451 452 /* Assemble the matrix */ 453 PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY)); 454 PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY)); 455 456 PetscCall(PetscLogFlops(199.0 * mx * my)); 457 PetscFunctionReturn(PETSC_SUCCESS); 458 } 459 460 /* ------------------------------------------------------------------- */ 461 /* 462 MSA_BoundaryConditions - Calculates the boundary conditions for 463 the region. 464 465 Input Parameter: 466 . user - user-defined application context 467 468 Output Parameter: 469 . user - user-defined application context 470 */ 471 static PetscErrorCode MSA_BoundaryConditions(AppCtx *user) 472 { 473 PetscInt i, j, k, limit = 0; 474 PetscInt maxits = 5; 475 PetscInt mx = user->mx, my = user->my; 476 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 477 PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10; 478 PetscReal fnorm, det, hx, hy, xt = 0, yt = 0; 479 PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 480 PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5; 481 PetscReal *boundary; 482 483 PetscFunctionBeginUser; 484 bsize = mx + 2; 485 lsize = my + 2; 486 rsize = my + 2; 487 tsize = mx + 2; 488 489 PetscCall(PetscMalloc1(bsize, &user->bottom)); 490 PetscCall(PetscMalloc1(tsize, &user->top)); 491 PetscCall(PetscMalloc1(lsize, &user->left)); 492 PetscCall(PetscMalloc1(rsize, &user->right)); 493 494 hx = (r - l) / (mx + 1); 495 hy = (t - b) / (my + 1); 496 497 for (j = 0; j < 4; j++) { 498 if (j == 0) { 499 yt = b; 500 xt = l; 501 limit = bsize; 502 boundary = user->bottom; 503 } else if (j == 1) { 504 yt = t; 505 xt = l; 506 limit = tsize; 507 boundary = user->top; 508 } else if (j == 2) { 509 yt = b; 510 xt = l; 511 limit = lsize; 512 boundary = user->left; 513 } else { /* if (j==3) */ 514 yt = b; 515 xt = r; 516 limit = rsize; 517 boundary = user->right; 518 } 519 520 for (i = 0; i < limit; i++) { 521 u1 = xt; 522 u2 = -yt; 523 for (k = 0; k < maxits; k++) { 524 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 525 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 526 fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2); 527 if (fnorm <= tol) break; 528 njac11 = one + u2 * u2 - u1 * u1; 529 njac12 = two * u1 * u2; 530 njac21 = -two * u1 * u2; 531 njac22 = -one - u1 * u1 + u2 * u2; 532 det = njac11 * njac22 - njac21 * njac12; 533 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 534 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 535 } 536 537 boundary[i] = u1 * u1 - u2 * u2; 538 if (j == 0 || j == 1) { 539 xt = xt + hx; 540 } else { /* if (j==2 || j==3) */ 541 yt = yt + hy; 542 } 543 } 544 } 545 PetscFunctionReturn(PETSC_SUCCESS); 546 } 547 548 /* ------------------------------------------------------------------- */ 549 /* 550 MSA_InitialPoint - Calculates the initial guess in one of three ways. 551 552 Input Parameters: 553 . user - user-defined application context 554 . X - vector for initial guess 555 556 Output Parameters: 557 . X - newly computed initial guess 558 */ 559 static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X) 560 { 561 PetscInt start = -1, i, j; 562 PetscReal zero = 0.0; 563 PetscBool flg; 564 565 PetscFunctionBeginUser; 566 PetscCall(VecSet(X, zero)); 567 PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg)); 568 569 if (flg && start == 0) { /* The zero vector is reasonable */ 570 PetscCall(VecSet(X, zero)); 571 } else { /* Take an average of the boundary conditions */ 572 PetscInt row; 573 PetscInt mx = user->mx, my = user->my; 574 PetscReal *x; 575 576 /* Get pointers to vector data */ 577 PetscCall(VecGetArray(X, &x)); 578 /* Perform local computations */ 579 for (j = 0; j < my; j++) { 580 for (i = 0; i < mx; i++) { 581 row = (j)*mx + (i); 582 x[row] = (((j + 1) * user->bottom[i + 1] + (my - j + 1) * user->top[i + 1]) / (my + 2) + ((i + 1) * user->left[j + 1] + (mx - i + 1) * user->right[j + 1]) / (mx + 2)) / 2.0; 583 } 584 } 585 /* Restore vectors */ 586 PetscCall(VecRestoreArray(X, &x)); 587 } 588 PetscFunctionReturn(PETSC_SUCCESS); 589 } 590 591 /*TEST 592 593 build: 594 requires: !complex 595 596 test: 597 args: -tao_smonitor -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4 598 requires: !single 599 600 test: 601 suffix: 2 602 args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3 603 requires: !single 604 605 test: 606 suffix: 3 607 args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3 608 requires: !single 609 610 test: 611 suffix: 4 612 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 613 requires: !single 614 615 test: 616 suffix: 4_ew 617 args: -tao_smonitor -tao_type bntr -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4 618 requires: !single 619 620 test: 621 suffix: 5 622 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 623 requires: !single 624 625 test: 626 suffix: 5_ew 627 args: -tao_smonitor -tao_type bntl -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4 628 requires: !single 629 630 test: 631 suffix: 6 632 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4 633 requires: !single 634 635 test: 636 suffix: 6_ew 637 args: -tao_smonitor -tao_type bnls -tao_ksp_ew -tao_bnk_ksp_ew_version 3 -mx 10 -my 8 -tao_gatol 1.e-4 638 requires: !single 639 640 test: 641 suffix: 7 642 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 643 requires: !single 644 645 test: 646 suffix: 8 647 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 648 requires: !single 649 650 test: 651 suffix: 9 652 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 653 requires: !single 654 655 test: 656 suffix: 10 657 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian 658 659 test: 660 suffix: 11 661 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian 662 requires: !single 663 664 test: 665 suffix: 12 666 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian 667 requires: !single 668 669 TEST*/ 670