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, NULL, 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_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 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 matrix used to compute the preconditioner 280 281 */ 282 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) 283 { 284 AppCtx *user = (AppCtx *)ptr; 285 286 PetscFunctionBeginUser; 287 /* Evaluate the Hessian entries*/ 288 PetscCall(QuadraticH(user, X, H)); 289 PetscFunctionReturn(PETSC_SUCCESS); 290 } 291 292 /* ------------------------------------------------------------------- */ 293 /* 294 QuadraticH - Evaluates the Hessian matrix. 295 296 Input Parameters: 297 . user - user-defined context, as set by TaoSetHessian() 298 . X - input vector 299 300 Output Parameter: 301 . H - Hessian matrix 302 */ 303 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) 304 { 305 PetscInt i, j, k, row; 306 PetscInt mx = user->mx, my = user->my; 307 PetscInt col[7]; 308 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy; 309 PetscReal rhx = mx + 1, rhy = my + 1; 310 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 311 PetscReal hl, hr, ht, hb, hc, htl, hbr; 312 const PetscScalar *x; 313 PetscReal v[7]; 314 315 PetscFunctionBeginUser; 316 /* Get pointers to vector data */ 317 PetscCall(VecGetArrayRead(X, &x)); 318 319 /* Initialize matrix entries to zero */ 320 PetscCall(MatZeroEntries(Hessian)); 321 322 /* Set various matrix options */ 323 PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE)); 324 325 /* Compute Hessian over the locally owned part of the mesh */ 326 for (i = 0; i < mx; i++) { 327 for (j = 0; j < my; j++) { 328 row = (j)*mx + (i); 329 330 xc = x[row]; 331 xlt = xrb = xl = xr = xb = xt = xc; 332 333 /* Left side */ 334 if (i == 0) { 335 xl = user->left[j + 1]; 336 xlt = user->left[j + 2]; 337 } else { 338 xl = x[row - 1]; 339 } 340 341 if (j == 0) { 342 xb = user->bottom[i + 1]; 343 xrb = user->bottom[i + 2]; 344 } else { 345 xb = x[row - mx]; 346 } 347 348 if (i + 1 == mx) { 349 xr = user->right[j + 1]; 350 xrb = user->right[j]; 351 } else { 352 xr = x[row + 1]; 353 } 354 355 if (j + 1 == my) { 356 xt = user->top[i + 1]; 357 xlt = user->top[i]; 358 } else { 359 xt = x[row + mx]; 360 } 361 362 if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx]; 363 if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx]; 364 365 d1 = (xc - xl) * rhx; 366 d2 = (xc - xr) * rhx; 367 d3 = (xc - xt) * rhy; 368 d4 = (xc - xb) * rhy; 369 d5 = (xrb - xr) * rhy; 370 d6 = (xrb - xb) * rhx; 371 d7 = (xlt - xl) * rhy; 372 d8 = (xlt - xt) * rhx; 373 374 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 375 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 376 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 377 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 378 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 379 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 380 381 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 382 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 383 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 384 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 385 386 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 387 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 388 389 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); 390 391 hl *= 0.5; 392 hr *= 0.5; 393 ht *= 0.5; 394 hb *= 0.5; 395 hbr *= 0.5; 396 htl *= 0.5; 397 hc *= 0.5; 398 399 k = 0; 400 if (j > 0) { 401 v[k] = hb; 402 col[k] = row - mx; 403 k++; 404 } 405 406 if (j > 0 && i < mx - 1) { 407 v[k] = hbr; 408 col[k] = row - mx + 1; 409 k++; 410 } 411 412 if (i > 0) { 413 v[k] = hl; 414 col[k] = row - 1; 415 k++; 416 } 417 418 v[k] = hc; 419 col[k] = row; 420 k++; 421 422 if (i < mx - 1) { 423 v[k] = hr; 424 col[k] = row + 1; 425 k++; 426 } 427 428 if (i > 0 && j < my - 1) { 429 v[k] = htl; 430 col[k] = row + mx - 1; 431 k++; 432 } 433 434 if (j < my - 1) { 435 v[k] = ht; 436 col[k] = row + mx; 437 k++; 438 } 439 440 /* 441 Set matrix values using local numbering, which was defined 442 earlier, in the main routine. 443 */ 444 PetscCall(MatSetValues(Hessian, 1, &row, k, col, v, INSERT_VALUES)); 445 } 446 } 447 448 /* Restore vectors */ 449 PetscCall(VecRestoreArrayRead(X, &x)); 450 451 /* Assemble the matrix */ 452 PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY)); 453 PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY)); 454 455 PetscCall(PetscLogFlops(199.0 * mx * my)); 456 PetscFunctionReturn(PETSC_SUCCESS); 457 } 458 459 /* ------------------------------------------------------------------- */ 460 /* 461 MSA_BoundaryConditions - Calculates the boundary conditions for 462 the region. 463 464 Input Parameter: 465 . user - user-defined application context 466 467 Output Parameter: 468 . user - user-defined application context 469 */ 470 static PetscErrorCode MSA_BoundaryConditions(AppCtx *user) 471 { 472 PetscInt i, j, k, limit = 0; 473 PetscInt maxits = 5; 474 PetscInt mx = user->mx, my = user->my; 475 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 476 PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10; 477 PetscReal fnorm, det, hx, hy, xt = 0, yt = 0; 478 PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 479 PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5; 480 PetscReal *boundary; 481 482 PetscFunctionBeginUser; 483 bsize = mx + 2; 484 lsize = my + 2; 485 rsize = my + 2; 486 tsize = mx + 2; 487 488 PetscCall(PetscMalloc1(bsize, &user->bottom)); 489 PetscCall(PetscMalloc1(tsize, &user->top)); 490 PetscCall(PetscMalloc1(lsize, &user->left)); 491 PetscCall(PetscMalloc1(rsize, &user->right)); 492 493 hx = (r - l) / (mx + 1); 494 hy = (t - b) / (my + 1); 495 496 for (j = 0; j < 4; j++) { 497 if (j == 0) { 498 yt = b; 499 xt = l; 500 limit = bsize; 501 boundary = user->bottom; 502 } else if (j == 1) { 503 yt = t; 504 xt = l; 505 limit = tsize; 506 boundary = user->top; 507 } else if (j == 2) { 508 yt = b; 509 xt = l; 510 limit = lsize; 511 boundary = user->left; 512 } else { /* if (j==3) */ 513 yt = b; 514 xt = r; 515 limit = rsize; 516 boundary = user->right; 517 } 518 519 for (i = 0; i < limit; i++) { 520 u1 = xt; 521 u2 = -yt; 522 for (k = 0; k < maxits; k++) { 523 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 524 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 525 fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2); 526 if (fnorm <= tol) break; 527 njac11 = one + u2 * u2 - u1 * u1; 528 njac12 = two * u1 * u2; 529 njac21 = -two * u1 * u2; 530 njac22 = -one - u1 * u1 + u2 * u2; 531 det = njac11 * njac22 - njac21 * njac12; 532 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 533 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 534 } 535 536 boundary[i] = u1 * u1 - u2 * u2; 537 if (j == 0 || j == 1) { 538 xt = xt + hx; 539 } else { /* if (j==2 || j==3) */ 540 yt = yt + hy; 541 } 542 } 543 } 544 PetscFunctionReturn(PETSC_SUCCESS); 545 } 546 547 /* ------------------------------------------------------------------- */ 548 /* 549 MSA_InitialPoint - Calculates the initial guess in one of three ways. 550 551 Input Parameters: 552 . user - user-defined application context 553 . X - vector for initial guess 554 555 Output Parameters: 556 . X - newly computed initial guess 557 */ 558 static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X) 559 { 560 PetscInt start = -1, i, j; 561 PetscReal zero = 0.0; 562 PetscBool flg; 563 564 PetscFunctionBeginUser; 565 PetscCall(VecSet(X, zero)); 566 PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg)); 567 568 if (flg && start == 0) { /* The zero vector is reasonable */ 569 PetscCall(VecSet(X, zero)); 570 } else { /* Take an average of the boundary conditions */ 571 PetscInt row; 572 PetscInt mx = user->mx, my = user->my; 573 PetscReal *x; 574 575 /* Get pointers to vector data */ 576 PetscCall(VecGetArray(X, &x)); 577 /* Perform local computations */ 578 for (j = 0; j < my; j++) { 579 for (i = 0; i < mx; i++) { 580 row = (j)*mx + (i); 581 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; 582 } 583 } 584 /* Restore vectors */ 585 PetscCall(VecRestoreArray(X, &x)); 586 } 587 PetscFunctionReturn(PETSC_SUCCESS); 588 } 589 590 /*TEST 591 592 build: 593 requires: !complex 594 595 test: 596 args: -tao_monitor_short -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4 597 requires: !single 598 599 test: 600 suffix: 2 601 args: -tao_monitor_short -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3 602 requires: !single 603 604 test: 605 suffix: 3 606 args: -tao_monitor_short -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3 607 requires: !single 608 609 test: 610 suffix: 4 611 args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 612 requires: !single 613 614 test: 615 suffix: 4_ew 616 args: -tao_monitor_short -tao_type bntr -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4 617 requires: !single 618 619 test: 620 suffix: 5 621 args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 622 requires: !single 623 624 test: 625 suffix: 5_ew 626 args: -tao_monitor_short -tao_type bntl -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4 627 requires: !single 628 629 test: 630 suffix: 6 631 args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4 632 requires: !single 633 634 test: 635 suffix: 6_ew 636 args: -tao_monitor_short -tao_type bnls -tao_ksp_ew -tao_bnk_ksp_ew_version 3 -mx 10 -my 8 -tao_gatol 1.e-4 637 requires: !single 638 639 test: 640 suffix: 7 641 args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short 642 requires: !single 643 644 test: 645 suffix: 8 646 args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short 647 requires: !single 648 649 test: 650 suffix: 9 651 args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short 652 requires: !single 653 654 test: 655 suffix: 10 656 args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian -tao_bnk_cg_tao_monitor_short 657 658 test: 659 suffix: 11 660 args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian 661 requires: !single 662 663 test: 664 suffix: 12 665 args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian 666 requires: !single 667 668 TEST*/ 669