Actual source code: eptorsion2f.F

      1: ! "$Id: eptorsion2f.F 1.52 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
      2: !
      3: !  Program usage: mpirun -np <proc> eptorsion2f [all TAO options]
      4: !
      5: !  Description:  This example demonstrates use of the TAO package to solve
      6: !  unconstrained minimization problems in parallel.  This example is based
      7: !  on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.
      8: !  The command line options are:
      9: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
     10: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
     11: !    -par <param>, where <param> = angle of twist per unit length
     12: !
     13: !/*T
     14: !   Concepts: TAO - Solving an unconstrained minimization problem
     15: !   Routines: TaoInitialize(); TaoFinalize();
     16: !   Routines: TaoCreate(); TaoDestroy();
     17: !   Routines: TaoApplicationCreate(); TaoAppDestroy();
     18: !   Routines: TaoAppSetObjectiveAndGradientRoutine();
     19: !   Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
     20: !   Routines: TaoSetApplication(); TaoSetOptions();
     21: !   Routines: TaoAppSetInitialSolutionVec(); TaoSolveApplication(); TaoDestroy();
     22: !   Routines: TaoGetSolutionStatus();
     23: !   Processors: n
     24: !T*/
     25: !
     26: ! ----------------------------------------------------------------------
     27: !
     28: !  Elastic-plastic torsion problem.
     29: !
     30: !  The elastic plastic torsion problem arises from the determination
     31: !  of the stress field on an infinitely long cylindrical bar, which is
     32: !  equivalent to the solution of the following problem:
     33: !     min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
     34: !  where C is the torsion angle per unit length.
     35: !
     36: !  The C version of this code is eptorsion2.c
     37: !
     38: ! ----------------------------------------------------------------------

     40:       implicit none
     41: #include "eptorsion2f.h"

     43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     44: !                   Variable declarations
     45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     46: !
     47: !  See additional variable declarations in the file eptorsion2f.h
     48: !
     49:       integer          info           ! used to check for functions returning nonzeros
     50:       Vec              x              ! solution vector
     51:       Mat              H              ! hessian matrix
     52:       integer          Nx, Ny         ! number of processes in x- and y- directions
     53:       TAO_SOLVER       tao            ! TAO_SOLVER solver context
     54:       TAO_APPLICATION  torsionapp     ! TAO application context (PETSc)
     55:       TaoTerminateReason reason
     56:       PetscTruth       flg
     57:       integer          iter           ! iteration information
     58:       PetscScalar      ff,gnorm,cnorm,xdiff
     59:

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

     64:       external FormInitialGuess,FormFunctionGradient,ComputeHessian

     66: !     Initialize TAO, PETSc  contexts
     67:       call PetscInitialize(PETSC_NULL_CHARACTER,info)
     68:       call TaoInitialize(PETSC_NULL_CHARACTER,info)

     70: !     Specify default parameters
     71:       param = 5.0d0
     72:       mx = 10
     73:       my = 10
     74:       Nx = PETSC_DECIDE
     75:       Ny = PETSC_DECIDE

     77: !     Check for any command line arguments that might override defaults
     78:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mx",mx,flg,info)
     79:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,info)
     80:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-par",              &
     81:      &                         param,flg,info)

     83:
     84: !     Set up distributed array and vectors
     85:       call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,      &
     86:      &     mx,my,Nx,Ny,1,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,        &
     87:      &     da,info)

     89: !     Create vectors
     90:       call DACreateGlobalVector(da,x,info)
     91:       call DACreateLocalVector(da,localX,info)

     93: !     Create Hessian
     94:       call DAGetMatrix(da,MATAIJ,H,info)
     95:       call MatSetOption(H,MAT_SYMMETRIC,info)

     97: !     The TAO code begins here

     99: !     Create TAO solver
    100:       call TaoCreate(MPI_COMM_WORLD,'tao_cg_fr',tao,info)
    101:       call TaoApplicationCreate(MPI_COMM_WORLD,torsionapp,info)

    103: !     Set routines for function and gradient evaluation

    105: !     TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
    106:       call TaoAppSetObjectiveAndGradientRo(torsionapp,                       &
    107:      &     FormFunctionGradient,TAO_NULL_OBJECT,info)
    108:       call TaoAppSetHessianMat(torsionapp,H,H,info)
    109:       call TaoAppSetHessianRoutine(torsionapp,ComputeHessian,                &
    110:      &    TAO_NULL_OBJECT,info)

    112: !     Set initial guess
    113:       call FormInitialGuess(x,info)
    114:       call TaoAppSetInitialSolutionVec(torsionapp,x,info)

    116: !     Check for any TAO command line options
    117:       call TaoSetOptions(torsionapp, tao,info)

    119: !     Now that the PETSc application is all set, attach to tao context
    120:       call TaoSetApplication(tao,torsionapp,info)

    123: !     SOLVE THE APPLICATION
    124:       call TaoSolveApplication(torsionapp,tao,info)

    126: !     Get information on termination
    127:       call TaoGetSolutionStatus(tao,iter,ff,gnorm,cnorm,xdiff,           &
    128:      &                          reason,info)
    129:       if (reason .lt. 0) then
    130:          print *,'TAO did not terminate successfully'
    131:          call TaoView(tao,info)
    132:       endif
    133:
    134: !     Free TAO data structures
    135:       call TaoDestroy(tao,info)
    136:       call TaoAppDestroy(torsionapp,info);

    138:
    139: !     Free PETSc data structures
    140:       call VecDestroy(x,info)
    141:       call VecDestroy(localX,info)
    142:       call MatDestroy(H,info)
    143:       call DADestroy(da,info)

    146: !     Finalize TAO and PETSc
    147:       call PetscFinalize(info)
    148:       call TaoFinalize(info)

    150:       end

    153: ! ---------------------------------------------------------------------
    154: !
    155: !   FormInitialGuess - Computes an initial approximation to the solution.
    156: !
    157: !   Input Parameters:
    158: !   X    - vector
    159: !
    160: !   Output Parameters:
    161: !   X    - vector
    162: !   info - error code
    163: !
    164:       subroutine FormInitialGuess(X,info)
    165:       implicit none

    167: ! mx, my are defined in eptorsion2f.h
    168: #include "eptorsion2f.h"

    170: !  Input/output variables:
    171:       Vec              X
    172:       integer          info

    174: !  Local variables:
    175:       integer          i, j, k, xe, ye
    176:       PetscScalar      temp, val, hx, hy
    177:       integer          xs, ys, xm, ym, gxm, gym, gxs, gys

    179:       hx = 1.0d0/(mx + 1)
    180:       hy = 1.0d0/(my + 1)

    182: !  Get corner information
    183:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    184:      &                  PETSC_NULL_INTEGER,info)
    185:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,              &
    186:      &                   gxm,gym,PETSC_NULL_INTEGER,info)

    190: !  Compute initial guess over locally owned part of mesh
    191:       xe = xs+xm
    192:       ye = ys+ym
    193:       do j=ys,ye-1
    194:          temp = min(j+1,my-j)*hy
    195:          do i=xs,xe-1
    196:             k   = (j-gys)*gxm + i-gxs
    197:             val = min((min(i+1,mx-i))*hx,temp)
    198:             call VecSetValuesLocal(X,1,k,val,ADD_VALUES,info)
    199:          end do
    200:       end do

    202:       return
    203:       end

    206: ! ---------------------------------------------------------------------
    207: !
    208: !  FormFunctionGradient - Evaluates gradient G(X).
    209: !
    210: !  Input Parameters:
    211: !  tao   - the TAO_SOLVER context
    212: !  X     - input vector
    213: !  dummy - optional user-defined context (not used here)
    214: !
    215: !  Output Parameters:
    216: !  f     - the function value at X
    217: !  G     - vector containing the newly evaluated gradient
    218: !  info  - error code
    219: !
    220: !  Notes:
    221: !  This routine serves as a wrapper for the lower-level routine
    222: !  "ApplicationGradient", where the actual computations are
    223: !  done using the standard Fortran style of treating the local
    224: !  input vector data as an array over the local mesh.
    225: !
    226:       subroutine FormFunctionGradient(taoapp,X,f,G,dummy,info)
    227:       implicit none

    229: ! da, mx, my, param, localX declared in eptorsion2f.h
    230: #include "eptorsion2f.h"

    232: !  Input/output variables:
    233:       TAO_APPLICATION  taoapp
    234:       Vec              X, G
    235:       PetscScalar      f
    236:       integer          dummy, info

    238: !  Declarations for use with local array:

    241: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    242: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
    243: ! will return an array of doubles referenced by x_array offset by x_index.
    244: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    245: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    246:       PetscScalar      lx_v(0:1)
    247:       PetscOffset      lx_i

    249: !  Local variables:
    250:       PetscScalar      zero, p5, area, cdiv3, val, flin, fquad,floc
    251:       PetscScalar      v, vb, vl, vr, vt, dvdx, dvdy, hx, hy
    252:       integer          xe, ye, xsm, ysm, xep, yep, i, j, k, ind
    253:       integer          xs, ys, xm, ym, gxs, gys, gxm, gym

    255:       info  = 0
    256:       cdiv3 = param/3.0d0
    257:       zero = 0.0d0
    258:       p5   = 0.5d0
    259:       hx = 1.0d0/(mx + 1)
    260:       hy = 1.0d0/(my + 1)
    261:       fquad = zero
    262:       flin = zero

    264: !  Initialize gradient to zero
    265:       call VecSet(G,zero,info)

    267: !  Scatter ghost points to local vector
    268:       call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
    269:       call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)

    272: !  Get corner information
    273:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    274:      &                  PETSC_NULL_INTEGER,info)
    275:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,              &
    276:      &                   gxm,gym,PETSC_NULL_INTEGER,info)

    278: !  Get pointer to vector data.
    279:       call VecGetArray(localX,lx_v,lx_i,info)

    282: !  Set local loop dimensions
    283:       xe = xs+xm
    284:       ye = ys+ym
    285:       if (xs .eq. 0) then
    286:          xsm = xs-1
    287:       else
    288:          xsm = xs
    289:       endif
    290:       if (ys .eq. 0) then
    291:          ysm = ys-1
    292:       else
    293:          ysm = ys
    294:       endif
    295:       if (xe .eq. mx) then
    296:          xep = xe+1
    297:       else
    298:          xep = xe
    299:       endif
    300:       if (ye .eq. my) then
    301:          yep = ye+1
    302:       else
    303:          yep = ye
    304:       endif

    306: !     Compute local gradient contributions over the lower triangular elements
    307:
    308:       do j = ysm, ye-1
    309:          do i = xsm, xe-1
    310:             k  = (j-gys)*gxm + i-gxs
    311:             v  = zero
    312:             vr = zero
    313:             vt = zero
    314:             if (i .ge. 0 .and. j .ge. 0)      v = lx_v(lx_i+k)
    315:             if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
    316:             if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
    317:             dvdx = (vr-v)/hx
    318:             dvdy = (vt-v)/hy
    319:             if (i .ne. -1 .and. j .ne. -1) then
    320:                ind = k
    321:                val = - dvdx/hx - dvdy/hy - cdiv3
    322:                call VecSetValuesLocal(G,1,k,val,ADD_VALUES,info)
    323:             endif
    324:             if (i .ne. mx-1 .and. j .ne. -1) then
    325:                ind = k+1
    326:                val =  dvdx/hx - cdiv3
    327:                call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
    328:             endif
    329:             if (i .ne. -1 .and. j .ne. my-1) then
    330:               ind = k+gxm
    331:               val = dvdy/hy - cdiv3
    332:               call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
    333:             endif
    334:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
    335:             flin = flin - cdiv3 * (v+vr+vt)
    336:          end do
    337:       end do

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

    341:       do j = ys, yep-1
    342:          do i = xs, xep-1
    343:             k  = (j-gys)*gxm + i-gxs
    344:             vb = zero
    345:             vl = zero
    346:             v  = zero
    347:             if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
    348:             if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
    349:             if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
    350:             dvdx = (v-vl)/hx
    351:             dvdy = (v-vb)/hy
    352:             if (i .ne. mx .and. j .ne. 0) then
    353:                ind = k-gxm
    354:                val = - dvdy/hy - cdiv3
    355:                call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
    356:             endif
    357:             if (i .ne. 0 .and. j .ne. my) then
    358:                ind = k-1
    359:                val =  - dvdx/hx - cdiv3
    360:                call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
    361:             endif
    362:             if (i .ne. mx .and. j .ne. my) then
    363:                ind = k
    364:                val =  dvdx/hx + dvdy/hy - cdiv3
    365:                call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
    366:             endif
    367:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
    368:             flin = flin - cdiv3*(vb + vl + v)
    369:          end do
    370:       end do

    372: !  Restore vector
    373:       call VecRestoreArray(localX,lx_v,lx_i,info)

    375: !  Assemble gradient vector
    376:       call VecAssemblyBegin(G,info)
    377:       call VecAssemblyEnd(G,info)

    379: ! Scale the gradient
    380:       area = p5*hx*hy
    381:       floc = area *(p5*fquad+flin)
    382:       call VecScale(G,area,info)

    385: !  Sum function contributions from all processes
    386:       call MPI_Allreduce(floc,f,1,MPI_DOUBLE_PRECISION,MPI_SUM,           &
    387:      &                   MPI_COMM_WORLD,info)

    389:       call PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16,       &
    390:      &                    info)

    393:       return
    394:       end

    399:       subroutine ComputeHessian(taoapp, X, H, Hpre, flag, dummy, info)
    400:       implicit none
    401: #include "eptorsion2f.h"
    402:       TAO_APPLICATION taoapp
    403:       Vec             X
    404:       Mat             H,Hpre
    405:       MatStructure    flag
    406:       integer         dummy,info

    408:
    409:       integer i,j,k
    410:       integer col(0:4),row
    411:       integer xs,xm,gxs,gxm,ys,ym,gys,gym
    412:       PetscScalar v(0:4)

    414: !     Get local grid boundaries
    415:       call DAGetCorners(da,xs,ys,TAO_NULL_INTEGER,xm,ym,                &
    416:      &                TAO_NULL_INTEGER,info)
    417:       call DAGetGhostCorners(da,gxs,gys,TAO_NULL_INTEGER,gxm,gym,        &
    418:      &                TAO_NULL_INTEGER,info)

    420:       do j=ys,ys+ym-1
    421:          do i=xs,xs+xm-1
    422:             row = (j-gys)*gxm + (i-gxs)

    424:             k = 0
    425:             if (j .gt. gys) then
    426:                v(k) = -1.0d0
    427:                col(k) = row-gxm
    428:                k = k + 1
    429:             endif

    431:             if (i .gt. gxs) then
    432:                v(k) = -1.0d0
    433:                col(k) = row - 1
    434:                k = k +1
    435:             endif

    437:             v(k) = 4.0d0
    438:             col(k) = row
    439:             k = k + 1

    441:             if (i+1 .lt. gxs + gxm) then
    442:                v(k) = -1.0d0
    443:                col(k) = row + 1
    444:                k = k + 1
    445:             endif

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

    453:             call MatSetValuesLocal(H,1,row,k,col,v,INSERT_VALUES,info)
    454:          enddo
    455:       enddo

    457:
    458: !     Assemble matrix
    459:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
    460:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)

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

    466:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,info)
    467:       call MatSetOption(H,MAT_SYMMETRIC,info)

    470:       call PetscLogFlops(9*xm*ym + 49*xm,info)

    472:       info = 0
    473:       return
    474:       end
    475:
    476:
    477: