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