Actual source code: plate2f.F

tao-2.1-p0 2012-07-24
```  1: !  Program usage: mpirun -np <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: TaoInitialize(); TaoFinalize();
18: !   Routines: TaoCreate();
19: !   Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
20: !   Routines: TaoSetHessianRoutine();
21: !   Routines: TaoSetVariableBoundsRoutine();
22: !   Routines: TaoSetInitialVector();
23: !   Routines: TaoSetFromOptions();
24: !   Routines: TaoSolve();
25: !   Routines: TaoDestroy();
26: !   Processors: n
27: !T*/

31:       implicit none

33: #include "plate2f.h"

35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: !                   Variable declarations
37: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: !
39: !  Variables:
40: !    (common from plate2f.h):
41: !    Nx, Ny           number of processors in x- and y- directions
42: !    mx, my           number of grid points in x,y directions
43: !    N    global dimension of vector

45:       PetscErrorCode   ierr          ! used to check for functions returning nonzeros
46:       Vec              x             ! solution vector
47:       Vec              xl, xu        ! lower and upper bounds vectorsp
48:       PetscInt         m             ! number of local elements in vector
49:       TaoSolver        tao           ! TaoSolver solver context
50:       Mat              H             ! Hessian matrix
51:       ISLocalToGlobalMapping isltog  ! local to global mapping object
52:       PetscBool        flg
53:       PetscInt         i1,i3,i7

56:       external FormFunctionGradient
57:       external FormHessian
58:       external MSA_BoundaryConditions
59:       external MSA_Plate
60:       external MSA_InitialPoint
61: ! Initialize Tao

63:       i1=1
64:       i3=3
65:       i7=7

68:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
69:       CHKERRQ(ierr)
70:       call TaoInitialize(PETSC_NULL_CHARACTER,ierr)
71:

73: ! Specify default dimensions of the problem
74:       mx = 10
75:       my = 10
76:       bheight = 0.1

78: ! Check for any command line arguments that override defaults
79:
80:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mx",mx,flg,ierr)
81:       CHKERRQ(ierr)
82:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,ierr)
83:       CHKERRQ(ierr)
84:
85:       bmx = mx/2
86:       bmy = my/2

88:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-bmx",bmx,flg,ierr)
89:       CHKERRQ(ierr)
90:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-bmy",bmy,flg,ierr)
91:       CHKERRQ(ierr)
92:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bheight",bheight,   &
93:      &      flg,ierr)
94:       CHKERRQ(ierr)
95:

97: ! Calculate any derived values from parameters
98:       N = mx*my

100: ! Let Petsc determine the dimensions of the local vectors
101:       Nx = PETSC_DECIDE
102:       NY = PETSC_DECIDE

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

108:       call DMDACreate2d(MPI_COMM_WORLD,DMDA_BOUNDARY_NONE,                    &
109:      &     DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX,                              &
110:      &     mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,           &
111:      &     dm,ierr)
112:       CHKERRQ(ierr)
113:

115: ! Extract global and local vectors from DM; The local vectors are
116: ! used solely as work space for the evaluation of the function,
117: ! gradient, and Hessian.  Duplicate for remaining vectors that are
118: ! the same types.

120:       call DMCreateGlobalVector(dm,x,ierr)
121:       CHKERRQ(ierr)
122:       call DMCreateLocalVector(dm,localX,ierr)
123:       CHKERRQ(ierr)
124:       call VecDuplicate(localX,localV,ierr)
125:       CHKERRQ(ierr)

127: ! Create a matrix data structure to store the Hessian.
128: ! Here we (optionally) also associate the local numbering scheme
129: ! with the matrix so that later we can use local indices for matrix
130: ! assembly

132:       call VecGetLocalSize(x,m,ierr)
133:       CHKERRQ(ierr)
134:       call MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER,   &
135:      &     i3,PETSC_NULL_INTEGER,H,ierr)
136:       CHKERRQ(ierr)

138:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
139:       CHKERRQ(ierr)
140:       call DMGetLocalToGlobalMapping(dm,isltog,ierr)
141:       CHKERRQ(ierr)
142:       call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
143:       CHKERRQ(ierr)
144:

146: ! The Tao code begins here
147: ! Create TAO solver and set desired solution method.
148: ! This problems uses bounded variables, so the
149: ! method must either be 'tao_tron' or 'tao_blmvm'

151:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
152:       CHKERRQ(ierr)
153:       call TaoSetType(tao,'tao_blmvm',ierr)
154:       CHKERRQ(ierr)

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

158:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
159:      &     FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
160:       CHKERRQ(ierr)

162:       call TaoSetHessianRoutine(tao,H,H,FormHessian,                    &
163:      &     PETSC_NULL_OBJECT, ierr)
164:       CHKERRQ(ierr)

166: ! Set Variable bounds
167:       call MSA_BoundaryConditions(ierr)
168:       CHKERRQ(ierr)
169:       call TaoSetVariableBoundsRoutine(tao,MSA_Plate,                   &
170:      &     PETSC_NULL_OBJECT,ierr)
171:       CHKERRQ(ierr)

173: ! Set the initial solution guess
174:       call MSA_InitialPoint(x, ierr)
175:       CHKERRQ(ierr)
176:       call TaoSetInitialVector(tao,x,ierr)
177:       CHKERRQ(ierr)

179: ! Check for any tao command line options
180:       call TaoSetFromOptions(tao,ierr)
181:       CHKERRQ(ierr)

183: ! Solve the application
184:       call TaoSolve(tao,ierr)
185:       CHKERRQ(ierr)

187: ! Free TAO data structures
188:       call TaoDestroy(tao,ierr)
189:       CHKERRQ(ierr)

191: ! Free PETSc data structures
192:       call VecDestroy(x,ierr)
193:       call VecDestroy(Top,ierr)
194:       call VecDestroy(Bottom,ierr)
195:       call VecDestroy(Left,ierr)
196:       call VecDestroy(Right,ierr)
197:       call MatDestroy(H,ierr)
198:       call VecDestroy(localX,ierr)
199:       call VecDestroy(localV,ierr)
200:       call DMDestroy(dm,ierr)

202: ! Finalize TAO

204:       call TaoFinalize(ierr)
205:       call PetscFinalize(ierr)

207:       end

209: ! ---------------------------------------------------------------------
210: !
211: !  FormFunctionGradient - Evaluates function f(X).
212: !
213: !  Input Parameters:
214: !  tao   - the TaoSolver context
215: !  X     - the input vector
216: !  dummy - optional user-defined context, as set by TaoSetFunction()
217: !          (not used here)
218: !
219: !  Output Parameters:
220: !  fcn     - the newly evaluated function
221: !  G       - the gradient vector
222: !  info  - error code
223: !

226:       subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
227:       implicit none

229: ! dm, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
230: #include "plate2f.h"
231:
232: ! Input/output variables

234:       TaoSolver        tao
235:       PetscReal      fcn
236:       Vec              X, G
237:       PetscErrorCode   ierr
238:       PetscInt         dummy
239:
240:       PetscInt         i,j,row
241:       PetscInt         xs, xm
242:       PetscInt         gxs, gxm
243:       PetscInt         ys, ym
244:       PetscInt         gys, gym
245:       PetscReal      ft,zero,hx,hy,hydhx,hxdhy
246:       PetscReal      area,rhx,rhy
247:       PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3
248:       PetscReal      d4,d5,d6,d7,d8
249:       PetscReal      df1dxc,df2dxc,df3dxc,df4dxc
250:       PetscReal      df5dxc,df6dxc
251:       PetscReal      xc,xl,xr,xt,xb,xlt,xrb

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

265:       ft = 0.0d0
266:       zero = 0.0d0
267:       hx = 1.0d0/(mx + 1)
268:       hy = 1.0d0/(my + 1)
269:       hydhx = hy/hx
270:       hxdhy = hx/hy
271:       area = 0.5d0 * hx * hy
272:       rhx = mx + 1.0d0
273:       rhy = my + 1.0d0

276: ! Get local mesh boundaries
277:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
278:      &                  PETSC_NULL_INTEGER,ierr)
279:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,             &
280:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

282: ! Scatter ghost points to local vector
283:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
284:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

286: ! Initialize the vector to zero
287:       call VecSet(localV,zero,ierr)

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

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

309:             if (i .eq. 0) then !left side
310:                xl = left_v(j - ys + 1 + left_i)
311:                xlt = left_v(j - ys + 2 + left_i)
312:             else
313:                xl = x_v(row - 1 + x_i)
314:             endif

316:             if (j .eq. 0) then !bottom side
317:                xb = bottom_v(i - xs + 1 + bottom_i)
318:                xrb = bottom_v(i - xs + 2 + bottom_i)
319:             else
320:                xb = x_v(row - gxm + x_i)
321:             endif

323:             if (i + 1 .eq. gxs + gxm) then !right side
324:                xr = right_v(j - ys + 1 + right_i)
325:                xrb = right_v(j - ys + right_i)
326:             else
327:                xr = x_v(row + 1 + x_i)
328:             endif

330:             if (j + 1 .eq. gys + gym) then !top side
331:                xt = top_v(i - xs + 1 + top_i)
332:                xlt = top_v(i - xs + top_i)
333:             else
334:                xt = x_v(row + gxm + x_i)
335:             endif

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

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

345:             d1 = xc-xl
346:             d2 = xc-xr
347:             d3 = xc-xt
348:             d4 = xc-xb
349:             d5 = xr-xrb
350:             d6 = xrb-xb
351:             d7 = xlt-xl
352:             d8 = xt-xlt

354:             df1dxc = d1 * hydhx
355:             df2dxc = d1 * hydhx + d4 * hxdhy
356:             df3dxc = d3 * hxdhy
357:             df4dxc = d2 * hydhx + d3 * hxdhy
358:             df5dxc = d2 * hydhx
359:             df6dxc = d4 * hxdhy

361:             d1 = d1 * rhx
362:             d2 = d2 * rhx
363:             d3 = d3 * rhy
364:             d4 = d4 * rhy
365:             d5 = d5 * rhy
366:             d6 = d6 * rhx
367:             d7 = d7 * rhy
368:             d8 = d8 * rhx

370:             f1 = sqrt(1.0d0 + d1*d1 + d7*d7)
371:             f2 = sqrt(1.0d0 + d1*d1 + d4*d4)
372:             f3 = sqrt(1.0d0 + d3*d3 + d8*d8)
373:             f4 = sqrt(1.0d0 + d3*d3 + d2*d2)
374:             f5 = sqrt(1.0d0 + d2*d2 + d5*d5)
375:             f6 = sqrt(1.0d0 + d4*d4 + d6*d6)

377:             ft = ft + f2 + f4

379:             df1dxc = df1dxc / f1
380:             df2dxc = df2dxc / f2
381:             df3dxc = df3dxc / f3
382:             df4dxc = df4dxc / f4
383:             df5dxc = df5dxc / f5
384:             df6dxc = df6dxc / f6

386:             g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc +  &
387:      &                              df5dxc + df6dxc)
388:          enddo
389:       enddo

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

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

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

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

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

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

443:       ft = ft * area
444:       call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR,                            &
445:      &             MPI_SUM,MPI_COMM_WORLD,ierr)

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

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

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

463:       return
464:       end  !FormFunctionGradient
465:

470: ! ----------------------------------------------------------------------------
471: !
472: !/*
473: !   FormHessian - Evaluates Hessian matrix.
474: !
475: !   Input Parameters:
476: !.  tao  - the TaoSolver context
477: !.  X    - input vector
478: !.  dummy  - not used
479: !
480: !   Output Parameters:
481: !.  Hessian    - Hessian matrix
482: !.  Hpc    - optionally different preconditioning matrix
483: !.  flag - flag indicating matrix structure
484: !
485: !   Notes:
486: !   Due to mesh point reordering with DMs, we must always work
487: !   with the local mesh points, and then transform them to the new
488: !   global numbering with the local-to-global mapping.  We cannot work
489: !   directly with the global numbers for the original uniprocessor mesh!
490: !
491: !   Two methods are available for imposing this transformation
492: !   when setting matrix entries:
493: !     (A) MatSetValuesLocal(), using the local ordering (including
494: !         ghost points!)
495: !         - Do the following two steps once, before calling TaoSolve()
496: !           - Use DMDAGetISLocalToGlobalMapping() to extract the
497: !             local-to-global map from the DM
498: !           - Associate this map with the matrix by calling
499: !             MatSetLocalToGlobalMapping()
500: !         - Then set matrix entries using the local ordering
501: !           by calling MatSetValuesLocal()
502: !     (B) MatSetValues(), using the global ordering
503: !         - Use DMDAGetGlobalIndices() to extract the local-to-global map
504: !         - Then apply this map explicitly yourself
505: !         - Set matrix entries using the global ordering by calling
506: !           MatSetValues()
507: !   Option (A) seems cleaner/easier in many cases, and is the procedure
508: !   used in this example.
509: */
510:       subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, ierr)
511:       implicit none

513: ! dm,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
514: #include "plate2f.h"
515:
516:       TaoSolver     tao
517:       Vec            X
518:       Mat            Hessian,Hpc
519:       MatStructure   flg
520:       PetscInt       dummy
521:       PetscErrorCode ierr

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

533: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
534: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
535: ! will return an array of doubles referenced by x_array offset by x_index.
536: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
537: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
538:       PetscReal   right_v(0:1),left_v(0:1)
539:       PetscReal   bottom_v(0:1),top_v(0:1)
540:       PetscReal   x_v(0:1)
541:       PetscOffset   x_i,right_i,left_i
542:       PetscOffset   bottom_i,top_i
543:       PetscReal   v(0:6)
544:       PetscBool     assembled
545:       PetscInt      i1
546:
547:       i1=1
548:
549: ! Set various matrix options
550:       call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,              &
551:      &                  PETSC_TRUE,ierr)

553: ! Get local mesh boundaries
554:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
555:      &                  PETSC_NULL_INTEGER,ierr)
556:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
557:      &                       PETSC_NULL_INTEGER,ierr)

559: ! Scatter ghost points to local vectors
560:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
561:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

563: ! Get pointers to vector data (see note on Fortran arrays above)
564:       call VecGetArray(localX,x_v,x_i,ierr)
565:       call VecGetArray(Top,top_v,top_i,ierr)
566:       call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
567:       call VecGetArray(Left,left_v,left_i,ierr)
568:       call VecGetArray(Right,right_v,right_i,ierr)

570: ! Initialize matrix entries to zero
571:       call MatAssembled(Hessian,assembled,ierr)
572:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)

575:       rhx = mx + 1.0
576:       rhy = my + 1.0
577:       hx = 1.0/rhx
578:       hy = 1.0/rhy
579:       hydhx = hy/hx
580:       hxdhy = hx/hy
581: ! compute Hessian over the locally owned part of the mesh

583:       do  i=xs,xs+xm-1
584:          do  j=ys,ys+ym-1
585:             row = (j-gys)*gxm + (i-gxs)
586:
587:             xc = x_v(row + x_i)
588:             xt = xc
589:             xb = xc
590:             xr = xc
591:             xl = xc
592:             xrb = xc
593:             xlt = xc

595:             if (i .eq. gxs) then   ! Left side
596:                xl = left_v(left_i + j - ys + 1)
597:                xlt = left_v(left_i + j - ys + 2)
598:             else
599:                xl = x_v(x_i + row -1 )
600:             endif

602:             if (j .eq. gys) then ! bottom side
603:                xb = bottom_v(bottom_i + i - xs + 1)
604:                xrb = bottom_v(bottom_i + i - xs + 2)
605:             else
606:                xb = x_v(x_i + row - gxm)
607:             endif

609:             if (i+1 .eq. gxs + gxm) then !right side
610:                xr = right_v(right_i + j - ys + 1)
611:                xrb = right_v(right_i + j - ys)
612:             else
613:                xr = x_v(x_i + row + 1)
614:             endif

616:             if (j+1 .eq. gym+gys) then !top side
617:                xt = top_v(top_i +i - xs + 1)
618:                xlt = top_v(top_i + i - xs)
619:             else
620:                xt = x_v(x_i + row + gxm)
621:             endif

623:             if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
624:                xlt = x_v(x_i + row - 1 + gxm)
625:             endif
626:
627:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
628:                xrb = x_v(x_i + row + 1 - gxm)
629:             endif

631:             d1 = (xc-xl)*rhx
632:             d2 = (xc-xr)*rhx
633:             d3 = (xc-xt)*rhy
634:             d4 = (xc-xb)*rhy
635:             d5 = (xrb-xr)*rhy
636:             d6 = (xrb-xb)*rhx
637:             d7 = (xlt-xl)*rhy
638:             d8 = (xlt-xt)*rhx
639:
640:             f1 = sqrt( 1.0d0 + d1*d1 + d7*d7)
641:             f2 = sqrt( 1.0d0 + d1*d1 + d4*d4)
642:             f3 = sqrt( 1.0d0 + d3*d3 + d8*d8)
643:             f4 = sqrt( 1.0d0 + d3*d3 + d2*d2)
644:             f5 = sqrt( 1.0d0 + d2*d2 + d5*d5)
645:             f6 = sqrt( 1.0d0 + d4*d4 + d6*d6)
646:
647:
648:             hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+                 &
649:      &              (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)

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

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

657:             hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+                 &
658:      &              (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
659:
660:             hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
661:             htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
662:
663:             hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) +                         &
664:      &              hxdhy*(1.0+d8*d8)/(f3*f3*f3) +                      &
665:      &              hydhx*(1.0+d5*d5)/(f5*f5*f5) +                      &
666:      &              hxdhy*(1.0+d6*d6)/(f6*f6*f6) +                      &
667:      &              (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-               &
668:      &              2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+          &
669:      &              hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
670:
671:             hl = hl * 0.5
672:             hr = hr * 0.5
673:             ht = ht * 0.5
674:             hb = hb * 0.5
675:             hbr = hbr * 0.5
676:             htl = htl * 0.5
677:             hc = hc * 0.5

679:             k = 0

681:             if (j .gt. 0) then
682:                v(k) = hb
683:                col(k) = row - gxm
684:                k=k+1
685:             endif

687:             if ((j .gt. 0) .and. (i .lt. mx-1)) then
688:                v(k) = hbr
689:                col(k) = row-gxm+1
690:                k=k+1
691:             endif

693:             if (i .gt. 0) then
694:                v(k) = hl
695:                col(k) = row - 1
696:                k = k+1
697:             endif

699:             v(k) = hc
700:             col(k) = row
701:             k=k+1

703:             if (i .lt. mx-1) then
704:                v(k) = hr
705:                col(k) = row + 1
706:                k=k+1
707:             endif

709:             if ((i .gt. 0) .and. (j .lt. my-1)) then
710:                v(k) = htl
711:                col(k) = row + gxm - 1
712:                k=k+1
713:             endif

715:             if (j .lt. my-1) then
716:                v(k) = ht
717:                col(k) = row + gxm
718:                k=k+1
719:             endif

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

725:

727:          enddo
728:       enddo
729:
730: ! restore vectors
731:       call VecRestoreArray(localX,x_v,x_i,ierr)
732:       call VecRestoreArray(Left,left_v,left_i,ierr)
733:       call VecRestoreArray(Right,right_v,right_i,ierr)
734:       call VecRestoreArray(Top,top_v,top_i,ierr)
735:       call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)

738: ! Assemble the matrix
739:       call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
740:       call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)

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

744:       return
745:       end
746:
747:

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

753: ! ----------------------------------------------------------------------------
754: !
755: !/*
756: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
757: !
758: !
759: !*/

761:       subroutine MSA_BoundaryConditions(ierr)
762:       implicit none

764: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
765: #include "plate2f.h"

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

783:       limit=0
784:       maxits = 5
785:       tol=1e-10
786:       b=-0.5d0
787:       t= 0.5d0
788:       l=-0.5d0
789:       r= 0.5d0
790:       xt=0
791:       yt=0
792:       one=1.0d0
793:       two=2.0d0
794:       three=3.0d0

797:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
798:      &                  PETSC_NULL_INTEGER,ierr)
799:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
800:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

802:       bsize = xm + 2
803:       lsize = ym + 2
804:       rsize = ym + 2
805:       tsize = xm + 2
806:

808:       call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
809:       call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
810:       call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
811:       call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)

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

816:       do j=0,3
817:
818:          if (j.eq.0) then
819:             yt=b
820:             xt=l+hx*xs
821:             limit=bsize
822:             call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
823:

825:          elseif (j.eq.1) then
826:             yt=t
827:             xt=l+hx*xs
828:             limit=tsize
829:             call VecGetArray(Top,boundary_v,boundary_i,ierr)

831:          elseif (j.eq.2) then
832:             yt=b+hy*ys
833:             xt=l
834:             limit=lsize
835:             call VecGetArray(Left,boundary_v,boundary_i,ierr)

837:          elseif (j.eq.3) then
838:             yt=b+hy*ys
839:             xt=r
840:             limit=rsize
841:             call VecGetArray(Right,boundary_v,boundary_i,ierr)
842:          endif
843:

845:          do i=0,limit-1
846:
847:             u1=xt
848:             u2=-yt
849:             k = 0
850:             exitloop = .false.
851:             do while (k .lt. maxits .and. (.not. exitloop) )

853:                nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
854:                nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
855:                fnorm=sqrt(nf1*nf1+nf2*nf2)
856:                if (fnorm .gt. tol) then
857:                   njac11=one+u2*u2-u1*u1
858:                   njac12=two*u1*u2
859:                   njac21=-two*u1*u2
860:                   njac22=-one - u1*u1 + u2*u2
861:                   det = njac11*njac22-njac21*njac12
862:                   u1 = u1-(njac22*nf1-njac12*nf2)/det
863:                   u2 = u2-(njac11*nf2-njac21*nf1)/det
864:                else
865:                   exitloop = .true.
866:                endif
867:                k=k+1
868:             enddo

870:             boundary_v(i + boundary_i) = u1*u1-u2*u2
871:             if ((j .eq. 0) .or. (j .eq. 1)) then
872:                xt = xt + hx
873:             else
874:                yt = yt + hy
875:             endif

877:          enddo
878:

880:          if (j.eq.0) then
881:             call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
882:          elseif (j.eq.1) then
883:             call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
884:          elseif (j.eq.2) then
885:             call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
886:          elseif (j.eq.3) then
887:             call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
888:          endif
889:
890:       enddo

893: ! Scale the boundary if desired
894:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom",            &
895:      &                         scl,flg,ierr)
896:       if (flg .eqv. PETSC_TRUE) then
897:          call VecScale(scl,Bottom,ierr)
898:       endif

900:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top",               &
901:      &                         scl,flg,ierr)
902:       if (flg .eqv. PETSC_TRUE) then
903:          call VecScale(scl,Top,ierr)
904:       endif

906:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right",             &
907:      &                         scl,flg,ierr)
908:       if (flg .eqv. PETSC_TRUE) then
909:          call VecScale(scl,Right,ierr)
910:       endif

912:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left",              &
913:      &                         scl,flg,ierr)
914:       if (flg .eqv. PETSC_TRUE) then
915:          call VecScale(scl,Left,ierr)
916:       endif
917:
918:
919:       return
920:       end

922: ! ----------------------------------------------------------------------------
923: !
924: !/*
925: !     MSA_Plate - Calculates an obstacle for surface to stretch over
926: !
927: !     Output Parameter:
928: !.    xl - lower bound vector
929: !.    xu - upper bound vector
930: !
931: !*/

933:       subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
934:       implicit none
935: ! mx,my,bmx,bmy,dm,bheight defined in plate2f.h
936: #include "plate2f.h"
937:       TaoSolver        tao
938:       Vec              xl,xu
939:       PetscErrorCode   ierr
940:       PetscInt         i,j,row
941:       PetscInt         xs, xm, ys, ym
942:       PetscReal      lb,ub
943:       PetscInt         dummy

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

953:       print *,'msa_plate'

955:       lb = -1.0d300
956:       ub = 1.0d300

958:       if (bmy .lt. 0) bmy = 0
959:       if (bmy .gt. my) bmy = my
960:       if (bmx .lt. 0) bmx = 0
961:       if (bmx .gt. mx) bmx = mx
962:

964:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
965:      &             PETSC_NULL_INTEGER,ierr)

967:       call VecSet(xl,lb,ierr)
968:       call VecSet(xu,ub,ierr)

970:       call VecGetArray(xl,xl_v,xl_i,ierr)
971:

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

975:          do j=ys,ys+ym-1
976:
977:             row=(j-ys)*xm + (i-xs)

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

983:             endif

985:          enddo
986:       enddo

989:       call VecRestoreArray(xl,xl_v,xl_i,ierr)
990:
991:       return
992:       end

996:
997:
998: ! ----------------------------------------------------------------------------
999: !
1000: !/*
1001: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
1002: !
1003: !     Output Parameter:
1004: !.    X - vector for initial guess
1005: !
1006: !*/

1008:       subroutine MSA_InitialPoint(X, ierr)
1009:       implicit none

1011: ! mx,my,localX,dm,Top,Left,Bottom,Right defined in plate2f.h
1012: #include "plate2f.h"
1013:       Vec               X
1014:       PetscErrorCode    ierr
1015:       PetscInt          start,i,j
1016:       PetscInt          row
1017:       PetscInt          xs,xm,gxs,gxm
1018:       PetscInt          ys,ym,gys,gym
1019:       PetscReal       zero, np5

1021: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1022: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
1023: ! will return an array of doubles referenced by x_array offset by x_index.
1024: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
1025: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1026:       PetscReal   left_v(0:1),right_v(0:1)
1027:       PetscReal   bottom_v(0:1),top_v(0:1)
1028:       PetscReal   x_v(0:1)
1029:       PetscOffset   left_i, right_i, top_i
1030:       PetscOffset   bottom_i,x_i
1031:       PetscBool     flg
1032:       PetscRandom   rctx
1033:
1034:       zero = 0.0d0
1035:       np5 = -0.5d0

1037:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start",            &
1038:      &                        start,flg,ierr)

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

1043:       elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then  ! random start -0.5 < xi < 0.5
1044:          call PetscRandomCreate(MPI_COMM_WORLD,rctx,ierr)
1045:          do i=0,start-1
1046:             call VecSetRandom(X,rctx,ierr)
1047:          enddo

1049:          call PetscRandomDestroy(rctx,ierr)
1050:          call VecShift(X,np5,ierr)

1052:       else   ! average of boundary conditions
1053:
1054: !        Get Local mesh boundaries
1055:          call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,             &
1056:      &                     PETSC_NULL_INTEGER,ierr)
1057:          call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,    &
1058:      &                     PETSC_NULL_INTEGER,ierr)

1062: !        Get pointers to vector data
1063:          call VecGetArray(Top,top_v,top_i,ierr)
1064:          call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1065:          call VecGetArray(Left,left_v,left_i,ierr)
1066:          call VecGetArray(Right,right_v,right_i,ierr)
1067:
1068:          call VecGetArray(localX,x_v,x_i,ierr)
1069:
1070: !        Perform local computations
1071:          do  j=ys,ys+ym-1
1072:             do i=xs,xs+xm-1
1073:                row = (j-gys)*gxm  + (i-gxs)
1074:                x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my        &
1075:      &             + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) +                  &
1076:      &              (i+1)*left_v(left_i+j-ys+1)/mx       +                  &
1077:      &              (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1078:             enddo
1079:          enddo

1081: !        Restore vectors
1082:          call VecRestoreArray(localX,x_v,x_i,ierr)

1084:          call VecRestoreArray(Left,left_v,left_i,ierr)
1085:          call VecRestoreArray(Top,top_v,top_i,ierr)
1086:          call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1087:          call VecRestoreArray(Right,right_v,right_i,ierr)

1089:          call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1090:          call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)

1092:       endif

1094:       return
1095:       end
```