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 26 module eptorsion2fmodule 27#include "petsc/finclude/petsctao.h" 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 end module 37 38 use eptorsion2fmodule 39 implicit none 40! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 41! Variable declarations 42! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 43! 44! See additional variable declarations in the file eptorsion2f.h 45! 46 PetscErrorCode ierr ! used to check for functions returning nonzeros 47 type(tVec) x ! solution vector 48 type(tMat) H ! hessian matrix 49 PetscInt Nx, Ny ! number of processes in x- and y- directions 50 type(tTao) tao ! Tao solver context 51 PetscBool flg 52 PetscInt i1 53 PetscInt dummy 54 55! Note: Any user-defined Fortran routines (such as FormGradient) 56! MUST be declared as external. 57 58 external FormInitialGuess,FormFunctionGradient,ComputeHessian 59 external Monitor,ConvergenceTest 60 61 i1 = 1 62 63! Initialize TAO, PETSc contexts 64 PetscCallA(PetscInitialize(ierr)) 65 66! Specify default parameters 67 param = 5.0 68 mx = 10 69 my = 10 70 Nx = PETSC_DECIDE 71 Ny = PETSC_DECIDE 72 73! Check for any command line arguments that might override defaults 74 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)) 75 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)) 76 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',param,flg,ierr)) 77 78! Set up distributed array and vectors 79 PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,dm,ierr)) 80 PetscCallA(DMSetFromOptions(dm,ierr)) 81 PetscCallA(DMSetUp(dm,ierr)) 82 83! Create vectors 84 PetscCallA(DMCreateGlobalVector(dm,x,ierr)) 85 PetscCallA(DMCreateLocalVector(dm,localX,ierr)) 86 87! Create Hessian 88 PetscCallA(DMCreateMatrix(dm,H,ierr)) 89 PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)) 90 91! The TAO code begins here 92 93! Create TAO solver 94 PetscCallA(TaoCreate(PETSC_COMM_WORLD,tao,ierr)) 95 PetscCallA(TaoSetType(tao,TAOCG,ierr)) 96 97! Set routines for function and gradient evaluation 98 99 PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr)) 100 PetscCallA(TaoSetHessian(tao,H,H,ComputeHessian,0,ierr)) 101 102! Set initial guess 103 PetscCallA(FormInitialGuess(x,ierr)) 104 PetscCallA(TaoSetSolution(tao,x,ierr)) 105 106 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testmonitor',flg,ierr)) 107 if (flg) then 108 PetscCallA(TaoMonitorSet(tao,Monitor,dummy,PETSC_NULL_FUNCTION,ierr)) 109 endif 110 111 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testconvergence',flg, ierr)) 112 if (flg) then 113 PetscCallA(TaoSetConvergenceTest(tao,ConvergenceTest,dummy,ierr)) 114 endif 115 116! Check for any TAO command line options 117 PetscCallA(TaoSetFromOptions(tao,ierr)) 118 119! SOLVE THE APPLICATION 120 PetscCallA(TaoSolve(tao,ierr)) 121 122! Free TAO data structures 123 PetscCallA(TaoDestroy(tao,ierr)) 124 125! Free PETSc data structures 126 PetscCallA(VecDestroy(x,ierr)) 127 PetscCallA(VecDestroy(localX,ierr)) 128 PetscCallA(MatDestroy(H,ierr)) 129 PetscCallA(DMDestroy(dm,ierr)) 130 131! Finalize TAO and PETSc 132 PetscCallA(PetscFinalize(ierr)) 133 end 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! 146 subroutine 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)) 182 end 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! 204 subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr) 205 use eptorsion2fmodule 206 implicit none 207 208! Input/output variables: 209 type(tTao) tao 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(VecGetArrayReadF90(localX,lx_v,ierr)) 253 254! Set local loop dimensions 255 xe = xs+xm 256 ye = ys+ym 257 if (xs .eq. 0) then 258 xsm = xs-1 259 else 260 xsm = xs 261 endif 262 if (ys .eq. 0) then 263 ysm = ys-1 264 else 265 ysm = ys 266 endif 267 if (xe .eq. mx) then 268 xep = xe+1 269 else 270 xep = xe 271 endif 272 if (ye .eq. my) then 273 yep = ye+1 274 else 275 yep = ye 276 endif 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 .ge. 0 .and. j .ge. 0) v = lx_v(k+1) 287 if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(k+2) 288 if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(k+1+gxm) 289 dvdx = (vr-v)/hx 290 dvdy = (vt-v)/hy 291 if (i .ne. -1 .and. j .ne. -1) then 292 ind = k 293 val = - dvdx/hx - dvdy/hy - cdiv3 294 PetscCall(VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)) 295 endif 296 if (i .ne. mx-1 .and. j .ne. -1) then 297 ind = k+1 298 val = dvdx/hx - cdiv3 299 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 300 endif 301 if (i .ne. -1 .and. j .ne. my-1) then 302 ind = k+gxm 303 val = dvdy/hy - cdiv3 304 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 305 endif 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 .lt. mx .and. j .gt. 0) vb = lx_v(k+1-gxm) 320 if (i .gt. 0 .and. j .lt. my) vl = lx_v(k) 321 if (i .lt. mx .and. j .lt. my) v = lx_v(1+k) 322 dvdx = (v-vl)/hx 323 dvdy = (v-vb)/hy 324 if (i .ne. mx .and. j .ne. 0) then 325 ind = k-gxm 326 val = - dvdy/hy - cdiv3 327 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 328 endif 329 if (i .ne. 0 .and. j .ne. my) then 330 ind = k-1 331 val = - dvdx/hx - cdiv3 332 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 333 endif 334 if (i .ne. mx .and. j .ne. my) then 335 ind = k 336 val = dvdx/hx + dvdy/hy - cdiv3 337 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 338 endif 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(VecRestoreArrayReadF90(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)) 359 end 360 361 subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr) 362 use eptorsion2fmodule 363 implicit none 364 365 type(tTao) tao 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 .gt. gys) then 390 v(k) = -1.0 391 col(k) = row-gxm 392 k = k + 1 393 endif 394 395 if (i .gt. gxs) then 396 v(k) = -1.0 397 col(k) = row - 1 398 k = k +1 399 endif 400 401 v(k) = 4.0 402 col(k) = row 403 k = k + 1 404 405 if (i+1 .lt. gxs + gxm) then 406 v(k) = -1.0 407 col(k) = row + 1 408 k = k + 1 409 endif 410 411 if (j+1 .lt. gys + gym) then 412 v(k) = -1.0 413 col(k) = row + gxm 414 k = k + 1 415 endif 416 417 PetscCall(MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)) 418 enddo 419 enddo 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 434 end 435 436 subroutine Monitor(tao, dummy, ierr) 437 use eptorsion2fmodule 438 implicit none 439 440 type(tTao) tao 441 PetscInt dummy 442 PetscErrorCode ierr 443 444 PetscInt its 445 PetscReal f,gnorm,cnorm,xdiff 446 TaoConvergedReason reason 447 448 PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr)) 449 if (mod(its,5) .ne. 0) then 450 PetscCall(PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',ierr)) 451 endif 452 453 ierr = 0 454 455 end 456 457 subroutine ConvergenceTest(tao, dummy, ierr) 458 use eptorsion2fmodule 459 implicit none 460 461 type(tTao) tao 462 PetscInt dummy 463 PetscErrorCode ierr 464 465 PetscInt its 466 PetscReal f,gnorm,cnorm,xdiff 467 TaoConvergedReason reason 468 469 PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr)) 470 if (its .eq. 7) then 471 PetscCall(TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)) 472 endif 473 474 ierr = 0 475 476 end 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