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 Vec localX 33 DM 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 Vec x ! solution vector 48 Mat H ! hessian matrix 49 PetscInt Nx, Ny ! number of processes in x- and y- directions 50 Tao 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 Tao tao 211 Vec X, G 212 PetscReal f 213 PetscErrorCode ierr 214 PetscInt dummy 215 216! Declarations for use with local array: 217 218! PETSc's VecGetArray acts differently in Fortran than it does in C. 219! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 220! will return an array of doubles referenced by x_array offset by x_index. 221! i.e., to reference the kth element of X, use x_array(k + x_index). 222! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 223 PetscReal lx_v(0:1) 224 PetscOffset lx_i 225 226! Local variables: 227 PetscReal zero, p5, area, cdiv3 228 PetscReal val, flin, fquad,floc 229 PetscReal v, vb, vl, vr, vt, dvdx 230 PetscReal dvdy, hx, hy 231 PetscInt xe, ye, xsm, ysm 232 PetscInt xep, yep, i, j, k, ind 233 PetscInt xs, ys, xm, ym 234 PetscInt gxs, gys, gxm, gym 235 PetscInt i1 236 237 i1 = 1 238 ierr = 0 239 cdiv3 = param/3.0 240 zero = 0.0 241 p5 = 0.5 242 hx = 1.0/real(mx + 1) 243 hy = 1.0/real(my + 1) 244 fquad = zero 245 flin = zero 246 247! Initialize gradient to zero 248 PetscCall(VecSet(G,zero,ierr)) 249 250! Scatter ghost points to local vector 251 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)) 252 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)) 253 254! Get corner information 255 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 256 PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)) 257 258! Get pointer to vector data. 259 PetscCall(VecGetArray(localX,lx_v,lx_i,ierr)) 260 261! Set local loop dimensions 262 xe = xs+xm 263 ye = ys+ym 264 if (xs .eq. 0) then 265 xsm = xs-1 266 else 267 xsm = xs 268 endif 269 if (ys .eq. 0) then 270 ysm = ys-1 271 else 272 ysm = ys 273 endif 274 if (xe .eq. mx) then 275 xep = xe+1 276 else 277 xep = xe 278 endif 279 if (ye .eq. my) then 280 yep = ye+1 281 else 282 yep = ye 283 endif 284 285! Compute local gradient contributions over the lower triangular elements 286 287 do j = ysm, ye-1 288 do i = xsm, xe-1 289 k = (j-gys)*gxm + i-gxs 290 v = zero 291 vr = zero 292 vt = zero 293 if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k) 294 if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1) 295 if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm) 296 dvdx = (vr-v)/hx 297 dvdy = (vt-v)/hy 298 if (i .ne. -1 .and. j .ne. -1) then 299 ind = k 300 val = - dvdx/hx - dvdy/hy - cdiv3 301 PetscCall(VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)) 302 endif 303 if (i .ne. mx-1 .and. j .ne. -1) then 304 ind = k+1 305 val = dvdx/hx - cdiv3 306 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 307 endif 308 if (i .ne. -1 .and. j .ne. my-1) then 309 ind = k+gxm 310 val = dvdy/hy - cdiv3 311 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 312 endif 313 fquad = fquad + dvdx*dvdx + dvdy*dvdy 314 flin = flin - cdiv3 * (v+vr+vt) 315 end do 316 end do 317 318! Compute local gradient contributions over the upper triangular elements 319 320 do j = ys, yep-1 321 do i = xs, xep-1 322 k = (j-gys)*gxm + i-gxs 323 vb = zero 324 vl = zero 325 v = zero 326 if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm) 327 if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1) 328 if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k) 329 dvdx = (v-vl)/hx 330 dvdy = (v-vb)/hy 331 if (i .ne. mx .and. j .ne. 0) then 332 ind = k-gxm 333 val = - dvdy/hy - cdiv3 334 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 335 endif 336 if (i .ne. 0 .and. j .ne. my) then 337 ind = k-1 338 val = - dvdx/hx - cdiv3 339 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 340 endif 341 if (i .ne. mx .and. j .ne. my) then 342 ind = k 343 val = dvdx/hx + dvdy/hy - cdiv3 344 PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)) 345 endif 346 fquad = fquad + dvdx*dvdx + dvdy*dvdy 347 flin = flin - cdiv3*(vb + vl + v) 348 end do 349 end do 350 351! Restore vector 352 PetscCall(VecRestoreArray(localX,lx_v,lx_i,ierr)) 353 354! Assemble gradient vector 355 PetscCall(VecAssemblyBegin(G,ierr)) 356 PetscCall(VecAssemblyEnd(G,ierr)) 357 358! Scale the gradient 359 area = p5*hx*hy 360 floc = area *(p5*fquad+flin) 361 PetscCall(VecScale(G,area,ierr)) 362 363! Sum function contributions from all processes 364 PetscCallMPI(MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr)) 365 PetscCall(PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+16.0d0*(xep-xs)*(yep-ys),ierr)) 366 return 367 end 368 369 subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr) 370 use eptorsion2fmodule 371 implicit none 372 373 Tao tao 374 Vec X 375 Mat H,Hpre 376 PetscErrorCode ierr 377 PetscInt dummy 378 379 PetscInt i,j,k 380 PetscInt col(0:4),row 381 PetscInt xs,xm,gxs,gxm 382 PetscInt ys,ym,gys,gym 383 PetscReal v(0:4) 384 PetscInt i1 385 386 i1 = 1 387 388! Get local grid boundaries 389 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 390 PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)) 391 392 do j=ys,ys+ym-1 393 do i=xs,xs+xm-1 394 row = (j-gys)*gxm + (i-gxs) 395 396 k = 0 397 if (j .gt. gys) then 398 v(k) = -1.0 399 col(k) = row-gxm 400 k = k + 1 401 endif 402 403 if (i .gt. gxs) then 404 v(k) = -1.0 405 col(k) = row - 1 406 k = k +1 407 endif 408 409 v(k) = 4.0 410 col(k) = row 411 k = k + 1 412 413 if (i+1 .lt. gxs + gxm) then 414 v(k) = -1.0 415 col(k) = row + 1 416 k = k + 1 417 endif 418 419 if (j+1 .lt. gys + gym) then 420 v(k) = -1.0 421 col(k) = row + gxm 422 k = k + 1 423 endif 424 425 PetscCall(MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)) 426 enddo 427 enddo 428 429! Assemble matrix 430 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)) 431 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)) 432 433! Tell the matrix we will never add a new nonzero location to the 434! matrix. If we do it will generate an error. 435 436 PetscCall(MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 437 PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)) 438 439 PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)) 440 441 ierr = 0 442 return 443 end 444 445 subroutine Monitor(tao, dummy, ierr) 446 use eptorsion2fmodule 447 implicit none 448 449 Tao tao 450 PetscInt dummy 451 PetscErrorCode ierr 452 453 PetscInt its 454 PetscReal f,gnorm,cnorm,xdiff 455 TaoConvergedReason reason 456 457 PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr)) 458 if (mod(its,5) .ne. 0) then 459 PetscCall(PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',ierr)) 460 endif 461 462 ierr = 0 463 464 return 465 end 466 467 subroutine ConvergenceTest(tao, dummy, ierr) 468 use eptorsion2fmodule 469 implicit none 470 471 Tao tao 472 PetscInt dummy 473 PetscErrorCode ierr 474 475 PetscInt its 476 PetscReal f,gnorm,cnorm,xdiff 477 TaoConvergedReason reason 478 479 PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr)) 480 if (its .eq. 7) then 481 PetscCall(TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)) 482 endif 483 484 ierr = 0 485 486 return 487 end 488 489!/*TEST 490! 491! build: 492! requires: !complex 493! 494! test: 495! args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2 496! 497! test: 498! suffix: 2 499! nsize: 2 500! args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2 501! 502! test: 503! suffix: 3 504! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2 505!TEST*/ 506