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