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

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

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

 64:       external FormFunctionGradient
 65:       external FormHessian
 66:       external MSA_BoundaryConditions
 67:       external MSA_Plate
 68:       external MSA_InitialPoint
 69: ! Initialize Tao

 71:       i1=1
 72:       i3=3
 73:       i7=7

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

 81: ! Specify default dimensions of the problem
 82:       mx = 10
 83:       my = 10
 84:       bheight = 0.1

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

 88:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,     &
 89:      &                        '-mx',mx,flg,ierr)
 90:       CHKERRA(ierr)
 91:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,     &
 92:      &                        '-my',my,flg,ierr)
 93:       CHKERRA(ierr)

 95:       bmx = mx/2
 96:       bmy = my/2

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

108: ! Calculate any derived values from parameters
109:       N = mx*my

111: ! Let Petsc determine the dimensions of the local vectors
112:       Nx = PETSC_DECIDE
113:       NY = PETSC_DECIDE

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

119:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,                    &
120:      &     DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,                              &
121:      &     mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,           &
122:      &     dm,ierr)
123:       CHKERRA(ierr)
124:       call DMSetFromOptions(dm,ierr)
125:       call DMSetUp(dm,ierr)

127: ! Extract global and local vectors from DM; The local vectors are
128: ! used solely as work space for the evaluation of the function,
129: ! gradient, and Hessian.  Duplicate for remaining vectors that are
130: ! the same types.

132:       call DMCreateGlobalVector(dm,x,ierr)
133:       CHKERRA(ierr)
134:       call DMCreateLocalVector(dm,localX,ierr)
135:       CHKERRA(ierr)
136:       call VecDuplicate(localX,localV,ierr)
137:       CHKERRA(ierr)

139: ! Create a matrix data structure to store the Hessian.
140: ! Here we (optionally) also associate the local numbering scheme
141: ! with the matrix so that later we can use local indices for matrix
142: ! assembly

144:       call VecGetLocalSize(x,m,ierr)
145:       CHKERRA(ierr)
146:       call MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER,   &
147:      &     i3,PETSC_NULL_INTEGER,H,ierr)
148:       CHKERRA(ierr)

150:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
151:       CHKERRA(ierr)
152:       call DMGetLocalToGlobalMapping(dm,isltog,ierr)
153:       CHKERRA(ierr)
154:       call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
155:       CHKERRA(ierr)

157: ! The Tao code begins here
158: ! Create TAO solver and set desired solution method.
159: ! This problems uses bounded variables, so the
160: ! method must either be 'tao_tron' or 'tao_blmvm'

162:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
163:       CHKERRA(ierr)
164:       call TaoSetType(tao,TAOBLMVM,ierr)
165:       CHKERRA(ierr)

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

169:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
170:      &     FormFunctionGradient,0,ierr)
171:       CHKERRA(ierr)

173:       call TaoSetHessianRoutine(tao,H,H,FormHessian,                    &
174:      &     0, ierr)
175:       CHKERRA(ierr)

177: ! Set Variable bounds
178:       call MSA_BoundaryConditions(ierr)
179:       CHKERRA(ierr)
180:       call TaoSetVariableBoundsRoutine(tao,MSA_Plate,                   &
181:      &     0,ierr)
182:       CHKERRA(ierr)

184: ! Set the initial solution guess
185:       call MSA_InitialPoint(x, ierr)
186:       CHKERRA(ierr)
187:       call TaoSetInitialVector(tao,x,ierr)
188:       CHKERRA(ierr)

190: ! Check for any tao command line options
191:       call TaoSetFromOptions(tao,ierr)
192:       CHKERRA(ierr)

194: ! Solve the application
195:       call TaoSolve(tao,ierr)
196:       CHKERRA(ierr)

198: ! Free TAO data structures
199:       call TaoDestroy(tao,ierr)
200:       CHKERRA(ierr)

