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