Actual source code: eptorsion2f.F

      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: !/*T
     12: !   Concepts: TAO^Solving an unconstrained minimization problem
     13: !   Routines: TaoCreate(); TaoSetType();
     14: !   Routines: TaoSetInitialVector();
     15: !   Routines: TaoSetObjectiveAndGradientRoutine();
     16: !   Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
     17: !   Routines: TaoSetMonitor(); TaoSetConvergenceTest()
     18: !   Routines: TaoSolve(); TaoGetSolutionStatus()
     19: !   Routines: TaoDestroy();

     21: !   Processors: n
     22: !T*/
     23: !
     24: ! ----------------------------------------------------------------------
     25: !
     26: !  Elastic-plastic torsion problem.
     27: !
     28: !  The elastic plastic torsion problem arises from the deconverged
     29: !  of the stress field on an infinitely long cylindrical bar, which is
     30: !  equivalent to the solution of the following problem:
     31: !     min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
     32: !  where C is the torsion angle per unit length.
     33: !
     34: !  The C version of this code is eptorsion2.c
     35: !
     36: ! ----------------------------------------------------------------------

     38:       module mymodule
     39: #include "petsc/finclude/petsctao.h"
     40:       use petscdmda
     41:       use petsctao
     42:       implicit none

     44:       Vec              localX
     45:       DM               dm
     46:       PetscReal      param
     47:       PetscInt         mx, my
     48:       end module

     50:       use mymodule
     51:       implicit none
     52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     53: !                   Variable declarations
     54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     55: !
     56: !  See additional variable declarations in the file eptorsion2f.h
     57: !
     58:       PetscErrorCode   ierr           ! used to check for functions returning nonzeros
     59:       Vec              x              ! solution vector
     60:       Mat              H              ! hessian matrix
     61:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
     62:       Tao        tao            ! Tao solver context
     63:       PetscBool        flg
     64:       PetscInt         i1
     65:       PetscInt         dummy

     68: !  Note: Any user-defined Fortran routines (such as FormGradient)
     69: !  MUST be declared as external.

     71:       external FormInitialGuess,FormFunctionGradient,ComputeHessian
     72:       external Monitor,ConvergenceTest

     74:       i1 = 1

     76: !     Initialize TAO, PETSc  contexts
     77:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     78:       if (ierr .ne. 0) then
     79:          print*,'Unable to initialize PETSc'
     80:          stop
     81:       endif

     83: !     Specify default parameters
     84:       param = 5.0
     85:       mx = 10
     86:       my = 10
     87:       Nx = PETSC_DECIDE
     88:       Ny = PETSC_DECIDE

     90: !     Check for any command line arguments that might override defaults
     91:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,    &
     92:      &                        '-mx',mx,flg,ierr)
     93:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,    &
     94:      &                        '-my',my,flg,ierr)
     95:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
     96:      &                         '-par',param,flg,ierr)

     99: !     Set up distributed array and vectors
    100:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,               &
    101:      &     DM_BOUNDARY_NONE,                                             &
    102:      &     DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,        &
    103:      &     PETSC_NULL_INTEGER,dm,ierr)
    104:       call DMSetFromOptions(dm,ierr)
    105:       call DMSetUp(dm,ierr)

    107: !     Create vectors
    108:       call DMCreateGlobalVector(dm,x,ierr)
    109:       call DMCreateLocalVector(dm,localX,ierr)

    111: !     Create Hessian
    112:       call DMCreateMatrix(dm,H,ierr)
    113:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

    115: !     The TAO code begins here

    117: !     Create TAO solver
    118:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
    119:       call TaoSetType(tao,TAOCG,ierr)

    121: !     Set routines for function and gradient evaluation

    123:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
    124:      &     FormFunctionGradient,0,ierr)
    125:       call TaoSetHessianRoutine(tao,H,H,ComputeHessian,                 &
    126:      &     0,ierr)

    128: !     Set initial guess
    129:       call FormInitialGuess(x,ierr)
    130:       call TaoSetInitialVector(tao,x,ierr)

    132:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
    133:      &                         '-testmonitor',flg,ierr)
    134:       if (flg) then
    135:          call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION,      &
    136:      &        ierr)
    137:       endif

    139:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
    140:      &                         '-testconvergence',flg, ierr)
    141:       if (flg) then
    142:          call TaoSetConvergenceTest(tao,ConvergenceTest,dummy,          &
    143:      &        ierr)
    144:       endif

    146: !     Check for any TAO command line options
    147:       call TaoSetFromOptions(tao,ierr)

    150: !     SOLVE THE APPLICATION
    151:       call TaoSolve(tao,ierr)

    153: !     Free TAO data structures
    154:       call TaoDestroy(tao,ierr)

    157: !     Free PETSc data structures
    158:       call VecDestroy(x,ierr)
    159:       call VecDestroy(localX,ierr)
    160:       call MatDestroy(H,ierr)
    161:       call DMDestroy(dm,ierr)

    164: !     Finalize TAO and PETSc
    165:       call PetscFinalize(ierr)
    166:       end

    169: ! ---------------------------------------------------------------------
    170: !
    171: !   FormInitialGuess - Computes an initial approximation to the solution.
    172: !
    173: !   Input Parameters:
    174: !   X    - vector
    175: !
    176: !   Output Parameters:
    177: !   X    - vector
    178: !   ierr - error code
    179: !
    180:       subroutine FormInitialGuess(X,ierr)
    181:       use mymodule
    182:       implicit none

    184: !  Input/output variables:
    185:       Vec              X
    186:       PetscErrorCode   ierr

    188: !  Local variables:
    189:       PetscInt         i, j, k, xe, ye
    190:       PetscReal      temp, val, hx, hy
    191:       PetscInt         xs, ys, xm, ym
    192:       PetscInt         gxm, gym, gxs, gys
    193:       PetscInt         i1

    195:       i1 = 1
    196:       hx = 1.0/real(mx + 1)
    197:       hy = 1.0/real(my + 1)

    199: !  Get corner information
    200:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    201:      &                  PETSC_NULL_INTEGER,ierr)
    202:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
    203:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)

    207: !  Compute initial guess over locally owned part of mesh
    208:       xe = xs+xm
    209:       ye = ys+ym
    210:       do j=ys,ye-1
    211:          temp = min(j+1,my-j)*hy
    212:          do i=xs,xe-1
    213:             k   = (j-gys)*gxm + i-gxs
    214:             val = min((min(i+1,mx-i))*hx,temp)
    215:             call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
    216:          end do
    217:       end do
    218:       call VecAssemblyBegin(X,ierr)
    219:       call VecAssemblyEnd(X,ierr)
    220:       return
    221:       end

    224: ! ---------------------------------------------------------------------
    225: !
    226: !  FormFunctionGradient - Evaluates gradient G(X).
    227: !
    228: !  Input Parameters:
    229: !  tao   - the Tao context
    230: !  X     - input vector
    231: !  dummy - optional user-defined context (not used here)
    232: !
    233: !  Output Parameters:
    234: !  f     - the function value at X
    235: !  G     - vector containing the newly evaluated gradient
    236: !  ierr  - error code
    237: !
    238: !  Notes:
    239: !  This routine serves as a wrapper for the lower-level routine
    240: !  "ApplicationGradient", where the actual computations are
    241: !  done using the standard Fortran style of treating the local
    242: !  input vector data as an array over the local mesh.
    243: !
    244:       subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
    245:       use mymodule
    246:       implicit none

    248: !  Input/output variables:
    249:       Tao        tao
    250:       Vec              X, G
    251:       PetscReal      f
    252:       PetscErrorCode   ierr
    253:       PetscInt         dummy

    255: !  Declarations for use with local array:

    258: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    259: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
    260: ! will return an array of doubles referenced by x_array offset by x_index.
    261: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    262: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    263:       PetscReal      lx_v(0:1)
    264:       PetscOffset      lx_i

    266: !  Local variables:
    267:       PetscReal      zero, p5, area, cdiv3
    268:       PetscReal      val, flin, fquad,floc
    269:       PetscReal      v, vb, vl, vr, vt, dvdx
    270:       PetscReal      dvdy, hx, hy
    271:       PetscInt         xe, ye, xsm, ysm
    272:       PetscInt         xep, yep, i, j, k, ind
    273:       PetscInt         xs, ys, xm, ym
    274:       PetscInt         gxs, gys, gxm, gym
    275:       PetscInt         i1

    277:       i1 = 1
    278:       0
    279:       cdiv3 = param/3.0
    280:       zero = 0.0
    281:       p5   = 0.5
    282:       hx = 1.0/real(mx + 1)
    283:       hy = 1.0/real(my + 1)
    284:       fquad = zero
    285:       flin = zero

    287: !  Initialize gradient to zero
    288:       call VecSet(G,zero,ierr)

    290: !  Scatter ghost points to local vector
    291:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
    292:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

    295: !  Get corner information
    296:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    297:      &                  PETSC_NULL_INTEGER,ierr)
    298:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
    299:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)

    301: !  Get pointer to vector data.
    302:       call VecGetArray(localX,lx_v,lx_i,ierr)

    305: !  Set local loop dimensions
    306:       xe = xs+xm
    307:       ye = ys+ym
    308:       if (xs .eq. 0) then
    309:          xsm = xs-1
    310:       else
    311:          xsm = xs
    312:       endif
    313:       if (ys .eq. 0) then
    314:          ysm = ys-1
    315:       else
    316:          ysm = ys
    317:       endif
    318:       if (xe .eq. mx) then
    319:          xep = xe+1
    320:       else
    321:          xep = xe
    322:       endif
    323:       if (ye .eq. my) then
    324:          yep = ye+1
    325:       else
    326:          yep = ye
    327:       endif

    329: !     Compute local gradient contributions over the lower triangular elements

    331:       do j = ysm, ye-1
    332:          do i = xsm, xe-1
    333:             k  = (j-gys)*gxm + i-gxs
    334:             v  = zero
    335:             vr = zero
    336:             vt = zero
    337:             if (i .ge. 0 .and. j .ge. 0)      v = lx_v(lx_i+k)
    338:             if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
    339:             if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
    340:             dvdx = (vr-v)/hx
    341:             dvdy = (vt-v)/hy
    342:             if (i .ne. -1 .and. j .ne. -1) then
    343:                ind = k
    344:                val = - dvdx/hx - dvdy/hy - cdiv3
    345:                call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
    346:             endif
    347:             if (i .ne. mx-1 .and. j .ne. -1) then
    348:                ind = k+1
    349:                val =  dvdx/hx - cdiv3
    350:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
    351:             endif
    352:             if (i .ne. -1 .and. j .ne. my-1) then
    353:               ind = k+gxm
    354:               val = dvdy/hy - cdiv3
    355:               call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
    356:             endif
    357:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
    358:             flin = flin - cdiv3 * (v+vr+vt)
    359:          end do
    360:       end do

    362: !     Compute local gradient contributions over the upper triangular elements

    364:       do j = ys, yep-1
    365:          do i = xs, xep-1
    366:             k  = (j-gys)*gxm + i-gxs
    367:             vb = zero
    368:             vl = zero
    369:             v  = zero
    370:             if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
    371:             if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
    372:             if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
    373:             dvdx = (v-vl)/hx
    374:             dvdy = (v-vb)/hy
    375:             if (i .ne. mx .and. j .ne. 0) then
    376:                ind = k-gxm
    377:                val = - dvdy/hy - cdiv3
    378:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
    379:             endif
    380:             if (i .ne. 0 .and. j .ne. my) then
    381:                ind = k-1
    382:                val =  - dvdx/hx - cdiv3
    383:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
    384:             endif
    385:             if (i .ne. mx .and. j .ne. my) then
    386:                ind = k
    387:                val =  dvdx/hx + dvdy/hy - cdiv3
    388:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
    389:             endif
    390:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
    391:             flin = flin - cdiv3*(vb + vl + v)
    392:          end do
    393:       end do

    395: !  Restore vector
    396:       call VecRestoreArray(localX,lx_v,lx_i,ierr)

    398: !  Assemble gradient vector
    399:       call VecAssemblyBegin(G,ierr)
    400:       call VecAssemblyEnd(G,ierr)

    402: ! Scale the gradient
    403:       area = p5*hx*hy
    404:       floc = area *(p5*fquad+flin)
    405:       call VecScale(G,area,ierr)

    408: !  Sum function contributions from all processes
    409:       call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,                   &
    410:      &                   PETSC_COMM_WORLD,ierr)
    411:       call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+                        &
    412:      &                   16.0d0*(xep-xs)*(yep-ys),ierr)
    413:       return
    414:       end

    416:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
    417:       use mymodule
    418:       implicit none

    420:       Tao       tao
    421:       Vec             X
    422:       Mat             H,Hpre
    423:       PetscErrorCode  ierr
    424:       PetscInt        dummy

    427:       PetscInt i,j,k
    428:       PetscInt col(0:4),row
    429:       PetscInt xs,xm,gxs,gxm
    430:       PetscInt ys,ym,gys,gym
    431:       PetscReal v(0:4)
    432:       PetscInt i1

    434:       i1 = 1

    436: !     Get local grid boundaries
    437:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    438:      &                PETSC_NULL_INTEGER,ierr)
    439:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,      &
    440:      &                PETSC_NULL_INTEGER,ierr)

    442:       do j=ys,ys+ym-1
    443:          do i=xs,xs+xm-1
    444:             row = (j-gys)*gxm + (i-gxs)

    446:             k = 0
    447:             if (j .gt. gys) then
    448:                v(k) = -1.0
    449:                col(k) = row-gxm
    450:                k = k + 1
    451:             endif

    453:             if (i .gt. gxs) then
    454:                v(k) = -1.0
    455:                col(k) = row - 1
    456:                k = k +1
    457:             endif

    459:             v(k) = 4.0
    460:             col(k) = row
    461:             k = k + 1

    463:             if (i+1 .lt. gxs + gxm) then
    464:                v(k) = -1.0
    465:                col(k) = row + 1
    466:                k = k + 1
    467:             endif

    469:             if (j+1 .lt. gys + gym) then
    470:                v(k) = -1.0
    471:                col(k) = row + gxm
    472:                k = k + 1
    473:             endif

    475:             call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
    476:          enddo
    477:       enddo

    480: !     Assemble matrix
    481:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
    482:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)

    485: !     Tell the matrix we will never add a new nonzero location to the
    486: !     matrix.  If we do it will generate an error.

    488:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
    489:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

    492:       call PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)

    494:       0
    495:       return
    496:       end

    500:       subroutine Monitor(tao, dummy, ierr)
    501:       use mymodule
    502:       implicit none

    504:       Tao tao
    505:       PetscInt dummy
    506:       PetscErrorCode ierr

    508:       PetscInt its
    509:       PetscReal f,gnorm,cnorm,xdiff
    510:       TaoConvergedReason reason

    512:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,             &
    513:      &     reason,ierr)
    514:       if (mod(its,5) .ne. 0) then
    515:          call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',  &
    516:      &        ierr)
    517:       endif

    519:       0

    521:       return
    522:       end

    524:       subroutine ConvergenceTest(tao, dummy, ierr)
    525:       use mymodule
    526:       implicit none

    528:       Tao tao
    529:       PetscInt dummy
    530:       PetscErrorCode ierr

    532:       PetscInt its
    533:       PetscReal f,gnorm,cnorm,xdiff
    534:       TaoConvergedReason reason

    536:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,            &
    537:      &     reason,ierr)
    538:       if (its .eq. 7) then
    539:        call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
    540:       endif

    542:       0

    544:       return
    545:       end

    547: !/*TEST
    548: !
    549: !   build:
    550: !      requires: !complex
    551: !
    552: !   test:
    553: !      args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
    554: !
    555: !   test:
    556: !      suffix: 2
    557: !      nsize: 2
    558: !      args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
    559: !
    560: !   test:
    561: !      suffix: 3
    562: !      args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
    563: !TEST*/