202: ! Free PETSc data structures
203:       call VecDestroy(x,ierr)
204:       call VecDestroy(Top,ierr)
205:       call VecDestroy(Bottom,ierr)
206:       call VecDestroy(Left,ierr)
207:       call VecDestroy(Right,ierr)
208:       call MatDestroy(H,ierr)
209:       call VecDestroy(localX,ierr)
210:       call VecDestroy(localV,ierr)
211:       call DMDestroy(dm,ierr)

213: ! Finalize TAO

215:       call PetscFinalize(ierr)

217:       end

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

235:       subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
236:       use mymodule
237:       implicit none

239: ! Input/output variables

241:       Tao        tao
242:       PetscReal      fcn
243:       Vec              X, G
244:       PetscErrorCode   ierr
245:       PetscInt         dummy

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

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

271:       ft = 0.0
272:       zero = 0.0
273:       hx = 1.0/real(mx + 1)
274:       hy = 1.0/real(my + 1)
275:       hydhx = hy/hx
276:       hxdhy = hx/hy
277:       area = 0.5 * hx * hy
278:       rhx = real(mx) + 1.0
279:       rhy = real(my) + 1.0

281: ! Get local mesh boundaries
282:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
283:      &                  PETSC_NULL_INTEGER,ierr)
284:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,             &
285:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

287: ! Scatter ghost points to local vector
288:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
289:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

291: ! Initialize the vector to zero
292:       call VecSet(localV,zero,ierr)

294: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
295:       call VecGetArray(localX,x_v,x_i,ierr)
296:       call VecGetArray(localV,g_v,g_i,ierr)
297:       call VecGetArray(Top,top_v,top_i,ierr)
298:       call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
299:       call VecGetArray(Left,left_v,left_i,ierr)
300:       call VecGetArray(Right,right_v,right_i,ierr)

302: ! Compute function over the locally owned part of the mesh
303:       do j = ys,ys+ym-1
304:          do i = xs,xs+xm-1
305:             row = (j-gys)*gxm + (i-gxs)
306:             xc = x_v(row+x_i)
307:             xt = xc
308:             xb = xc
309:             xr = xc
310:             xl = xc
311:             xrb = xc
312:             xlt = xc

314:             if (i .eq. 0) then !left side
315:                xl = left_v(j - ys + 1 + left_i)
316:                xlt = left_v(j - ys + 2 + left_i)
317:             else
318:                xl = x_v(row - 1 + x_i)
319:             endif

321:             if (j .eq. 0) then !bottom side
322:                xb = bottom_v(i - xs + 1 + bottom_i)
323:                xrb = bottom_v(i - xs + 2 + bottom_i)
324:             else
325:                xb = x_v(row - gxm + x_i)
326:             endif

328:             if (i + 1 .eq. gxs + gxm) then !right side
329:                xr = right_v(j - ys + 1 + right_i)
330:                xrb = right_v(j - ys + right_i)
331:             else
332:                xr = x_v(row + 1 + x_i)
333:             endif

335:             if (j + 1 .eq. gys + gym) then !top side
336:                xt = top_v(i - xs + 1 + top_i)
337:                xlt = top_v(i - xs + top_i)
338:             else
339:                xt = x_v(row + gxm + x_i)
340:             endif

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

346:             if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
347:                xrb = x_v(row + 1 - gxm + x_i)
348:             endif

350:             d1 = xc-xl
351:             d2 = xc-xr
352:             d3 = xc-xt
353:             d4 = xc-xb
354:             d5 = xr-xrb
355:             d6 = xrb-xb
356:             d7 = xlt-xl
357:             d8 = xt-xlt

359:             df1dxc = d1 * hydhx
360:             df2dxc = d1 * hydhx + d4 * hxdhy
361:             df3dxc = d3 * hxdhy
362:             df4dxc = d2 * hydhx + d3 * hxdhy
363:             df5dxc = d2 * hydhx
364:             df6dxc = d4 * hxdhy

