1! Program usage: mpiexec -n <proc> eptorsion2f [all TAO options] 2! 3! Description: This example demonstrates use of the TAO package to solve 4! unconstrained minimization problems in parallel. This example is based 5! on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite. 6! The command line options are: 7! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction 8! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction 9! -par <param>, where <param> = angle of twist per unit length 10! 11! 12! ---------------------------------------------------------------------- 13! 14! Elastic-plastic torsion problem. 15! 16! The elastic plastic torsion problem arises from the deconverged 17! of the stress field on an infinitely long cylindrical bar, which is 18! equivalent to the solution of the following problem: 19! min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)} 20! where C is the torsion angle per unit length. 21! 22! The C version of this code is eptorsion2.c 23! 24! ---------------------------------------------------------------------- 25#include "petsc/finclude/petscdmda.h" 26#include "petsc/finclude/petsctao.h" 27module eptorsion2fmodule 28 use petscdmda 29 use petsctao 30 implicit none 31 32 type(tVec) localX 33 type(tDM) dm 34 PetscReal param 35 PetscInt mx, my 36end module 37 38use eptorsion2fmodule 39implicit none 40! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 41! Variable declarations 42! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 43! 44! See additional variable declarations in the file eptorsion2f.h 45! 46PetscErrorCode ierr ! used to check for functions returning nonzeros 47type(tVec) x ! solution vector 48type(tMat) H ! hessian matrix 49PetscInt Nx, Ny ! number of processes in x- and y- directions 50type(tTao) ta ! Tao solver context 51PetscBool flg 52PetscInt i1 53PetscInt dummy 54 55! Note: Any user-defined Fortran routines (such as FormGradient) 56! MUST be declared as external. 57 58external FormInitialGuess, FormFunctionGradient, ComputeHessian 59external Monitor, ConvergenceTest 60 61i1 = 1 62 63! Initialize TAO, PETSc contexts 64PetscCallA(PetscInitialize(ierr)) 65 66! Specify default parameters 67param = 5.0 68mx = 10 69my = 10 70Nx = PETSC_DECIDE 71Ny = PETSC_DECIDE 72 73! Check for any command line arguments that might override defaults 74PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr)) 75PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr)) 76PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', param, flg, ierr)) 77 78! Set up distributed array and vectors 79PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, Nx, Ny, i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr)) 80PetscCallA(DMSetFromOptions(dm, ierr)) 81PetscCallA(DMSetUp(dm, ierr)) 82 83! Create vectors 84PetscCallA(DMCreateGlobalVector(dm, x, ierr)) 85PetscCallA(DMCreateLocalVector(dm, localX, ierr)) 86 87! Create Hessian 88PetscCallA(DMCreateMatrix(dm, H, ierr)) 89PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr)) 90 91! The TAO code begins here 92 93! Create TAO solver 94PetscCallA(TaoCreate(PETSC_COMM_WORLD, ta, ierr)) 95PetscCallA(TaoSetType(ta, TAOCG, ierr)) 96 97! Set routines for function and gradient evaluation 98 99PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr)) 100PetscCallA(TaoSetHessian(ta, H, H, ComputeHessian, 0, ierr)) 101 102! Set initial guess 103PetscCallA(FormInitialGuess(x, ierr)) 104PetscCallA(TaoSetSolution(ta, x, ierr)) 105 106PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testmonitor', flg, ierr)) 107if (flg) then 108 PetscCallA(TaoMonitorSet(ta, Monitor, dummy, PETSC_NULL_FUNCTION, ierr)) 109end if 110 111PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testconvergence', flg, ierr)) 112if (flg) then 113 PetscCallA(TaoSetConvergenceTest(ta, ConvergenceTest, dummy, ierr)) 114end if 115 116! Check for any TAO command line options 117PetscCallA(TaoSetFromOptions(ta, ierr)) 118 119! SOLVE THE APPLICATION 120PetscCallA(TaoSolve(ta, ierr)) 121 122! Free TAO data structures 123PetscCallA(TaoDestroy(ta, ierr)) 124 125! Free PETSc data structures 126PetscCallA(VecDestroy(x, ierr)) 127PetscCallA(VecDestroy(localX, ierr)) 128PetscCallA(MatDestroy(H, ierr)) 129PetscCallA(DMDestroy(dm, ierr)) 130 131! Finalize TAO and PETSc 132PetscCallA(PetscFinalize(ierr)) 133end 134 135! --------------------------------------------------------------------- 136! 137! FormInitialGuess - Computes an initial approximation to the solution. 138! 139! Input Parameters: 140! X - vector 141! 142! Output Parameters: 143! X - vector 144! ierr - error code 145! 146subroutine FormInitialGuess(X, ierr) 147 use eptorsion2fmodule 148 implicit none 149 150! Input/output variables: 151 Vec X 152 PetscErrorCode ierr 153 154! Local variables: 155 PetscInt i, j, k, xe, ye 156 PetscReal temp, val, hx, hy 157 PetscInt xs, ys, xm, ym 158 PetscInt gxm, gym, gxs, gys 159 PetscInt i1 160 161 i1 = 1 162 hx = 1.0/real(mx + 1) 163 hy = 1.0/real(my + 1) 164 165! Get corner information 166 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 167 PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)) 168 169! Compute initial guess over locally owned part of mesh 170 xe = xs + xm 171 ye = ys + ym 172 do j = ys, ye - 1 173 temp = min(j + 1, my - j)*hy 174 do i = xs, xe - 1 175 k = (j - gys)*gxm + i - gxs 176 val = min((min(i + 1, mx - i))*hx, temp) 177 PetscCall(VecSetValuesLocal(X, i1, [k], [val], ADD_VALUES, ierr)) 178 end do 179 end do 180 PetscCall(VecAssemblyBegin(X, ierr)) 181 PetscCall(VecAssemblyEnd(X, ierr)) 182end 183 184! --------------------------------------------------------------------- 185! 186! FormFunctionGradient - Evaluates gradient G(X). 187! 188! Input Parameters: 189! tao - the Tao context 190! X - input vector 191! dummy - optional user-defined context (not used here) 192! 193! Output Parameters: 194! f - the function value at X 195! G - vector containing the newly evaluated gradient 196! ierr - error code 197! 198! Notes: 199! This routine serves as a wrapper for the lower-level routine 200! "ApplicationGradient", where the actual computations are 201! done using the standard Fortran style of treating the local 202! input vector data as an array over the local mesh. 203! 204subroutine FormFunctionGradient(ta, X, f, G, dummy, ierr) 205 use eptorsion2fmodule 206 implicit none 207 208! Input/output variables: 209 type(tTao) ta 210 type(tVec) X, G 211 PetscReal f 212 PetscErrorCode ierr 213 PetscInt dummy 214 215! Declarations for use with local array: 216 217 PetscReal, pointer :: lx_v(:) 218 219! Local variables: 220 PetscReal zero, p5, area, cdiv3 221 PetscReal val, flin, fquad, floc 222 PetscReal v, vb, vl, vr, vt, dvdx 223 PetscReal dvdy, hx, hy 224 PetscInt xe, ye, xsm, ysm 225 PetscInt xep, yep, i, j, k, ind 226 PetscInt xs, ys, xm, ym 227 PetscInt gxs, gys, gxm, gym 228 PetscInt i1 229 230 i1 = 1 231 ierr = 0 232 cdiv3 = param/3.0 233 zero = 0.0 234 p5 = 0.5 235 hx = 1.0/real(mx + 1) 236 hy = 1.0/real(my + 1) 237 fquad = zero 238 flin = zero 239 240! Initialize gradient to zero 241 PetscCall(VecSet(G, zero, ierr)) 242 243! Scatter ghost points to local vector 244 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr)) 245 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr)) 246 247! Get corner information 248 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 249 PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)) 250 251! Get pointer to vector data. 252 PetscCall(VecGetArrayRead(localX, lx_v, ierr)) 253 254! Set local loop dimensions 255 xe = xs + xm 256 ye = ys + ym 257 if (xs == 0) then 258 xsm = xs - 1 259 else 260 xsm = xs 261 end if 262 if (ys == 0) then 263 ysm = ys - 1 264 else 265 ysm = ys 266 end if 267 if (xe == mx) then 268 xep = xe + 1 269 else 270 xep = xe 271 end if 272 if (ye == my) then 273 yep = ye + 1 274 else 275 yep = ye 276 end if 277 278! Compute local gradient contributions over the lower triangular elements 279 280 do j = ysm, ye - 1 281 do i = xsm, xe - 1 282 k = (j - gys)*gxm + i - gxs 283 v = zero 284 vr = zero 285 vt = zero 286 if (i >= 0 .and. j >= 0) v = lx_v(k + 1) 287 if (i < mx - 1 .and. j > -1) vr = lx_v(k + 2) 288 if (i > -1 .and. j < my - 1) vt = lx_v(k + 1 + gxm) 289 dvdx = (vr - v)/hx 290 dvdy = (vt - v)/hy 291 if (i /= -1 .and. j /= -1) then 292 ind = k 293 val = -dvdx/hx - dvdy/hy - cdiv3 294 PetscCall(VecSetValuesLocal(G, i1, [k], [val], ADD_VALUES, ierr)) 295 end if 296 if (i /= mx - 1 .and. j /= -1) then 297 ind = k + 1 298 val = dvdx/hx - cdiv3 299 PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr)) 300 end if 301 if (i /= -1 .and. j /= my - 1) then 302 ind = k + gxm 303 val = dvdy/hy - cdiv3 304 PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr)) 305 end if 306 fquad = fquad + dvdx*dvdx + dvdy*dvdy 307 flin = flin - cdiv3*(v + vr + vt) 308 end do 309 end do 310 311! Compute local gradient contributions over the upper triangular elements 312 313 do j = ys, yep - 1 314 do i = xs, xep - 1 315 k = (j - gys)*gxm + i - gxs 316 vb = zero 317 vl = zero 318 v = zero 319 if (i < mx .and. j > 0) vb = lx_v(k + 1 - gxm) 320 if (i > 0 .and. j < my) vl = lx_v(k) 321 if (i < mx .and. j < my) v = lx_v(1 + k) 322 dvdx = (v - vl)/hx 323 dvdy = (v - vb)/hy 324 if (i /= mx .and. j /= 0) then 325 ind = k - gxm 326 val = -dvdy/hy - cdiv3 327 PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr)) 328 end if 329 if (i /= 0 .and. j /= my) then 330 ind = k - 1 331 val = -dvdx/hx - cdiv3 332 PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr)) 333 end if 334 if (i /= mx .and. j /= my) then 335 ind = k 336 val = dvdx/hx + dvdy/hy - cdiv3 337 PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr)) 338 end if 339 fquad = fquad + dvdx*dvdx + dvdy*dvdy 340 flin = flin - cdiv3*(vb + vl + v) 341 end do 342 end do 343 344! Restore vector 345 PetscCall(VecRestoreArrayRead(localX, lx_v, ierr)) 346 347! Assemble gradient vector 348 PetscCall(VecAssemblyBegin(G, ierr)) 349 PetscCall(VecAssemblyEnd(G, ierr)) 350 351! Scale the gradient 352 area = p5*hx*hy 353 floc = area*(p5*fquad + flin) 354 PetscCall(VecScale(G, area, ierr)) 355 356! Sum function contributions from all processes 357 PetscCallMPI(MPI_Allreduce(floc, f, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD, ierr)) 358 PetscCall(PetscLogFlops(20.0d0*(ye - ysm)*(xe - xsm) + 16.0d0*(xep - xs)*(yep - ys), ierr)) 359end 360 361subroutine ComputeHessian(ta, X, H, Hpre, dummy, ierr) 362 use eptorsion2fmodule 363 implicit none 364 365 type(tTao) ta 366 type(tVec) X 367 type(tMat) H, Hpre 368 PetscErrorCode ierr 369 PetscInt dummy 370 371 PetscInt i, j, k 372 PetscInt col(0:4), row 373 PetscInt xs, xm, gxs, gxm 374 PetscInt ys, ym, gys, gym 375 PetscReal v(0:4) 376 PetscInt i1 377 378 i1 = 1 379 380! Get local grid boundaries 381 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 382 PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)) 383 384 do j = ys, ys + ym - 1 385 do i = xs, xs + xm - 1 386 row = (j - gys)*gxm + (i - gxs) 387 388 k = 0 389 if (j > gys) then 390 v(k) = -1.0 391 col(k) = row - gxm 392 k = k + 1 393 end if 394 395 if (i > gxs) then 396 v(k) = -1.0 397 col(k) = row - 1 398 k = k + 1 399 end if 400 401 v(k) = 4.0 402 col(k) = row 403 k = k + 1 404 405 if (i + 1 < gxs + gxm) then 406 v(k) = -1.0 407 col(k) = row + 1 408 k = k + 1 409 end if 410 411 if (j + 1 < gys + gym) then 412 v(k) = -1.0 413 col(k) = row + gxm 414 k = k + 1 415 end if 416 417 PetscCall(MatSetValuesLocal(H, i1, [row], k, col, v, INSERT_VALUES, ierr)) 418 end do 419 end do 420 421! Assemble matrix 422 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr)) 423 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr)) 424 425! Tell the matrix we will never add a new nonzero location to the 426! matrix. If we do it will generate an error. 427 428 PetscCall(MatSetOption(H, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr)) 429 PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr)) 430 431 PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm, ierr)) 432 433 ierr = 0 434end 435 436subroutine Monitor(ta, dummy, ierr) 437 use eptorsion2fmodule 438 implicit none 439 440 type(tTao) ta 441 PetscInt dummy 442 PetscErrorCode ierr 443 444 PetscInt its 445 PetscReal f, gnorm, cnorm, xdiff 446 TaoConvergedReason reason 447 448 PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr)) 449 if (mod(its, 5) /= 0) then 450 PetscCall(PetscPrintf(PETSC_COMM_WORLD, 'iteration multiple of 5\n', ierr)) 451 end if 452 453 ierr = 0 454 455end 456 457subroutine ConvergenceTest(ta, dummy, ierr) 458 use eptorsion2fmodule 459 implicit none 460 461 type(tTao) ta 462 PetscInt dummy 463 PetscErrorCode ierr 464 465 PetscInt its 466 PetscReal f, gnorm, cnorm, xdiff 467 TaoConvergedReason reason 468 469 PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr)) 470 if (its == 7) then 471 PetscCall(TaoSetConvergedReason(ta, TAO_DIVERGED_MAXITS, ierr)) 472 end if 473 474 ierr = 0 475 476end 477 478!/*TEST 479! 480! build: 481! requires: !complex 482! 483! test: 484! args: -tao_monitor_short -tao_type nls -tao_gttol 1.e-2 485! 486! test: 487! suffix: 2 488! nsize: 2 489! args: -tao_monitor_short -tao_type lmvm -tao_gttol 1.e-2 490! 491! test: 492! suffix: 3 493! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2 494!TEST*/ 495