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