366:             d1 = d1 * rhx
367:             d2 = d2 * rhx
368:             d3 = d3 * rhy
369:             d4 = d4 * rhy
370:             d5 = d5 * rhy
371:             d6 = d6 * rhx
372:             d7 = d7 * rhy
373:             d8 = d8 * rhx

375:             f1 = sqrt(1.0 + d1*d1 + d7*d7)
376:             f2 = sqrt(1.0 + d1*d1 + d4*d4)
377:             f3 = sqrt(1.0 + d3*d3 + d8*d8)
378:             f4 = sqrt(1.0 + d3*d3 + d2*d2)
379:             f5 = sqrt(1.0 + d2*d2 + d5*d5)
380:             f6 = sqrt(1.0 + d4*d4 + d6*d6)

382:             ft = ft + f2 + f4

384:             df1dxc = df1dxc / f1
385:             df2dxc = df2dxc / f2
386:             df3dxc = df3dxc / f3
387:             df4dxc = df4dxc / f4
388:             df5dxc = df5dxc / f5
389:             df6dxc = df6dxc / f6

391:             g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc +  &
392:      &                              df5dxc + df6dxc)
393:          enddo
394:       enddo

396: ! Compute triangular areas along the border of the domain.
397:       if (xs .eq. 0) then  ! left side
398:          do j=ys,ys+ym-1
399:             d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i))         &
400:      &                 * rhy
401:             d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i))        &
402:      &                 * rhx
403:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
404:          enddo
405:       endif

407:       if (ys .eq. 0) then !bottom side
408:          do i=xs,xs+xm-1
409:             d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i))    &
410:      &                    * rhx
411:             d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
412:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
413:          enddo
414:       endif

416:       if (xs + xm .eq. mx) then ! right side
417:          do j=ys,ys+ym-1
418:             d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
419:             d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
420:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
421:          enddo
422:       endif

424:       if (ys + ym .eq. my) then
425:          do i=xs,xs+xm-1
426:             d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
427:             d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
428:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
429:          enddo
430:       endif

432:       if ((ys .eq. 0) .and. (xs .eq. 0)) then
433:          d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
434:          d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
435:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
436:       endif

438:       if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
439:          d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
440:          d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
441:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
442:       endif

444:       ft = ft * area
445:       call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR,                            &
446:      &             MPIU_SUM,PETSC_COMM_WORLD,ierr)

448: ! Restore vectors
449:       call VecRestoreArray(localX,x_v,x_i,ierr)
450:       call VecRestoreArray(localV,g_v,g_i,ierr)
451:       call VecRestoreArray(Left,left_v,left_i,ierr)
452:       call VecRestoreArray(Top,top_v,top_i,ierr)
453:       call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
454:       call VecRestoreArray(Right,right_v,right_i,ierr)

456: ! Scatter values to global vector
457:       call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
458:       call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)

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

462:       return
463:       end  !FormFunctionGradient

465: ! ----------------------------------------------------------------------------
466: !
467: !
468: !   FormHessian - Evaluates Hessian matrix.
469: !
470: !   Input Parameters:
471: !.  tao  - the Tao context
472: !.  X    - input vector
473: !.  dummy  - not used
474: !
475: !   Output Parameters:
476: !.  Hessian    - Hessian matrix
477: !.  Hpc    - optionally different preconditioning matrix
478: !.  flag - flag indicating matrix structure
479: !
480: !   Notes:
481: !   Due to mesh point reordering with DMs, we must always work
482: !   with the local mesh points, and then transform them to the new
483: !   global numbering with the local-to-global mapping.  We cannot work
484: !   directly with the global numbers for the original uniprocessor mesh!
485: !
486: !      MatSetValuesLocal(), using the local ordering (including
487: !         ghost points!)
488: !         - Then set matrix entries using the local ordering
489: !           by calling MatSetValuesLocal()

491:       subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
492:       use mymodule
493:       implicit none

495:       Tao     tao
496:       Vec            X
497:       Mat            Hessian,Hpc
498:       PetscInt       dummy
499:       PetscErrorCode ierr

