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(TaoSetMonitor(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 return 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(tao,X,f,G,dummy,ierr) 206 use eptorsion2fmodule 207 implicit none 208 209! Input/output variables: 210 type(tTao) tao 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(VecGetArrayReadF90(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(VecRestoreArrayReadF90(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 return 361 end 362 363 subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr) 364 use eptorsion2fmodule 365 implicit none 366 367 type(tTao) tao 368 type(tVec) X 369 type(tMat) H,Hpre 370 PetscErrorCode ierr 371 PetscInt dummy 372 373 PetscInt i,j,k 374 PetscInt col(0:4),row 375 PetscInt xs,xm,gxs,gxm 376 PetscInt ys,ym,gys,gym 377 PetscReal v(0:4) 378 PetscInt i1 379 380 i1 = 1 381 382! Get local grid boundaries 383 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 384 PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)) 385 386 do j=ys,ys+ym-1 387 do i=xs,xs+xm-1 388 row = (j-gys)*gxm + (i-gxs) 389 390 k = 0 391 if (j .gt. gys) then 392 v(k) = -1.0 393 col(k) = row-gxm 394 k = k + 1 395 endif 396 397 if (i .gt. gxs) then 398 v(k) = -1.0 399 col(k) = row - 1 400 k = k +1 401 endif 402 403 v(k) = 4.0 404 col(k) = row 405 k = k + 1 406 407 if (i+1 .lt. gxs + gxm) then 408 v(k) = -1.0 409 col(k) = row + 1 410 k = k + 1 411 endif 412 413 if (j+1 .lt. gys + gym) then 414 v(k) = -1.0 415 col(k) = row + gxm 416 k = k + 1 417 endif 418 419 PetscCall(MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)) 420 enddo 421 enddo 422 423! Assemble matrix 424 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)) 425 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)) 426 427! Tell the matrix we will never add a new nonzero location to the 428! matrix. If we do it will generate an error. 429 430 PetscCall(MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 431 PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)) 432 433 PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)) 434 435 ierr = 0 436 return 437 end 438 439 subroutine Monitor(tao, dummy, ierr) 440 use eptorsion2fmodule 441 implicit none 442 443 type(tTao) tao 444 PetscInt dummy 445 PetscErrorCode ierr 446 447 PetscInt its 448 PetscReal f,gnorm,cnorm,xdiff 449 TaoConvergedReason reason 450 451 PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr)) 452 if (mod(its,5) .ne. 0) then 453 PetscCall(PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',ierr)) 454 endif 455 456 ierr = 0 457 458 return 459 end 460 461 subroutine ConvergenceTest(tao, dummy, ierr) 462 use eptorsion2fmodule 463 implicit none 464 465 type(tTao) tao 466 PetscInt dummy 467 PetscErrorCode ierr 468 469 PetscInt its 470 PetscReal f,gnorm,cnorm,xdiff 471 TaoConvergedReason reason 472 473 PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr)) 474 if (its .eq. 7) then 475 PetscCall(TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)) 476 endif 477 478 ierr = 0 479 480 return 481 end 482 483!/*TEST 484! 485! build: 486! requires: !complex 487! 488! test: 489! args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2 490! 491! test: 492! suffix: 2 493! nsize: 2 494! args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2 495! 496! test: 497! suffix: 3 498! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2 499!TEST*/ 500