1 #include <petsctao.h> 2 3 static char help[] = "This example demonstrates use of the TAO package to\n\ 4 solve an unconstrained system of equations. This example is based on a\n\ 5 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and\n\ 6 boundary values along the edges of the domain, the objective is to find the\n\ 7 surface with the minimal area that satisfies the boundary conditions.\n\ 8 This application solves this problem using complimentarity -- We are actually\n\ 9 solving the system (grad f)_i >= 0, if x_i == l_i \n\ 10 (grad f)_i = 0, if l_i < x_i < u_i \n\ 11 (grad f)_i <= 0, if x_i == u_i \n\ 12 where f is the function to be minimized. \n\ 13 \n\ 14 The command line options are:\n\ 15 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ 16 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ 17 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n"; 18 19 /* 20 User-defined application context - contains data needed by the 21 application-provided call-back routines, FormFunctionGradient(), 22 FormHessian(). 23 */ 24 typedef struct { 25 PetscInt mx, my; 26 PetscReal *bottom, *top, *left, *right; 27 } AppCtx; 28 29 /* -------- User-defined Routines --------- */ 30 31 static PetscErrorCode MSA_BoundaryConditions(AppCtx *); 32 static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec); 33 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *); 34 PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *); 35 36 int main(int argc, char **argv) { 37 Vec x; /* solution vector */ 38 Vec c; /* Constraints function vector */ 39 Vec xl, xu; /* Bounds on the variables */ 40 PetscBool flg; /* A return variable when checking for user options */ 41 Tao tao; /* TAO solver context */ 42 Mat J; /* Jacobian matrix */ 43 PetscInt N; /* Number of elements in vector */ 44 PetscScalar lb = PETSC_NINFINITY; /* lower bound constant */ 45 PetscScalar ub = PETSC_INFINITY; /* upper bound constant */ 46 AppCtx user; /* user-defined work context */ 47 48 /* Initialize PETSc, TAO */ 49 PetscFunctionBeginUser; 50 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 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 /* Calculate any derived values from parameters */ 61 N = user.mx * user.my; 62 63 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Minimum Surface Area Problem -----\n")); 64 PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx:%" PetscInt_FMT ", my:%" PetscInt_FMT "\n", user.mx, user.my)); 65 66 /* Create appropriate vectors and matrices */ 67 PetscCall(VecCreateSeq(MPI_COMM_SELF, N, &x)); 68 PetscCall(VecDuplicate(x, &c)); 69 PetscCall(MatCreateSeqAIJ(MPI_COMM_SELF, N, N, 7, NULL, &J)); 70 71 /* The TAO code begins here */ 72 73 /* Create TAO solver and set desired solution method */ 74 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao)); 75 PetscCall(TaoSetType(tao, TAOSSILS)); 76 77 /* Set data structure */ 78 PetscCall(TaoSetSolution(tao, x)); 79 80 /* Set routines for constraints function and Jacobian evaluation */ 81 PetscCall(TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user)); 82 PetscCall(TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user)); 83 84 /* Set the variable bounds */ 85 PetscCall(MSA_BoundaryConditions(&user)); 86 87 /* Set initial solution guess */ 88 PetscCall(MSA_InitialPoint(&user, x)); 89 90 /* Set Bounds on variables */ 91 PetscCall(VecDuplicate(x, &xl)); 92 PetscCall(VecDuplicate(x, &xu)); 93 PetscCall(VecSet(xl, lb)); 94 PetscCall(VecSet(xu, ub)); 95 PetscCall(TaoSetVariableBounds(tao, xl, xu)); 96 97 /* Check for any tao command line options */ 98 PetscCall(TaoSetFromOptions(tao)); 99 100 /* Solve the application */ 101 PetscCall(TaoSolve(tao)); 102 103 /* Free Tao data structures */ 104 PetscCall(TaoDestroy(&tao)); 105 106 /* Free PETSc data structures */ 107 PetscCall(VecDestroy(&x)); 108 PetscCall(VecDestroy(&xl)); 109 PetscCall(VecDestroy(&xu)); 110 PetscCall(VecDestroy(&c)); 111 PetscCall(MatDestroy(&J)); 112 113 /* Free user-created data structures */ 114 PetscCall(PetscFree(user.bottom)); 115 PetscCall(PetscFree(user.top)); 116 PetscCall(PetscFree(user.left)); 117 PetscCall(PetscFree(user.right)); 118 119 PetscCall(PetscFinalize()); 120 return 0; 121 } 122 123 /* -------------------------------------------------------------------- */ 124 125 /* FormConstraints - Evaluates gradient of f. 126 127 Input Parameters: 128 . tao - the TAO_APPLICATION context 129 . X - input vector 130 . ptr - optional user-defined context, as set by TaoSetConstraintsRoutine() 131 132 Output Parameters: 133 . G - vector containing the newly evaluated gradient 134 */ 135 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec G, void *ptr) { 136 AppCtx *user = (AppCtx *)ptr; 137 PetscInt i, j, row; 138 PetscInt mx = user->mx, my = user->my; 139 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy; 140 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 141 PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc; 142 PetscScalar zero = 0.0; 143 PetscScalar *g, *x; 144 145 PetscFunctionBegin; 146 /* Initialize vector to zero */ 147 PetscCall(VecSet(G, zero)); 148 149 /* Get pointers to vector data */ 150 PetscCall(VecGetArray(X, &x)); 151 PetscCall(VecGetArray(G, &g)); 152 153 /* Compute function over the locally owned part of the mesh */ 154 for (j = 0; j < my; j++) { 155 for (i = 0; i < mx; i++) { 156 row = j * mx + i; 157 158 xc = x[row]; 159 xlt = xrb = xl = xr = xb = xt = xc; 160 161 if (i == 0) { /* left side */ 162 xl = user->left[j + 1]; 163 xlt = user->left[j + 2]; 164 } else { 165 xl = x[row - 1]; 166 } 167 168 if (j == 0) { /* bottom side */ 169 xb = user->bottom[i + 1]; 170 xrb = user->bottom[i + 2]; 171 } else { 172 xb = x[row - mx]; 173 } 174 175 if (i + 1 == mx) { /* right side */ 176 xr = user->right[j + 1]; 177 xrb = user->right[j]; 178 } else { 179 xr = x[row + 1]; 180 } 181 182 if (j + 1 == 0 + my) { /* top side */ 183 xt = user->top[i + 1]; 184 xlt = user->top[i]; 185 } else { 186 xt = x[row + mx]; 187 } 188 189 if (i > 0 && j + 1 < my) { xlt = x[row - 1 + mx]; } 190 if (j > 0 && i + 1 < mx) { xrb = x[row + 1 - mx]; } 191 192 d1 = (xc - xl); 193 d2 = (xc - xr); 194 d3 = (xc - xt); 195 d4 = (xc - xb); 196 d5 = (xr - xrb); 197 d6 = (xrb - xb); 198 d7 = (xlt - xl); 199 d8 = (xt - xlt); 200 201 df1dxc = d1 * hydhx; 202 df2dxc = (d1 * hydhx + d4 * hxdhy); 203 df3dxc = d3 * hxdhy; 204 df4dxc = (d2 * hydhx + d3 * hxdhy); 205 df5dxc = d2 * hydhx; 206 df6dxc = d4 * hxdhy; 207 208 d1 /= hx; 209 d2 /= hx; 210 d3 /= hy; 211 d4 /= hy; 212 d5 /= hy; 213 d6 /= hx; 214 d7 /= hy; 215 d8 /= hx; 216 217 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 218 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 219 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 220 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 221 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 222 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 223 224 df1dxc /= f1; 225 df2dxc /= f2; 226 df3dxc /= f3; 227 df4dxc /= f4; 228 df5dxc /= f5; 229 df6dxc /= f6; 230 231 g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0; 232 } 233 } 234 235 /* Restore vectors */ 236 PetscCall(VecRestoreArray(X, &x)); 237 PetscCall(VecRestoreArray(G, &g)); 238 PetscCall(PetscLogFlops(67 * mx * my)); 239 PetscFunctionReturn(0); 240 } 241 242 /* ------------------------------------------------------------------- */ 243 /* 244 FormJacobian - Evaluates Jacobian matrix. 245 246 Input Parameters: 247 . tao - the TAO_APPLICATION context 248 . X - input vector 249 . ptr - optional user-defined context, as set by TaoSetJacobian() 250 251 Output Parameters: 252 . tH - Jacobian matrix 253 254 */ 255 PetscErrorCode FormJacobian(Tao tao, Vec X, Mat H, Mat tHPre, void *ptr) { 256 AppCtx *user = (AppCtx *)ptr; 257 PetscInt i, j, k, row; 258 PetscInt mx = user->mx, my = user->my; 259 PetscInt col[7]; 260 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy; 261 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 262 PetscReal hl, hr, ht, hb, hc, htl, hbr; 263 const PetscScalar *x; 264 PetscScalar v[7]; 265 PetscBool assembled; 266 267 /* Set various matrix options */ 268 PetscFunctionBegin; 269 PetscCall(MatSetOption(H, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE)); 270 PetscCall(MatAssembled(H, &assembled)); 271 if (assembled) PetscCall(MatZeroEntries(H)); 272 273 /* Get pointers to vector data */ 274 PetscCall(VecGetArrayRead(X, &x)); 275 276 /* Compute Jacobian over the locally owned part of the mesh */ 277 for (i = 0; i < mx; i++) { 278 for (j = 0; j < my; j++) { 279 row = j * mx + i; 280 281 xc = x[row]; 282 xlt = xrb = xl = xr = xb = xt = xc; 283 284 /* Left side */ 285 if (i == 0) { 286 xl = user->left[j + 1]; 287 xlt = user->left[j + 2]; 288 } else { 289 xl = x[row - 1]; 290 } 291 292 if (j == 0) { 293 xb = user->bottom[i + 1]; 294 xrb = user->bottom[i + 2]; 295 } else { 296 xb = x[row - mx]; 297 } 298 299 if (i + 1 == mx) { 300 xr = user->right[j + 1]; 301 xrb = user->right[j]; 302 } else { 303 xr = x[row + 1]; 304 } 305 306 if (j + 1 == my) { 307 xt = user->top[i + 1]; 308 xlt = user->top[i]; 309 } else { 310 xt = x[row + mx]; 311 } 312 313 if (i > 0 && j + 1 < my) { xlt = x[row - 1 + mx]; } 314 if (j > 0 && i + 1 < mx) { xrb = x[row + 1 - mx]; } 315 316 d1 = (xc - xl) / hx; 317 d2 = (xc - xr) / hx; 318 d3 = (xc - xt) / hy; 319 d4 = (xc - xb) / hy; 320 d5 = (xrb - xr) / hy; 321 d6 = (xrb - xb) / hx; 322 d7 = (xlt - xl) / hy; 323 d8 = (xlt - xt) / hx; 324 325 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 326 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 327 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 328 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 329 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 330 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 331 332 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 333 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 334 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 335 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 336 337 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 338 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 339 340 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); 341 342 hl /= 2.0; 343 hr /= 2.0; 344 ht /= 2.0; 345 hb /= 2.0; 346 hbr /= 2.0; 347 htl /= 2.0; 348 hc /= 2.0; 349 350 k = 0; 351 if (j > 0) { 352 v[k] = hb; 353 col[k] = row - mx; 354 k++; 355 } 356 357 if (j > 0 && i < mx - 1) { 358 v[k] = hbr; 359 col[k] = row - mx + 1; 360 k++; 361 } 362 363 if (i > 0) { 364 v[k] = hl; 365 col[k] = row - 1; 366 k++; 367 } 368 369 v[k] = hc; 370 col[k] = row; 371 k++; 372 373 if (i < mx - 1) { 374 v[k] = hr; 375 col[k] = row + 1; 376 k++; 377 } 378 379 if (i > 0 && j < my - 1) { 380 v[k] = htl; 381 col[k] = row + mx - 1; 382 k++; 383 } 384 385 if (j < my - 1) { 386 v[k] = ht; 387 col[k] = row + mx; 388 k++; 389 } 390 391 /* 392 Set matrix values using local numbering, which was defined 393 earlier, in the main routine. 394 */ 395 PetscCall(MatSetValues(H, 1, &row, k, col, v, INSERT_VALUES)); 396 } 397 } 398 399 /* Restore vectors */ 400 PetscCall(VecRestoreArrayRead(X, &x)); 401 402 /* Assemble the matrix */ 403 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 404 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 405 PetscCall(PetscLogFlops(199 * mx * my)); 406 PetscFunctionReturn(0); 407 } 408 409 /* ------------------------------------------------------------------- */ 410 /* 411 MSA_BoundaryConditions - Calculates the boundary conditions for 412 the region. 413 414 Input Parameter: 415 . user - user-defined application context 416 417 Output Parameter: 418 . user - user-defined application context 419 */ 420 static PetscErrorCode MSA_BoundaryConditions(AppCtx *user) { 421 PetscInt i, j, k, limit = 0, maxits = 5; 422 PetscInt mx = user->mx, my = user->my; 423 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 424 PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10; 425 PetscReal fnorm, det, hx, hy, xt = 0, yt = 0; 426 PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 427 PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5; 428 PetscReal *boundary; 429 430 PetscFunctionBegin; 431 bsize = mx + 2; 432 lsize = my + 2; 433 rsize = my + 2; 434 tsize = mx + 2; 435 436 PetscCall(PetscMalloc1(bsize, &user->bottom)); 437 PetscCall(PetscMalloc1(tsize, &user->top)); 438 PetscCall(PetscMalloc1(lsize, &user->left)); 439 PetscCall(PetscMalloc1(rsize, &user->right)); 440 441 hx = (r - l) / (mx + 1); 442 hy = (t - b) / (my + 1); 443 444 for (j = 0; j < 4; j++) { 445 if (j == 0) { 446 yt = b; 447 xt = l; 448 limit = bsize; 449 boundary = user->bottom; 450 } else if (j == 1) { 451 yt = t; 452 xt = l; 453 limit = tsize; 454 boundary = user->top; 455 } else if (j == 2) { 456 yt = b; 457 xt = l; 458 limit = lsize; 459 boundary = user->left; 460 } else { /* if (j==3) */ 461 yt = b; 462 xt = r; 463 limit = rsize; 464 boundary = user->right; 465 } 466 467 for (i = 0; i < limit; i++) { 468 u1 = xt; 469 u2 = -yt; 470 for (k = 0; k < maxits; k++) { 471 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 472 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 473 fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2); 474 if (fnorm <= tol) break; 475 njac11 = one + u2 * u2 - u1 * u1; 476 njac12 = two * u1 * u2; 477 njac21 = -two * u1 * u2; 478 njac22 = -one - u1 * u1 + u2 * u2; 479 det = njac11 * njac22 - njac21 * njac12; 480 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 481 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 482 } 483 484 boundary[i] = u1 * u1 - u2 * u2; 485 if (j == 0 || j == 1) { 486 xt = xt + hx; 487 } else { /* if (j==2 || j==3) */ 488 yt = yt + hy; 489 } 490 } 491 } 492 PetscFunctionReturn(0); 493 } 494 495 /* ------------------------------------------------------------------- */ 496 /* 497 MSA_InitialPoint - Calculates the initial guess in one of three ways. 498 499 Input Parameters: 500 . user - user-defined application context 501 . X - vector for initial guess 502 503 Output Parameters: 504 . X - newly computed initial guess 505 */ 506 static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X) { 507 PetscInt start = -1, i, j; 508 PetscScalar zero = 0.0; 509 PetscBool flg; 510 511 PetscFunctionBegin; 512 PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg)); 513 514 if (flg && start == 0) { /* The zero vector is reasonable */ 515 PetscCall(VecSet(X, zero)); 516 } else { /* Take an average of the boundary conditions */ 517 PetscInt row; 518 PetscInt mx = user->mx, my = user->my; 519 PetscScalar *x; 520 521 /* Get pointers to vector data */ 522 PetscCall(VecGetArray(X, &x)); 523 524 /* Perform local computations */ 525 for (j = 0; j < my; j++) { 526 for (i = 0; i < mx; i++) { 527 row = (j)*mx + (i); 528 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; 529 } 530 } 531 532 /* Restore vectors */ 533 PetscCall(VecRestoreArray(X, &x)); 534 } 535 PetscFunctionReturn(0); 536 } 537 538 /*TEST 539 540 build: 541 requires: !complex 542 543 test: 544 args: -tao_monitor -tao_view -tao_type ssils -tao_gttol 1.e-5 545 requires: !single 546 547 test: 548 suffix: 2 549 args: -tao_monitor -tao_view -tao_type ssfls -tao_gttol 1.e-5 550 551 TEST*/ 552