501:       PetscInt       i,j,k,row
502:       PetscInt       xs,xm,gxs,gxm
503:       PetscInt       ys,ym,gys,gym
504:       PetscInt       col(0:6)
505:       PetscReal    hx,hy,hydhx,hxdhy,rhx,rhy
506:       PetscReal    f1,f2,f3,f4,f5,f6,d1,d2,d3
507:       PetscReal    d4,d5,d6,d7,d8
508:       PetscReal    xc,xl,xr,xt,xb,xlt,xrb
509:       PetscReal    hl,hr,ht,hb,hc,htl,hbr

511: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
512: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
513: ! will return an array of doubles referenced by x_array offset by x_index.
514: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
515: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
516:       PetscReal   right_v(0:1),left_v(0:1)
517:       PetscReal   bottom_v(0:1),top_v(0:1)
518:       PetscReal   x_v(0:1)
519:       PetscOffset   x_i,right_i,left_i
520:       PetscOffset   bottom_i,top_i
521:       PetscReal   v(0:6)
522:       PetscBool     assembled
523:       PetscInt      i1

525:       i1=1

527: ! Set various matrix options
528:       call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,              &
529:      &                  PETSC_TRUE,ierr)

531: ! Get local mesh boundaries
532:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
533:      &                  PETSC_NULL_INTEGER,ierr)
534:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
535:      &                       PETSC_NULL_INTEGER,ierr)

537: ! Scatter ghost points to local vectors
538:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
539:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

541: ! Get pointers to vector data (see note on Fortran arrays above)
542:       call VecGetArray(localX,x_v,x_i,ierr)
543:       call VecGetArray(Top,top_v,top_i,ierr)
544:       call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
545:       call VecGetArray(Left,left_v,left_i,ierr)
546:       call VecGetArray(Right,right_v,right_i,ierr)

548: ! Initialize matrix entries to zero
549:       call MatAssembled(Hessian,assembled,ierr)
550:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)

552:       rhx = real(mx) + 1.0
553:       rhy = real(my) + 1.0
554:       hx = 1.0/rhx
555:       hy = 1.0/rhy
556:       hydhx = hy/hx
557:       hxdhy = hx/hy
558: ! compute Hessian over the locally owned part of the mesh

560:       do  i=xs,xs+xm-1
561:          do  j=ys,ys+ym-1
562:             row = (j-gys)*gxm + (i-gxs)

564:             xc = x_v(row + x_i)
565:             xt = xc
566:             xb = xc
567:             xr = xc
568:             xl = xc
569:             xrb = xc
570:             xlt = xc

572:             if (i .eq. gxs) then   ! Left side
573:                xl = left_v(left_i + j - ys + 1)
574:                xlt = left_v(left_i + j - ys + 2)
575:             else
576:                xl = x_v(x_i + row -1)
577:             endif

579:             if (j .eq. gys) then ! bottom side
580:                xb = bottom_v(bottom_i + i - xs + 1)
581:                xrb = bottom_v(bottom_i + i - xs + 2)
582:             else
583:                xb = x_v(x_i + row - gxm)
584:             endif

586:             if (i+1 .eq. gxs + gxm) then !right side
587:                xr = right_v(right_i + j - ys + 1)
588:                xrb = right_v(right_i + j - ys)
589:             else
590:                xr = x_v(x_i + row + 1)
591:             endif

593:             if (j+1 .eq. gym+gys) then !top side
594:                xt = top_v(top_i +i - xs + 1)
595:                xlt = top_v(top_i + i - xs)
596:             else
597:                xt = x_v(x_i + row + gxm)
598:             endif

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

604:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
605:                xrb = x_v(x_i + row + 1 - gxm)
606:             endif

608:             d1 = (xc-xl)*rhx
609:             d2 = (xc-xr)*rhx
610:             d3 = (xc-xt)*rhy
611:             d4 = (xc-xb)*rhy
612:             d5 = (xrb-xr)*rhy
613:             d6 = (xrb-xb)*rhx
614:             d7 = (xlt-xl)*rhy
615:             d8 = (xlt-xt)*rhx

