Actual source code: plate2f.F90

      1: !  Program usage: mpiexec -n <proc> plate2f [all TAO options]
      2: !
      3: !  This example demonstrates use of the TAO package to solve a bound constrained
      4: !  minimization problem.  This example is based on a problem from the
      5: !  MINPACK-2 test suite.  Given a rectangular 2-D domain and boundary values
      6: !  along the edges of the domain, the objective is to find the surface
      7: !  with the minimal area that satisfies the boundary conditions.
      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: !    -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
     12: !    -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
     13: !    -bheight <ht>, where <ht> = height of the plate
     14: !
     15: !!/*T
     16: !   Concepts: TAO^Solving a bound constrained minimization problem
     17: !   Routines: TaoCreate();
     18: !   Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
     19: !   Routines: TaoSetHessianRoutine();
     20: !   Routines: TaoSetVariableBoundsRoutine();
     21: !   Routines: TaoSetInitialVector();
     22: !   Routines: TaoSetFromOptions();
     23: !   Routines: TaoSolve();
     24: !   Routines: TaoDestroy();
     25: !   Processors: n
     26: !T*/

     28:       module mymodule
     29: #include "petsc/finclude/petscdmda.h"
     30: #include "petsc/finclude/petsctao.h"
     31:       use petscdmda
     32:       use petsctao

     34:       Vec              localX, localV
     35:       Vec              Top, Left
     36:       Vec              Right, Bottom
     37:       DM               dm
     38:       PetscReal      bheight
     39:       PetscInt         bmx, bmy
     40:       PetscInt         mx, my, Nx, Ny, N
     41:       end module

     44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     45: !                   Variable declarations
     46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     47: !
     48: !  Variables:
     49: !    (common from plate2f.h):
     50: !    Nx, Ny           number of processors in x- and y- directions
     51: !    mx, my           number of grid points in x,y directions
     52: !    N    global dimension of vector
     53:       use mymodule
     54:       implicit none

     56:       PetscErrorCode   ierr          ! used to check for functions returning nonzeros
     57:       Vec              x             ! solution vector
     58:       PetscInt         m             ! number of local elements in vector
     59:       Tao              tao           ! Tao solver context
     60:       Mat              H             ! Hessian matrix
     61:       ISLocalToGlobalMapping isltog  ! local to global mapping object
     62:       PetscBool        flg
     63:       PetscInt         i1,i3,i7

     66:       external FormFunctionGradient
     67:       external FormHessian
     68:       external MSA_BoundaryConditions
     69:       external MSA_Plate
     70:       external MSA_InitialPoint
     71: ! Initialize Tao

     73:       i1=1
     74:       i3=3
     75:       i7=7

     78:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     79:       if (ierr .ne. 0) then
     80:         print*,'Unable to initialize PETSc'
     81:         stop
     82:       endif

     84: ! Specify default dimensions of the problem
     85:       mx = 10
     86:       my = 10
     87:       bheight = 0.1

     89: ! Check for any command line arguments that override defaults

     91:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,     &
     92:      &                        '-mx',mx,flg,ierr)
     93:       CHKERRA(ierr)
     94:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,     &
     95:      &                        '-my',my,flg,ierr)
     96:       CHKERRA(ierr)

     98:       bmx = mx/2
     99:       bmy = my/2

    101:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,     &
    102:      &                        '-bmx',bmx,flg,ierr)
    103:       CHKERRA(ierr)
    104:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,     &
    105:      &                        '-bmy',bmy,flg,ierr)
    106:       CHKERRA(ierr)
    107:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,    &
    108:      &                         '-bheight',bheight,flg,ierr)
    109:       CHKERRA(ierr)

    112: ! Calculate any derived values from parameters
    113:       N = mx*my

    115: ! Let Petsc determine the dimensions of the local vectors
    116:       Nx = PETSC_DECIDE
    117:       NY = PETSC_DECIDE

    119: ! A two dimensional distributed array will help define this problem, which
    120: ! derives from an elliptic PDE on a two-dimensional domain.  From the
    121: ! distributed array, create the vectors

    123:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,                    &
    124:      &     DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,                              &
    125:      &     mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,           &
    126:      &     dm,ierr)
    127:       CHKERRA(ierr)
    128:       call DMSetFromOptions(dm,ierr)
    129:       call DMSetUp(dm,ierr)

    131: ! Extract global and local vectors from DM; The local vectors are
    132: ! used solely as work space for the evaluation of the function,
    133: ! gradient, and Hessian.  Duplicate for remaining vectors that are
    134: ! the same types.

    136:       call DMCreateGlobalVector(dm,x,ierr)
    137:       CHKERRA(ierr)
    138:       call DMCreateLocalVector(dm,localX,ierr)
    139:       CHKERRA(ierr)
    140:       call VecDuplicate(localX,localV,ierr)
    141:       CHKERRA(ierr)

    143: ! Create a matrix data structure to store the Hessian.
    144: ! Here we (optionally) also associate the local numbering scheme
    145: ! with the matrix so that later we can use local indices for matrix
    146: ! assembly

    148:       call VecGetLocalSize(x,m,ierr)
    149:       CHKERRA(ierr)
    150:       call MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER,   &
    151:      &     i3,PETSC_NULL_INTEGER,H,ierr)
    152:       CHKERRA(ierr)

    154:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
    155:       CHKERRA(ierr)
    156:       call DMGetLocalToGlobalMapping(dm,isltog,ierr)
    157:       CHKERRA(ierr)
    158:       call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
    159:       CHKERRA(ierr)

    162: ! The Tao code begins here
    163: ! Create TAO solver and set desired solution method.
    164: ! This problems uses bounded variables, so the
    165: ! method must either be 'tao_tron' or 'tao_blmvm'

    167:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
    168:       CHKERRA(ierr)
    169:       call TaoSetType(tao,TAOBLMVM,ierr)
    170:       CHKERRA(ierr)

    172: !     Set minimization function and gradient, hessian evaluation functions

    174:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
    175:      &     FormFunctionGradient,0,ierr)
    176:       CHKERRA(ierr)

    178:       call TaoSetHessianRoutine(tao,H,H,FormHessian,                    &
    179:      &     0, ierr)
    180:       CHKERRA(ierr)

    182: ! Set Variable bounds
    183:       call MSA_BoundaryConditions(ierr)
    184:       CHKERRA(ierr)
    185:       call TaoSetVariableBoundsRoutine(tao,MSA_Plate,                   &
    186:      &     0,ierr)
    187:       CHKERRA(ierr)

    189: ! Set the initial solution guess
    190:       call MSA_InitialPoint(x, ierr)
    191:       CHKERRA(ierr)
    192:       call TaoSetInitialVector(tao,x,ierr)
    193:       CHKERRA(ierr)

    195: ! Check for any tao command line options
    196:       call TaoSetFromOptions(tao,ierr)
    197:       CHKERRA(ierr)

    199: ! Solve the application
    200:       call TaoSolve(tao,ierr)
    201:       CHKERRA(ierr)

    203: ! Free TAO data structures
    204:       call TaoDestroy(tao,ierr)
    205:       CHKERRA(ierr)

    207: ! Free PETSc data structures
    208:       call VecDestroy(x,ierr)
    209:       call VecDestroy(Top,ierr)
    210:       call VecDestroy(Bottom,ierr)
    211:       call VecDestroy(Left,ierr)
    212:       call VecDestroy(Right,ierr)
    213:       call MatDestroy(H,ierr)
    214:       call VecDestroy(localX,ierr)
    215:       call VecDestroy(localV,ierr)
    216:       call DMDestroy(dm,ierr)

    218: ! Finalize TAO

    220:       call PetscFinalize(ierr)

    222:       end

    224: ! ---------------------------------------------------------------------
    225: !
    226: !  FormFunctionGradient - Evaluates function f(X).
    227: !
    228: !  Input Parameters:
    229: !  tao   - the Tao context
    230: !  X     - the input vector
    231: !  dummy - optional user-defined context, as set by TaoSetFunction()
    232: !          (not used here)
    233: !
    234: !  Output Parameters:
    235: !  fcn     - the newly evaluated function
    236: !  G       - the gradient vector
    237: !  info  - error code
    238: !

    241:       subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
    242:       use mymodule
    243:       implicit none

    245: ! Input/output variables

    247:       Tao        tao
    248:       PetscReal      fcn
    249:       Vec              X, G
    250:       PetscErrorCode   ierr
    251:       PetscInt         dummy

    253:       PetscInt         i,j,row
    254:       PetscInt         xs, xm
    255:       PetscInt         gxs, gxm
    256:       PetscInt         ys, ym
    257:       PetscInt         gys, gym
    258:       PetscReal      ft,zero,hx,hy,hydhx,hxdhy
    259:       PetscReal      area,rhx,rhy
    260:       PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3
    261:       PetscReal      d4,d5,d6,d7,d8
    262:       PetscReal      df1dxc,df2dxc,df3dxc,df4dxc
    263:       PetscReal      df5dxc,df6dxc
    264:       PetscReal      xc,xl,xr,xt,xb,xlt,xrb

    267: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    268: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
    269: ! will return an array of doubles referenced by x_array offset by x_index.
    270: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    271: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    272:       PetscReal      g_v(0:1),x_v(0:1)
    273:       PetscReal      top_v(0:1),left_v(0:1)
    274:       PetscReal      right_v(0:1),bottom_v(0:1)
    275:       PetscOffset      g_i,left_i,right_i
    276:       PetscOffset      bottom_i,top_i,x_i

    278:       ft = 0.0
    279:       zero = 0.0
    280:       hx = 1.0/real(mx + 1)
    281:       hy = 1.0/real(my + 1)
    282:       hydhx = hy/hx
    283:       hxdhy = hx/hy
    284:       area = 0.5 * hx * hy
    285:       rhx = real(mx) + 1.0
    286:       rhy = real(my) + 1.0

    289: ! Get local mesh boundaries
    290:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
    291:      &                  PETSC_NULL_INTEGER,ierr)
    292:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,             &
    293:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

    295: ! Scatter ghost points to local vector
    296:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
    297:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

    299: ! Initialize the vector to zero
    300:       call VecSet(localV,zero,ierr)

    302: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
    303:       call VecGetArray(localX,x_v,x_i,ierr)
    304:       call VecGetArray(localV,g_v,g_i,ierr)
    305:       call VecGetArray(Top,top_v,top_i,ierr)
    306:       call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
    307:       call VecGetArray(Left,left_v,left_i,ierr)
    308:       call VecGetArray(Right,right_v,right_i,ierr)

    310: ! Compute function over the locally owned part of the mesh
    311:       do j = ys,ys+ym-1
    312:          do i = xs,xs+xm-1
    313:             row = (j-gys)*gxm + (i-gxs)
    314:             xc = x_v(row+x_i)
    315:             xt = xc
    316:             xb = xc
    317:             xr = xc
    318:             xl = xc
    319:             xrb = xc
    320:             xlt = xc

    322:             if (i .eq. 0) then !left side
    323:                xl = left_v(j - ys + 1 + left_i)
    324:                xlt = left_v(j - ys + 2 + left_i)
    325:             else
    326:                xl = x_v(row - 1 + x_i)
    327:             endif

    329:             if (j .eq. 0) then !bottom side
    330:                xb = bottom_v(i - xs + 1 + bottom_i)
    331:                xrb = bottom_v(i - xs + 2 + bottom_i)
    332:             else
    333:                xb = x_v(row - gxm + x_i)
    334:             endif

    336:             if (i + 1 .eq. gxs + gxm) then !right side
    337:                xr = right_v(j - ys + 1 + right_i)
    338:                xrb = right_v(j - ys + right_i)
    339:             else
    340:                xr = x_v(row + 1 + x_i)
    341:             endif

    343:             if (j + 1 .eq. gys + gym) then !top side
    344:                xt = top_v(i - xs + 1 + top_i)
    345:                xlt = top_v(i - xs + top_i)
    346:             else
    347:                xt = x_v(row + gxm + x_i)
    348:             endif

    350:             if ((i .gt. gxs) .and. (j + 1 .lt. gys + gym)) then
    351:                xlt = x_v(row - 1 + gxm + x_i)
    352:             endif

    354:             if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
    355:                xrb = x_v(row + 1 - gxm + x_i)
    356:             endif

    358:             d1 = xc-xl
    359:             d2 = xc-xr
    360:             d3 = xc-xt
    361:             d4 = xc-xb
    362:             d5 = xr-xrb
    363:             d6 = xrb-xb
    364:             d7 = xlt-xl
    365:             d8 = xt-xlt

    367:             df1dxc = d1 * hydhx
    368:             df2dxc = d1 * hydhx + d4 * hxdhy
    369:             df3dxc = d3 * hxdhy
    370:             df4dxc = d2 * hydhx + d3 * hxdhy
    371:             df5dxc = d2 * hydhx
    372:             df6dxc = d4 * hxdhy

    374:             d1 = d1 * rhx
    375:             d2 = d2 * rhx
    376:             d3 = d3 * rhy
    377:             d4 = d4 * rhy
    378:             d5 = d5 * rhy
    379:             d6 = d6 * rhx
    380:             d7 = d7 * rhy
    381:             d8 = d8 * rhx

    383:             f1 = sqrt(1.0 + d1*d1 + d7*d7)
    384:             f2 = sqrt(1.0 + d1*d1 + d4*d4)
    385:             f3 = sqrt(1.0 + d3*d3 + d8*d8)
    386:             f4 = sqrt(1.0 + d3*d3 + d2*d2)
    387:             f5 = sqrt(1.0 + d2*d2 + d5*d5)
    388:             f6 = sqrt(1.0 + d4*d4 + d6*d6)

    390:             ft = ft + f2 + f4

    392:             df1dxc = df1dxc / f1
    393:             df2dxc = df2dxc / f2
    394:             df3dxc = df3dxc / f3
    395:             df4dxc = df4dxc / f4
    396:             df5dxc = df5dxc / f5
    397:             df6dxc = df6dxc / f6

    399:             g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc +  &
    400:      &                              df5dxc + df6dxc)
    401:          enddo
    402:       enddo

    404: ! Compute triangular areas along the border of the domain.
    405:       if (xs .eq. 0) then  ! left side
    406:          do j=ys,ys+ym-1
    407:             d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i))         &
    408:      &                 * rhy
    409:             d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i))        &
    410:      &                 * rhx
    411:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
    412:          enddo
    413:       endif

    416:       if (ys .eq. 0) then !bottom side
    417:          do i=xs,xs+xm-1
    418:             d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i))    &
    419:      &                    * rhx
    420:             d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
    421:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
    422:          enddo
    423:       endif

    426:       if (xs + xm .eq. mx) then ! right side
    427:          do j=ys,ys+ym-1
    428:             d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
    429:             d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
    430:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
    431:          enddo
    432:       endif

    435:       if (ys + ym .eq. my) then
    436:          do i=xs,xs+xm-1
    437:             d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
    438:             d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
    439:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
    440:          enddo
    441:       endif

    444:       if ((ys .eq. 0) .and. (xs .eq. 0)) then
    445:          d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
    446:          d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
    447:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
    448:       endif

    450:       if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
    451:          d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
    452:          d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
    453:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
    454:       endif

    456:       ft = ft * area
    457:       call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR,                            &
    458:      &             MPIU_SUM,PETSC_COMM_WORLD,ierr)

    461: ! Restore vectors
    462:       call VecRestoreArray(localX,x_v,x_i,ierr)
    463:       call VecRestoreArray(localV,g_v,g_i,ierr)
    464:       call VecRestoreArray(Left,left_v,left_i,ierr)
    465:       call VecRestoreArray(Top,top_v,top_i,ierr)
    466:       call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
    467:       call VecRestoreArray(Right,right_v,right_i,ierr)

    469: ! Scatter values to global vector
    470:       call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
    471:       call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)

    473:       call PetscLogFlops(70.0d0*xm*ym,ierr)

    475:       return
    476:       end  !FormFunctionGradient

    482: ! ----------------------------------------------------------------------------
    483: !
    484: !
    485: !   FormHessian - Evaluates Hessian matrix.
    486: !
    487: !   Input Parameters:
    488: !.  tao  - the Tao context
    489: !.  X    - input vector
    490: !.  dummy  - not used
    491: !
    492: !   Output Parameters:
    493: !.  Hessian    - Hessian matrix
    494: !.  Hpc    - optionally different preconditioning matrix
    495: !.  flag - flag indicating matrix structure
    496: !
    497: !   Notes:
    498: !   Due to mesh point reordering with DMs, we must always work
    499: !   with the local mesh points, and then transform them to the new
    500: !   global numbering with the local-to-global mapping.  We cannot work
    501: !   directly with the global numbers for the original uniprocessor mesh!
    502: !
    503: !      MatSetValuesLocal(), using the local ordering (including
    504: !         ghost points!)
    505: !         - Then set matrix entries using the local ordering
    506: !           by calling MatSetValuesLocal()

    508:       subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
    509:       use mymodule
    510:       implicit none

    512:       Tao     tao
    513:       Vec            X
    514:       Mat            Hessian,Hpc
    515:       PetscInt       dummy
    516:       PetscErrorCode ierr

    518:       PetscInt       i,j,k,row
    519:       PetscInt       xs,xm,gxs,gxm
    520:       PetscInt       ys,ym,gys,gym
    521:       PetscInt       col(0:6)
    522:       PetscReal    hx,hy,hydhx,hxdhy,rhx,rhy
    523:       PetscReal    f1,f2,f3,f4,f5,f6,d1,d2,d3
    524:       PetscReal    d4,d5,d6,d7,d8
    525:       PetscReal    xc,xl,xr,xt,xb,xlt,xrb
    526:       PetscReal    hl,hr,ht,hb,hc,htl,hbr

    528: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    529: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
    530: ! will return an array of doubles referenced by x_array offset by x_index.
    531: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    532: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    533:       PetscReal   right_v(0:1),left_v(0:1)
    534:       PetscReal   bottom_v(0:1),top_v(0:1)
    535:       PetscReal   x_v(0:1)
    536:       PetscOffset   x_i,right_i,left_i
    537:       PetscOffset   bottom_i,top_i
    538:       PetscReal   v(0:6)
    539:       PetscBool     assembled
    540:       PetscInt      i1

    542:       i1=1

    544: ! Set various matrix options
    545:       call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,              &
    546:      &                  PETSC_TRUE,ierr)

    548: ! Get local mesh boundaries
    549:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
    550:      &                  PETSC_NULL_INTEGER,ierr)
    551:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
    552:      &                       PETSC_NULL_INTEGER,ierr)

    554: ! Scatter ghost points to local vectors
    555:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
    556:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

    558: ! Get pointers to vector data (see note on Fortran arrays above)
    559:       call VecGetArray(localX,x_v,x_i,ierr)
    560:       call VecGetArray(Top,top_v,top_i,ierr)
    561:       call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
    562:       call VecGetArray(Left,left_v,left_i,ierr)
    563:       call VecGetArray(Right,right_v,right_i,ierr)

    565: ! Initialize matrix entries to zero
    566:       call MatAssembled(Hessian,assembled,ierr)
    567:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)

    570:       rhx = real(mx) + 1.0
    571:       rhy = real(my) + 1.0
    572:       hx = 1.0/rhx
    573:       hy = 1.0/rhy
    574:       hydhx = hy/hx
    575:       hxdhy = hx/hy
    576: ! compute Hessian over the locally owned part of the mesh

    578:       do  i=xs,xs+xm-1
    579:          do  j=ys,ys+ym-1
    580:             row = (j-gys)*gxm + (i-gxs)

    582:             xc = x_v(row + x_i)
    583:             xt = xc
    584:             xb = xc
    585:             xr = xc
    586:             xl = xc
    587:             xrb = xc
    588:             xlt = xc

    590:             if (i .eq. gxs) then   ! Left side
    591:                xl = left_v(left_i + j - ys + 1)
    592:                xlt = left_v(left_i + j - ys + 2)
    593:             else
    594:                xl = x_v(x_i + row -1)
    595:             endif

    597:             if (j .eq. gys) then ! bottom side
    598:                xb = bottom_v(bottom_i + i - xs + 1)
    599:                xrb = bottom_v(bottom_i + i - xs + 2)
    600:             else
    601:                xb = x_v(x_i + row - gxm)
    602:             endif

    604:             if (i+1 .eq. gxs + gxm) then !right side
    605:                xr = right_v(right_i + j - ys + 1)
    606:                xrb = right_v(right_i + j - ys)
    607:             else
    608:                xr = x_v(x_i + row + 1)
    609:             endif

    611:             if (j+1 .eq. gym+gys) then !top side
    612:                xt = top_v(top_i +i - xs + 1)
    613:                xlt = top_v(top_i + i - xs)
    614:             else
    615:                xt = x_v(x_i + row + gxm)
    616:             endif

    618:             if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
    619:                xlt = x_v(x_i + row - 1 + gxm)
    620:             endif

    622:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
    623:                xrb = x_v(x_i + row + 1 - gxm)
    624:             endif

    626:             d1 = (xc-xl)*rhx
    627:             d2 = (xc-xr)*rhx
    628:             d3 = (xc-xt)*rhy
    629:             d4 = (xc-xb)*rhy
    630:             d5 = (xrb-xr)*rhy
    631:             d6 = (xrb-xb)*rhx
    632:             d7 = (xlt-xl)*rhy
    633:             d8 = (xlt-xt)*rhx

    635:             f1 = sqrt(1.0 + d1*d1 + d7*d7)
    636:             f2 = sqrt(1.0 + d1*d1 + d4*d4)
    637:             f3 = sqrt(1.0 + d3*d3 + d8*d8)
    638:             f4 = sqrt(1.0 + d3*d3 + d2*d2)
    639:             f5 = sqrt(1.0 + d2*d2 + d5*d5)
    640:             f6 = sqrt(1.0 + d4*d4 + d6*d6)

    643:             hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+                 &
    644:      &              (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)

    646:             hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+                 &
    647:      &            (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)

    649:             ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+                 &
    650:      &                (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)

    652:             hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+                 &
    653:      &              (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)

    655:             hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
    656:             htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)

    658:             hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) +                         &
    659:      &              hxdhy*(1.0+d8*d8)/(f3*f3*f3) +                      &
    660:      &              hydhx*(1.0+d5*d5)/(f5*f5*f5) +                      &
    661:      &              hxdhy*(1.0+d6*d6)/(f6*f6*f6) +                      &
    662:      &              (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-               &
    663:      &              2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+          &
    664:      &              hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)

    666:             hl = hl * 0.5
    667:             hr = hr * 0.5
    668:             ht = ht * 0.5
    669:             hb = hb * 0.5
    670:             hbr = hbr * 0.5
    671:             htl = htl * 0.5
    672:             hc = hc * 0.5

    674:             k = 0

    676:             if (j .gt. 0) then
    677:                v(k) = hb
    678:                col(k) = row - gxm
    679:                k=k+1
    680:             endif

    682:             if ((j .gt. 0) .and. (i .lt. mx-1)) then
    683:                v(k) = hbr
    684:                col(k) = row-gxm+1
    685:                k=k+1
    686:             endif

    688:             if (i .gt. 0) then
    689:                v(k) = hl
    690:                col(k) = row - 1
    691:                k = k+1
    692:             endif

    694:             v(k) = hc
    695:             col(k) = row
    696:             k=k+1

    698:             if (i .lt. mx-1) then
    699:                v(k) = hr
    700:                col(k) = row + 1
    701:                k=k+1
    702:             endif

    704:             if ((i .gt. 0) .and. (j .lt. my-1)) then
    705:                v(k) = htl
    706:                col(k) = row + gxm - 1
    707:                k=k+1
    708:             endif

    710:             if (j .lt. my-1) then
    711:                v(k) = ht
    712:                col(k) = row + gxm
    713:                k=k+1
    714:             endif

    716: ! Set matrix values using local numbering, defined earlier in main routine
    717:             call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES,      &
    718:      &                              ierr)

    722:          enddo
    723:       enddo

    725: ! restore vectors
    726:       call VecRestoreArray(localX,x_v,x_i,ierr)
    727:       call VecRestoreArray(Left,left_v,left_i,ierr)
    728:       call VecRestoreArray(Right,right_v,right_i,ierr)
    729:       call VecRestoreArray(Top,top_v,top_i,ierr)
    730:       call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)

    733: ! Assemble the matrix
    734:       call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
    735:       call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)

    737:       call PetscLogFlops(199.0d0*xm*ym,ierr)

    739:       return
    740:       end

    746: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h

    748: ! ----------------------------------------------------------------------------
    749: !
    750: !/*
    751: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
    752: !
    753: !
    754: !*/

    756:       subroutine MSA_BoundaryConditions(ierr)
    757:       use mymodule
    758:       implicit none

    760:       PetscErrorCode   ierr
    761:       PetscInt         i,j,k,limit,maxits
    762:       PetscInt         xs, xm, gxs, gxm
    763:       PetscInt         ys, ym, gys, gym
    764:       PetscInt         bsize, lsize
    765:       PetscInt         tsize, rsize
    766:       PetscReal      one,two,three,tol
    767:       PetscReal      scl,fnorm,det,xt
    768:       PetscReal      yt,hx,hy,u1,u2,nf1,nf2
    769:       PetscReal      njac11,njac12,njac21,njac22
    770:       PetscReal      b, t, l, r
    771:       PetscReal      boundary_v(0:1)
    772:       PetscOffset      boundary_i
    773:       logical exitloop
    774:       PetscBool flg

    776:       limit=0
    777:       maxits = 5
    778:       tol=1e-10
    779:       b=-0.5
    780:       t= 0.5
    781:       l=-0.5
    782:       r= 0.5
    783:       xt=0
    784:       yt=0
    785:       one=1.0
    786:       two=2.0
    787:       three=3.0

    790:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    791:      &                  PETSC_NULL_INTEGER,ierr)
    792:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
    793:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

    795:       bsize = xm + 2
    796:       lsize = ym + 2
    797:       rsize = ym + 2
    798:       tsize = xm + 2

    801:       call VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
    802:       call VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
    803:       call VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
    804:       call VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)

    806:       hx= (r-l)/(mx+1)
    807:       hy= (t-b)/(my+1)

    809:       do j=0,3

    811:          if (j.eq.0) then
    812:             yt=b
    813:             xt=l+hx*xs
    814:             limit=bsize
    815:             call VecGetArray(Bottom,boundary_v,boundary_i,ierr)

    818:          elseif (j.eq.1) then
    819:             yt=t
    820:             xt=l+hx*xs
    821:             limit=tsize
    822:             call VecGetArray(Top,boundary_v,boundary_i,ierr)

    824:          elseif (j.eq.2) then
    825:             yt=b+hy*ys
    826:             xt=l
    827:             limit=lsize
    828:             call VecGetArray(Left,boundary_v,boundary_i,ierr)

    830:          elseif (j.eq.3) then
    831:             yt=b+hy*ys
    832:             xt=r
    833:             limit=rsize
    834:             call VecGetArray(Right,boundary_v,boundary_i,ierr)
    835:          endif

    838:          do i=0,limit-1

    840:             u1=xt
    841:             u2=-yt
    842:             k = 0
    843:             exitloop = .false.
    844:             do while (k .lt. maxits .and. (.not. exitloop))

    846:                nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
    847:                nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
    848:                fnorm=sqrt(nf1*nf1+nf2*nf2)
    849:                if (fnorm .gt. tol) then
    850:                   njac11=one+u2*u2-u1*u1
    851:                   njac12=two*u1*u2
    852:                   njac21=-two*u1*u2
    853:                   njac22=-one - u1*u1 + u2*u2
    854:                   det = njac11*njac22-njac21*njac12
    855:                   u1 = u1-(njac22*nf1-njac12*nf2)/det
    856:                   u2 = u2-(njac11*nf2-njac21*nf1)/det
    857:                else
    858:                   exitloop = .true.
    859:                endif
    860:                k=k+1
    861:             enddo

    863:             boundary_v(i + boundary_i) = u1*u1-u2*u2
    864:             if ((j .eq. 0) .or. (j .eq. 1)) then
    865:                xt = xt + hx
    866:             else
    867:                yt = yt + hy
    868:             endif

    870:          enddo

    873:          if (j.eq.0) then
    874:             call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
    875:          elseif (j.eq.1) then
    876:             call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
    877:          elseif (j.eq.2) then
    878:             call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
    879:          elseif (j.eq.3) then
    880:             call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
    881:          endif

    883:       enddo

    886: ! Scale the boundary if desired
    887:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
    888:      &                         '-bottom',scl,flg,ierr)
    889:       if (flg .eqv. PETSC_TRUE) then
    890:          call VecScale(Bottom,scl,ierr)
    891:       endif

    893:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
    894:      &                         '-top',scl,flg,ierr)
    895:       if (flg .eqv. PETSC_TRUE) then
    896:          call VecScale(Top,scl,ierr)
    897:       endif

    899:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
    900:      &                         '-right',scl,flg,ierr)
    901:       if (flg .eqv. PETSC_TRUE) then
    902:          call VecScale(Right,scl,ierr)
    903:       endif

    905:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
    906:      &                         '-left',scl,flg,ierr)
    907:       if (flg .eqv. PETSC_TRUE) then
    908:          call VecScale(Left,scl,ierr)
    909:       endif

    912:       return
    913:       end

    915: ! ----------------------------------------------------------------------------
    916: !
    917: !/*
    918: !     MSA_Plate - Calculates an obstacle for surface to stretch over
    919: !
    920: !     Output Parameter:
    921: !.    xl - lower bound vector
    922: !.    xu - upper bound vector
    923: !
    924: !*/

    926:       subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
    927:       use mymodule
    928:       implicit none

    930:       Tao        tao
    931:       Vec              xl,xu
    932:       PetscErrorCode   ierr
    933:       PetscInt         i,j,row
    934:       PetscInt         xs, xm, ys, ym
    935:       PetscReal      lb,ub
    936:       PetscInt         dummy

    938: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    939: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
    940: ! will return an array of doubles referenced by x_array offset by x_index.
    941: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    942: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    943:       PetscReal      xl_v(0:1)
    944:       PetscOffset      xl_i

    946:       lb = PETSC_NINFINITY
    947:       ub = PETSC_INFINITY

    949:       if (bmy .lt. 0) bmy = 0
    950:       if (bmy .gt. my) bmy = my
    951:       if (bmx .lt. 0) bmx = 0
    952:       if (bmx .gt. mx) bmx = mx

    955:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    956:      &             PETSC_NULL_INTEGER,ierr)

    958:       call VecSet(xl,lb,ierr)
    959:       call VecSet(xu,ub,ierr)

    961:       call VecGetArray(xl,xl_v,xl_i,ierr)

    964:       do i=xs,xs+xm-1

    966:          do j=ys,ys+ym-1

    968:             row=(j-ys)*xm + (i-xs)

    970:             if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and.           &
    971:      &          j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
    972:                xl_v(xl_i+row) = bheight

    974:             endif

    976:          enddo
    977:       enddo

    980:       call VecRestoreArray(xl,xl_v,xl_i,ierr)

    982:       return
    983:       end

    989: ! ----------------------------------------------------------------------------
    990: !
    991: !/*
    992: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
    993: !
    994: !     Output Parameter:
    995: !.    X - vector for initial guess
    996: !
    997: !*/

    999:       subroutine MSA_InitialPoint(X, ierr)
    1000:       use mymodule
    1001:       implicit none

    1003:       Vec               X
    1004:       PetscErrorCode    ierr
    1005:       PetscInt          start,i,j
    1006:       PetscInt          row
    1007:       PetscInt          xs,xm,gxs,gxm
    1008:       PetscInt          ys,ym,gys,gym
    1009:       PetscReal       zero, np5

    1011: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    1012: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
    1013: ! will return an array of doubles referenced by x_array offset by x_index.
    1014: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    1015: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    1016:       PetscReal   left_v(0:1),right_v(0:1)
    1017:       PetscReal   bottom_v(0:1),top_v(0:1)
    1018:       PetscReal   x_v(0:1)
    1019:       PetscOffset   left_i, right_i, top_i
    1020:       PetscOffset   bottom_i,x_i
    1021:       PetscBool     flg
    1022:       PetscRandom   rctx

    1024:       zero = 0.0
    1025:       np5 = -0.5

    1027:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,        &
    1028:      &                        '-start', start,flg,ierr)

    1030:       if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then  ! the zero vector is reasonable
    1031:          call VecSet(X,zero,ierr)

    1033:       elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then  ! random start -0.5 < xi < 0.5
    1034:          call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
    1035:          do i=0,start-1
    1036:             call VecSetRandom(X,rctx,ierr)
    1037:          enddo

    1039:          call PetscRandomDestroy(rctx,ierr)
    1040:          call VecShift(X,np5,ierr)

    1042:       else   ! average of boundary conditions

    1044: !        Get Local mesh boundaries
    1045:          call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,             &
    1046:      &                     PETSC_NULL_INTEGER,ierr)
    1047:          call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,    &
    1048:      &                     PETSC_NULL_INTEGER,ierr)

    1052: !        Get pointers to vector data
    1053:          call VecGetArray(Top,top_v,top_i,ierr)
    1054:          call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
    1055:          call VecGetArray(Left,left_v,left_i,ierr)
    1056:          call VecGetArray(Right,right_v,right_i,ierr)

    1058:          call VecGetArray(localX,x_v,x_i,ierr)

    1060: !        Perform local computations
    1061:          do  j=ys,ys+ym-1
    1062:             do i=xs,xs+xm-1
    1063:                row = (j-gys)*gxm  + (i-gxs)
    1064:                x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my        &
    1065:      &             + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) +                  &
    1066:      &              (i+1)*left_v(left_i+j-ys+1)/mx       +                  &
    1067:      &              (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
    1068:             enddo
    1069:          enddo

    1071: !        Restore vectors
    1072:          call VecRestoreArray(localX,x_v,x_i,ierr)

    1074:          call VecRestoreArray(Left,left_v,left_i,ierr)
    1075:          call VecRestoreArray(Top,top_v,top_i,ierr)
    1076:          call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
    1077:          call VecRestoreArray(Right,right_v,right_i,ierr)

    1079:          call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
    1080:          call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)

    1082:       endif

    1084:       return
    1085:       end

    1087: !
    1088: !/*TEST
    1089: !
    1090: !   build:
    1091: !      requires: !complex
    1092: !
    1093: !   test:
    1094: !      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
    1095: !      filter: sort -b
    1096: !      filter_output: sort -b
    1097: !      requires: !single
    1098: !
    1099: !   test:
    1100: !      suffix: 2
    1101: !      nsize: 2
    1102: !      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
    1103: !      filter: sort -b
    1104: !      filter_output: sort -b
    1105: !      requires: !single
    1106: !
    1107: !TEST*/