617:             f1 = sqrt(1.0 + d1*d1 + d7*d7)
618:             f2 = sqrt(1.0 + d1*d1 + d4*d4)
619:             f3 = sqrt(1.0 + d3*d3 + d8*d8)
620:             f4 = sqrt(1.0 + d3*d3 + d2*d2)
621:             f5 = sqrt(1.0 + d2*d2 + d5*d5)
622:             f6 = sqrt(1.0 + d4*d4 + d6*d6)

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

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

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

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

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

639:             hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) +                         &
640:      &              hxdhy*(1.0+d8*d8)/(f3*f3*f3) +                      &
641:      &              hydhx*(1.0+d5*d5)/(f5*f5*f5) +                      &
642:      &              hxdhy*(1.0+d6*d6)/(f6*f6*f6) +                      &
643:      &              (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-               &
644:      &              2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+          &
645:      &              hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)

647:             hl = hl * 0.5
648:             hr = hr * 0.5
649:             ht = ht * 0.5
650:             hb = hb * 0.5
651:             hbr = hbr * 0.5
652:             htl = htl * 0.5
653:             hc = hc * 0.5

655:             k = 0

657:             if (j .gt. 0) then
658:                v(k) = hb
659:                col(k) = row - gxm
660:                k=k+1
661:             endif

663:             if ((j .gt. 0) .and. (i .lt. mx-1)) then
664:                v(k) = hbr
665:                col(k) = row-gxm+1
666:                k=k+1
667:             endif

669:             if (i .gt. 0) then
670:                v(k) = hl
671:                col(k) = row - 1
672:                k = k+1
673:             endif

675:             v(k) = hc
676:             col(k) = row
677:             k=k+1

679:             if (i .lt. mx-1) then
680:                v(k) = hr
681:                col(k) = row + 1
682:                k=k+1
683:             endif

685:             if ((i .gt. 0) .and. (j .lt. my-1)) then
686:                v(k) = htl
687:                col(k) = row + gxm - 1
688:                k=k+1
689:             endif

691:             if (j .lt. my-1) then
692:                v(k) = ht
693:                col(k) = row + gxm
694:                k=k+1
695:             endif

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

701:          enddo
702:       enddo

704: ! restore vectors
705:       call VecRestoreArray(localX,x_v,x_i,ierr)
706:       call VecRestoreArray(Left,left_v,left_i,ierr)
707:       call VecRestoreArray(Right,right_v,right_i,ierr)
708:       call VecRestoreArray(Top,top_v,top_i,ierr)
709:       call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)

711: ! Assemble the matrix
712:       call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
713:       call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)

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

717:       return
718:       end

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

722: ! ----------------------------------------------------------------------------
723: !
724: !/*
725: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
726: !
727: !
728: !*/

730:       subroutine MSA_BoundaryConditions(ierr)
731:       use mymodule
732:       implicit none

734:       PetscErrorCode   ierr
735:       PetscInt         i,j,k,limit,maxits
736:       PetscInt         xs, xm, gxs, gxm
737:       PetscInt         ys, ym, gys, gym
738:       PetscInt         bsize, lsize
739:       PetscInt         tsize, rsize
740:       PetscReal      one,two,three,tol
741:       PetscReal      scl,fnorm,det,xt
742:       PetscReal      yt,hx,hy,u1,u2,nf1,nf2
743:       PetscReal      njac11,njac12,njac21,njac22
744:       PetscReal      b, t, l, r
745:       PetscReal      boundary_v(0:1)
746:       PetscOffset      boundary_i
747:       logical exitloop
748:       PetscBool flg

750:       limit=0
751:       maxits = 5
752:       tol=1e-10
753:       b=-0.5
754:       t= 0.5
755:       l=-0.5
756:       r= 0.5
757:       xt=0
758:       yt=0
759:       one=1.0
760:       two=2.0
761:       three=3.0

763:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
764:      &                  PETSC_NULL_INTEGER,ierr)
765:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
766:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

768:       bsize = xm + 2
769:       lsize = ym + 2
770:       rsize = ym + 2
771:       tsize = xm + 2

773:       call VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
774:       call VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
775:       call VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
776:       call VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)

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

781:       do j=0,3

783:          if (j.eq.0) then
784:             yt=b
785:             xt=l+hx*xs
786:             limit=bsize
787:             call VecGetArray(Bottom,boundary_v,boundary_i,ierr)

789:          elseif (j.eq.1) then
790:             yt=t
791:             xt=l+hx*xs
792:             limit=tsize
793:             call VecGetArray(Top,boundary_v,boundary_i,ierr)

795:          elseif (j.eq.2) then
796:             yt=b+hy*ys
797:             xt=l
798:             limit=lsize
799:             call VecGetArray(Left,boundary_v,boundary_i,ierr)

801:          elseif (j.eq.3) then
802:             yt=b+hy*ys
803:             xt=r
804:             limit=rsize
805:             call VecGetArray(Right,boundary_v,boundary_i,ierr)
806:          endif

808:          do i=0,limit-1

810:             u1=xt
811:             u2=-yt
812:             k = 0
813:             exitloop = .false.
814:             do while (k .lt. maxits .and. (.not. exitloop))

816:                nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
817:                nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
818:                fnorm=sqrt(nf1*nf1+nf2*nf2)
819:                if (fnorm .gt. tol) then
820:                   njac11=one+u2*u2-u1*u1
821:                   njac12=two*u1*u2
822:                   njac21=-two*u1*u2
823:                   njac22=-one - u1*u1 + u2*u2
824:                   det = njac11*njac22-njac21*njac12
825:                   u1 = u1-(njac22*nf1-njac12*nf2)/det
826:                   u2 = u2-(njac11*nf2-njac21*nf1)/det
827:                else
828:                   exitloop = .true.
829:                endif
830:                k=k+1
831:             enddo

833:             boundary_v(i + boundary_i) = u1*u1-u2*u2
834:             if ((j .eq. 0) .or. (j .eq. 1)) then
835:                xt = xt + hx
836:             else
837:                yt = yt + hy
838:             endif

840:          enddo

842:          if (j.eq.0) then
843:             call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
844:          elseif (j.eq.1) then
845:             call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
846:          elseif (j.eq.2) then
847:             call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
848:          elseif (j.eq.3) then
849:             call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
850:          endif

852:       enddo

854: ! Scale the boundary if desired
855:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
856:      &                         '-bottom',scl,flg,ierr)
857:       if (flg .eqv. PETSC_TRUE) then
858:          call VecScale(Bottom,scl,ierr)
859:       endif

861:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
862:      &                         '-top',scl,flg,ierr)
863:       if (flg .eqv. PETSC_TRUE) then
864:          call VecScale(Top,scl,ierr)
865:       endif

867:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
868:      &                         '-right',scl,flg,ierr)
869:       if (flg .eqv. PETSC_TRUE) then
870:          call VecScale(Right,scl,ierr)
871:       endif

873:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,          &
874:      &                         '-left',scl,flg,ierr)
875:       if (flg .eqv. PETSC_TRUE) then
876:          call VecScale(Left,scl,ierr)
877:       endif

879:       return
880:       end

882: ! ----------------------------------------------------------------------------
883: !
884: !/*
885: !     MSA_Plate - Calculates an obstacle for surface to stretch over
886: !
887: !     Output Parameter:
888: !.    xl - lower bound vector
889: !.    xu - upper bound vector
890: !
891: !*/

893:       subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
894:       use mymodule
895:       implicit none

897:       Tao        tao
898:       Vec              xl,xu
899:       PetscErrorCode   ierr
900:       PetscInt         i,j,row
901:       PetscInt         xs, xm, ys, ym
902:       PetscReal      lb,ub
903:       PetscInt         dummy

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

913:       lb = PETSC_NINFINITY
914:       ub = PETSC_INFINITY

916:       if (bmy .lt. 0) bmy = 0
917:       if (bmy .gt. my) bmy = my
918:       if (bmx .lt. 0) bmx = 0
919:       if (bmx .gt. mx) bmx = mx

921:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
922:      &             PETSC_NULL_INTEGER,ierr)

924:       call VecSet(xl,lb,ierr)
925:       call VecSet(xu,ub,ierr)

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

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

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

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

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

939:             endif

941:          enddo
942:       enddo

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

946:       return
947:       end

949: ! ----------------------------------------------------------------------------
950: !
951: !/*
952: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
953: !
954: !     Output Parameter:
955: !.    X - vector for initial guess
956: !
957: !*/

959:       subroutine MSA_InitialPoint(X, ierr)
960:       use mymodule
961:       implicit none

963:       Vec               X
964:       PetscErrorCode    ierr
965:       PetscInt          start,i,j
966:       PetscInt          row
967:       PetscInt          xs,xm,gxs,gxm
968:       PetscInt          ys,ym,gys,gym
969:       PetscReal       zero, np5

971: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
972: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
973: ! will return an array of doubles referenced by x_array offset by x_index.
974: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
975: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
976:       PetscReal   left_v(0:1),right_v(0:1)
977:       PetscReal   bottom_v(0:1),top_v(0:1)
978:       PetscReal   x_v(0:1)
979:       PetscOffset   left_i, right_i, top_i
980:       PetscOffset   bottom_i,x_i
981:       PetscBool     flg
982:       PetscRandom   rctx

984:       zero = 0.0
985:       np5 = -0.5

987:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,        &
988:      &                        '-start', start,flg,ierr)

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

993:       elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then  ! random start -0.5 < xi < 0.5
994:          call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
995:          do i=0,start-1
996:             call VecSetRandom(X,rctx,ierr)
997:          enddo

999:          call PetscRandomDestroy(rctx,ierr)
1000:          call VecShift(X,np5,ierr)

1002:       else   ! average of boundary conditions

1004: !        Get Local mesh boundaries
1005:          call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,             &
1006:      &                     PETSC_NULL_INTEGER,ierr)
1007:          call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,    &
1008:      &                     PETSC_NULL_INTEGER,ierr)

1010: !        Get pointers to vector data
1011:          call VecGetArray(Top,top_v,top_i,ierr)
1012:          call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1013:          call VecGetArray(Left,left_v,left_i,ierr)
1014:          call VecGetArray(Right,right_v,right_i,ierr)

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

1018: !        Perform local computations
1019:          do  j=ys,ys+ym-1
1020:             do i=xs,xs+xm-1
1021:                row = (j-gys)*gxm  + (i-gxs)
1022:                x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my        &
1023:      &             + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) +                  &
1024:      &              (i+1)*left_v(left_i+j-ys+1)/mx       +                  &
1025:      &              (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1026:             enddo
1027:          enddo

1029: !        Restore vectors
1030:          call VecRestoreArray(localX,x_v,x_i,ierr)

1032:          call VecRestoreArray(Left,left_v,left_i,ierr)
1033:          call VecRestoreArray(Top,top_v,top_i,ierr)
1034:          call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1035:          call VecRestoreArray(Right,right_v,right_i,ierr)

1037:          call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1038:          call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)

1040:       endif

1042:       return
1043:       end

1045: !
1046: !/*TEST
1047: !
1048: !   build:
1049: !      requires: !complex
1050: !
1051: !   test:
1052: !      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1053: !      filter: sort -b
1054: !      filter_output: sort -b
1055: !      requires: !single
1056: !
1057: !   test:
1058: !      suffix: 2
1059: !      nsize: 2
1060: !      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1061: !      filter: sort -b
1062: !      filter_output: sort -b
1063: !      requires: !single
1064: !
1065: !TEST*/