Actual source code: plate2f.F

  1: ! "$Id$";
  2: !
  3: !  Program usage: mpirun -np <proc> plate2f [all TAO options] 
  4: !
  5: !  This example demonstrates use of the TAO package to solve a bound constrained
  6: !  minimization problem.  This example is based on a problem from the
  7: !  MINPACK-2 test suite.  Given a rectangular 2-D domain and boundary values
  8: !  along the edges of the domain, the objective is to find the surface
  9: !  with the minimal area that satisfies the boundary conditions.
 10: !  The command line options are:
 11: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
 12: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
 13: !    -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
 14: !    -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
 15: !    -bheight <ht>, where <ht> = height of the plate
 16: !
 17: !/*T
 18: !   Concepts: TAO - Solving a bound constrained minimization problem
 19: !   Routines: TaoInitialize(); TaoFinalize();
 20: !   Routines: TaoCreate(); TaoDestroy();
 21: !   Routines: TaoAppSetObjectiveAndGradientRoutine(); 
 22: !   Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
 23: !   Routines: TaoAppSetInitialSolutionVec(); TaoAppSetVariableBounds();
 24: !   Routines: TaoSetOptions();
 25: !   Routines: TaoApplicationCreate(); TaoSolve();
 26: !   Routines: TaoView(); TaoAppDestroy();
 27: !   Processors: n
 28: !T*/



 32:       implicit none

 34: #include "plate2f.h"

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

 46:       integer          info          ! used to check for functions returning nonzeros
 47:       Vec              x             ! solution vector
 48:       Vec              xl, xu        ! lower and upper bounds vectorsp
 49:       PetscInt         m             ! number of local elements in vector
 50:       TAO_SOLVER       tao           ! TAO_SOLVER solver context
 51:       TAO_APPLICATION  plateapp      ! PETSc application
 52:       Mat              H             ! Hessian matrix
 53:       ISLocalToGlobalMapping isltog  ! local to global mapping object
 54:       PetscTruth       flg
 55:       PetscInt         i1,i3,i7


 58:       external FormFunctionGradient
 59:       external FormHessian
 60:       external MSA_BoundaryConditions
 61:       external MSA_Plate
 62:       external MSA_InitialPoint
 63: ! Initialize Tao

 65:       i1=1
 66:       i3=3
 67:       i7=7


 70:       call PetscInitialize(PETSC_NULL_CHARACTER,info)
 71:       call TaoInitialize(PETSC_NULL_CHARACTER,info)
 72:       

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

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

 87:       call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmx",bmx,flg,info)
 88:       call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmy",bmy,flg,info)
 89:       call PetscOptionsGetReal(TAO_NULL_CHARACTER,"-bheight",bheight,   &
 90:      &      flg,info)
 91:       

 93: ! Calculate any derived values from parameters
 94:       N = mx*my

 96: ! Let Petsc determine the dimensions of the local vectors 
 97:       Nx = PETSC_DECIDE
 98:       NY = PETSC_DECIDE

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

104:       call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,         &
105:      &     mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,           &
106:      &     da,info)
107:       

109: ! Extract global and local vectors from DA; The local vectors are
110: ! used solely as work space for the evaluation of the function, 
111: ! gradient, and Hessian.  Duplicate for remaining vectors that are 
112: ! the same types.

114:       call DACreateGlobalVector(da,x,info)
115:       call DACreateLocalVector(da,localX,info)
116:       call VecDuplicate(localX,localV,info)

118: ! Create a matrix data structure to store the Hessian.
119: ! Here we (optionally) also associate the local numbering scheme
120: ! with the matrix so that later we can use local indices for matrix
121: ! assembly

123:       call VecGetLocalSize(x,m,info)
124:       call MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER,  &
125:      &     i3,PETSC_NULL_INTEGER,H,info)

127:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)
128:       call DAGetISLocalToGlobalMapping(da,isltog,info)
129:       call MatSetLocalToGlobalMapping(H,isltog,info)
130:       

132: ! The Tao code begins here
133: ! Create TAO solver and set desired solution method.
134: ! This problems uses bounded variables, so the
135: ! method must either be 'tao_tron' or 'tao_blmvm'

137:       call TaoCreate(MPI_COMM_WORLD,'tao_blmvm',tao,info)
138:       call TaoApplicationCreate(MPI_COMM_WORLD,plateapp,info)
139:       

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

143: !     TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
144:       call TaoAppSetObjectiveAndGradientRo(plateapp,                    &
145:      &     FormFunctionGradient,PETSC_NULL_OBJECT,info)

147:       call TaoAppSetHessianMat(plateapp,H,H,info)
148:       call TaoAppSetHessianRoutine(plateapp,FormHessian,                 &
149:      &     PETSC_NULL_OBJECT, info)
150:       

152: ! Set Variable bounds 
153:       call MSA_BoundaryConditions(info)
154:       call VecDuplicate(x,xl,info)
155:       call VecDuplicate(x,xu,info)
156:       call MSA_Plate(xl,xu,info)
157:       call TaoAppSetVariableBounds(plateapp,xl,xu,info)

159: ! Set the initial solution guess
160:       call MSA_InitialPoint(x, info)
161:       call TaoAppSetInitialSolutionVec(plateapp,x,info)

163: ! Check for any tao command line options
164:       call TaoSetOptions(plateapp,tao,info)

166: ! Solve the application
167:       call TaoSolveApplication(plateapp,tao,info)

169: !     View TAO solver information
170: !      call TaoView(tao,info)

172: ! Free TAO data structures
173:       call TaoDestroy(tao,info)
174:       call TaoAppDestroy(plateapp,info)  

176: ! Free PETSc data structures
177:       call VecDestroy(x,info)
178:       call VecDestroy(xl,info)
179:       call VecDestroy(xu,info)
180:       call VecDestroy(Top,info)
181:       call VecDestroy(Bottom,info)
182:       call VecDestroy(Left,info)
183:       call VecDestroy(Right,info)
184:       call MatDestroy(H,info)
185:       call VecDestroy(localX,info)
186:       call VecDestroy(localV,info)
187:       call DADestroy(da,info)

189: ! Finalize TAO 

191:       call TaoFinalize(info)
192:       call PetscFinalize(info)

194:       end

196: ! ---------------------------------------------------------------------
197: !
198: !  FormFunctionGradient - Evaluates function f(X). 
199: !    
200: !  Input Parameters:
201: !  tao   - the TAO_SOLVER context
202: !  X     - the input vector 
203: !  dummy - optional user-defined context, as set by TaoSetFunction()
204: !          (not used here)
205: !
206: !  Output Parameters:
207: !  fcn     - the newly evaluated function
208: !  G       - the gradient vector
209: !  info  - error code
210: !


213:       subroutine FormFunctionGradient(tao,X,fcn,G,dummy,info)
214:       implicit none

216: ! da, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
217: #include "plate2f.h"
218:       
219: ! Input/output variables

221:       TAO_SOLVER       tao
222:       PetscScalar fcn
223:       Vec              X, G
224:       integer          info
225:       PetscInt         dummy
226:       
227:       PetscInt         i,j,row
228:       PetscInt         xs, xm
229:       PetscInt         gxs, gxm
230:       PetscInt         ys, ym
231:       PetscInt         gys, gym
232:       PetscScalar      ft,zero,hx,hy,hydhx,hxdhy
233:       PetscScalar      area,rhx,rhy
234:       PetscScalar      f1,f2,f3,f4,f5,f6,d1,d2,d3
235:       PetscScalar      d4,d5,d6,d7,d8
236:       PetscScalar      df1dxc,df2dxc,df3dxc,df4dxc
237:       PetscScalar      df5dxc,df6dxc
238:       PetscScalar      xc,xl,xr,xt,xb,xlt,xrb


241: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
242: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
243: ! will return an array of doubles referenced by x_array offset by x_index.
244: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
245: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
246:       PetscScalar      g_v(0:1),x_v(0:1)
247:       PetscScalar      top_v(0:1),left_v(0:1)
248:       PetscScalar      right_v(0:1),bottom_v(0:1)
249:       PetscOffset      g_i,left_i,right_i
250:       PetscOffset      bottom_i,top_i,x_i

252:       ft = 0.0d0
253:       zero = 0.0d0
254:       hx = 1.0d0/(mx + 1)
255:       hy = 1.0d0/(my + 1)
256:       hydhx = hy/hx
257:       hxdhy = hx/hy
258:       area = 0.5d0 * hx * hy
259:       rhx = mx + 1.0d0
260:       rhy = my + 1.0d0


263: ! Get local mesh boundaries
264:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
265:      &                  PETSC_NULL_INTEGER,info)
266:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,             &
267:      &                       gxm,gym,PETSC_NULL_INTEGER,info)

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

273: ! Initialize the vector to zero
274:       call VecSet(localV,zero,info)

276: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
277:       call VecGetArray(localX,x_v,x_i,info)
278:       call VecGetArray(localV,g_v,g_i,info)
279:       call VecGetArray(Top,top_v,top_i,info)
280:       call VecGetArray(Bottom,bottom_v,bottom_i,info)
281:       call VecGetArray(Left,left_v,left_i,info)
282:       call VecGetArray(Right,right_v,right_i,info)

284: ! Compute function over the locally owned part of the mesh       
285:       do j = ys,ys+ym-1
286:          do i = xs,xs+xm-1
287:             row = (j-gys)*gxm + (i-gxs)
288:             xc = x_v(row+x_i)
289:             xt = xc
290:             xb = xc
291:             xr = xc
292:             xl = xc
293:             xrb = xc
294:             xlt = xc

296:             if (i .eq. 0) then !left side
297:                xl = left_v(j - ys + 1 + left_i)
298:                xlt = left_v(j - ys + 2 + left_i)
299:             else
300:                xl = x_v(row - 1 + x_i)
301:             endif

303:             if (j .eq. 0) then !bottom side
304:                xb = bottom_v(i - xs + 1 + bottom_i)
305:                xrb = bottom_v(i - xs + 2 + bottom_i)
306:             else
307:                xb = x_v(row - gxm + x_i)
308:             endif

310:             if (i + 1 .eq. gxs + gxm) then !right side
311:                xr = right_v(j - ys + 1 + right_i)
312:                xrb = right_v(j - ys + right_i)
313:             else
314:                xr = x_v(row + 1 + x_i)
315:             endif

317:             if (j + 1 .eq. gys + gym) then !top side
318:                xt = top_v(i - xs + 1 + top_i)
319:                xlt = top_v(i - xs + top_i)
320:             else
321:                xt = x_v(row + gxm + x_i)
322:             endif

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

328:             if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
329:                xrb = x_v(row + 1 - gxm + x_i)
330:             endif

332:             d1 = xc-xl
333:             d2 = xc-xr
334:             d3 = xc-xt
335:             d4 = xc-xb
336:             d5 = xr-xrb
337:             d6 = xrb-xb
338:             d7 = xlt-xl
339:             d8 = xt-xlt

341:             df1dxc = d1 * hydhx
342:             df2dxc = d1 * hydhx + d4 * hxdhy
343:             df3dxc = d3 * hxdhy
344:             df4dxc = d2 * hydhx + d3 * hxdhy
345:             df5dxc = d2 * hydhx
346:             df6dxc = d4 * hxdhy

348:             d1 = d1 * rhx
349:             d2 = d2 * rhx
350:             d3 = d3 * rhy
351:             d4 = d4 * rhy
352:             d5 = d5 * rhy
353:             d6 = d6 * rhx
354:             d7 = d7 * rhy
355:             d8 = d8 * rhx

357:             f1 = sqrt(1.0d0 + d1*d1 + d7*d7)
358:             f2 = sqrt(1.0d0 + d1*d1 + d4*d4)
359:             f3 = sqrt(1.0d0 + d3*d3 + d8*d8)
360:             f4 = sqrt(1.0d0 + d3*d3 + d2*d2)
361:             f5 = sqrt(1.0d0 + d2*d2 + d5*d5)
362:             f6 = sqrt(1.0d0 + d4*d4 + d6*d6)

364:             ft = ft + f2 + f4

366:             df1dxc = df1dxc / f1
367:             df2dxc = df2dxc / f2
368:             df3dxc = df3dxc / f3
369:             df4dxc = df4dxc / f4
370:             df5dxc = df5dxc / f5
371:             df6dxc = df6dxc / f6

373:             g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc +  &
374:      &                              df5dxc + df6dxc)
375:          enddo
376:       enddo

378: ! Compute triangular areas along the border of the domain. 
379:       if (xs .eq. 0) then  ! left side
380:          do j=ys,ys+ym-1
381:             d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i))         &
382:      &                 * rhy
383:             d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i))        &
384:      &                 * rhx
385:             ft = ft + sqrt(1.0d0 + d3*d3 + d2*d2)
386:          enddo
387:       endif

389:       
390:       if (ys .eq. 0) then !bottom side
391:          do i=xs,xs+xm-1
392:             d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i))    &
393:      &                    * rhx
394:             d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
395:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
396:          enddo
397:       endif

399:       
400:       if (xs + xm .eq. mx) then ! right side
401:          do j=ys,ys+ym-1
402:             d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
403:             d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
404:             ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
405:          enddo
406:       endif

408:       
409:       if (ys + ym .eq. my) then
410:          do i=xs,xs+xm-1
411:             d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
412:             d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
413:             ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
414:          enddo
415:       endif

417:       
418:       if ((ys .eq. 0) .and. (xs .eq. 0)) then
419:          d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
420:          d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
421:          ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
422:       endif

424:       if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
425:          d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
426:          d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
427:          ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
428:       endif

430:       ft = ft * area
431:       call MPI_Allreduce(ft,fcn,1,MPI_DOUBLE_PRECISION,                  &
432:      &             MPI_SUM,MPI_COMM_WORLD,info)



436: ! Restore vectors
437:       call VecRestoreArray(localX,x_v,x_i,info)
438:       call VecRestoreArray(localV,g_v,g_i,info)
439:       call VecRestoreArray(Left,left_v,left_i,info)
440:       call VecRestoreArray(Top,top_v,top_i,info)
441:       call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
442:       call VecRestoreArray(Right,right_v,right_i,info)

444: ! Scatter values to global vector
445:       call DALocalToGlobal(da,localV,INSERT_VALUES,G,info)

447:       call PetscLogFlops(70*xm*ym,info)

449:       return
450:       end  !FormFunctionGradient
451:       




456: ! ----------------------------------------------------------------------------
457: ! 
458: !/*
459: !   FormHessian - Evaluates Hessian matrix.
460: !
461: !   Input Parameters:
462: !.  tao  - the TAO_SOLVER context
463: !.  X    - input vector
464: !.  dummy  - not used 
465: !
466: !   Output Parameters:
467: !.  Hessian    - Hessian matrix
468: !.  Hpc    - optionally different preconditioning matrix
469: !.  flag - flag indicating matrix structure
470: !
471: !   Notes:
472: !   Due to mesh point reordering with DAs, we must always work
473: !   with the local mesh points, and then transform them to the new
474: !   global numbering with the local-to-global mapping.  We cannot work
475: !   directly with the global numbers for the original uniprocessor mesh!  
476: !
477: !   Two methods are available for imposing this transformation
478: !   when setting matrix entries:
479: !     (A) MatSetValuesLocal(), using the local ordering (including
480: !         ghost points!)
481: !         - Do the following two steps once, before calling TaoSolve()
482: !           - Use DAGetISLocalToGlobalMapping() to extract the
483: !             local-to-global map from the DA
484: !           - Associate this map with the matrix by calling
485: !             MatSetLocalToGlobalMapping() 
486: !         - Then set matrix entries using the local ordering
487: !           by calling MatSetValuesLocal()
488: !     (B) MatSetValues(), using the global ordering 
489: !         - Use DAGetGlobalIndices() to extract the local-to-global map
490: !         - Then apply this map explicitly yourself
491: !         - Set matrix entries using the global ordering by calling
492: !           MatSetValues()
493: !   Option (A) seems cleaner/easier in many cases, and is the procedure
494: !   used in this example.
495: */
496:       subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, info)
497:       implicit none

499: ! da,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
500: #include "plate2f.h"
501:       
502:       TAO_SOLVER     tao
503:       Vec            X
504:       Mat            Hessian,Hpc
505:       MatStructure   flg
506:       PetscInt       dummy
507:       integer        info

509:       PetscInt       i,j,k,row
510:       PetscInt       xs,xm,gxs,gxm
511:       PetscInt       ys,ym,gys,gym
512:       PetscInt       col(0:6)
513:       PetscScalar    hx,hy,hydhx,hxdhy,rhx,rhy
514:       PetscScalar    f1,f2,f3,f4,f5,f6,d1,d2,d3
515:       PetscScalar    d4,d5,d6,d7,d8
516:       PetscScalar    xc,xl,xr,xt,xb,xlt,xrb
517:       PetscScalar    hl,hr,ht,hb,hc,htl,hbr

519: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
520: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
521: ! will return an array of doubles referenced by x_array offset by x_index.
522: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
523: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
524:       PetscScalar   right_v(0:1),left_v(0:1)
525:       PetscScalar   bottom_v(0:1),top_v(0:1)
526:       PetscScalar   x_v(0:1)
527:       PetscOffset   x_i,right_i,left_i
528:       PetscOffset   bottom_i,top_i
529:       PetscScalar   v(0:6)
530:       PetscTruth    assembled
531:       PetscInt      i1
532:  
533:       i1=1
534:       
535: ! Set various matrix options 
536:       call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,
537:      &                  PETSC_TRUE,info)

539: ! Get local mesh boundaries
540:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
541:      &                  PETSC_NULL_INTEGER,info)
542:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
543:      &                       PETSC_NULL_INTEGER,info)

545: ! Scatter ghost points to local vectors
546:       call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
547:       call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)

549: ! Get pointers to vector data (see note on Fortran arrays above)
550:       call VecGetArray(localX,x_v,x_i,info)
551:       call VecGetArray(Top,top_v,top_i,info)
552:       call VecGetArray(Bottom,bottom_v,bottom_i,info)
553:       call VecGetArray(Left,left_v,left_i,info)
554:       call VecGetArray(Right,right_v,right_i,info)

556: ! Initialize matrix entries to zero
557:       call MatAssembled(Hessian,assembled,info)
558:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,info)


561:       rhx = mx + 1.0
562:       rhy = my + 1.0
563:       hx = 1.0/rhx
564:       hy = 1.0/rhy
565:       hydhx = hy/hx
566:       hxdhy = hx/hy
567: ! compute Hessian over the locally owned part of the mesh

569:       do  i=xs,xs+xm-1
570:          do  j=ys,ys+ym-1
571:             row = (j-gys)*gxm + (i-gxs)
572:             
573:             xc = x_v(row + x_i)
574:             xt = xc
575:             xb = xc
576:             xr = xc
577:             xl = xc
578:             xrb = xc
579:             xlt = xc

581:             if (i .eq. gxs) then   ! Left side
582:                xl = left_v(left_i + j - ys + 1)
583:                xlt = left_v(left_i + j - ys + 2)
584:             else
585:                xl = x_v(x_i + row -1 )
586:             endif

588:             if (j .eq. gys) then ! bottom side
589:                xb = bottom_v(bottom_i + i - xs + 1)
590:                xrb = bottom_v(bottom_i + i - xs + 2)
591:             else
592:                xb = x_v(x_i + row - gxm)
593:             endif

595:             if (i+1 .eq. gxs + gxm) then !right side
596:                xr = right_v(right_i + j - ys + 1)
597:                xrb = right_v(right_i + j - ys)
598:             else
599:                xr = x_v(x_i + row + 1)
600:             endif

602:             if (j+1 .eq. gym+gys) then !top side
603:                xt = top_v(top_i +i - xs + 1)
604:                xlt = top_v(top_i + i - xs)
605:             else
606:                xt = x_v(x_i + row + gxm)
607:             endif

609:             if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
610:                xlt = x_v(x_i + row - 1 + gxm)
611:             endif
612:             
613:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
614:                xrb = x_v(x_i + row + 1 - gxm)
615:             endif

617:             d1 = (xc-xl)*rhx
618:             d2 = (xc-xr)*rhx
619:             d3 = (xc-xt)*rhy
620:             d4 = (xc-xb)*rhy
621:             d5 = (xrb-xr)*rhy
622:             d6 = (xrb-xb)*rhx
623:             d7 = (xlt-xl)*rhy
624:             d8 = (xlt-xt)*rhx
625:             
626:             f1 = sqrt( 1.0d0 + d1*d1 + d7*d7)
627:             f2 = sqrt( 1.0d0 + d1*d1 + d4*d4)
628:             f3 = sqrt( 1.0d0 + d3*d3 + d8*d8)
629:             f4 = sqrt( 1.0d0 + d3*d3 + d2*d2)
630:             f5 = sqrt( 1.0d0 + d2*d2 + d5*d5)
631:             f6 = sqrt( 1.0d0 + d4*d4 + d6*d6)
632:             
633:             
634:             hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+                 &
635:      &              (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)

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

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

643:             hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+                 &
644:      &              (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
645:             
646:             hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
647:             htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
648:             
649:             hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) +                         &
650:      &              hxdhy*(1.0+d8*d8)/(f3*f3*f3) +                      &
651:      &              hydhx*(1.0+d5*d5)/(f5*f5*f5) +                      &
652:      &              hxdhy*(1.0+d6*d6)/(f6*f6*f6) +                      &
653:      &              (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-               &
654:      &              2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+          &
655:      &              hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)               
656:             
657:             hl = hl * 0.5
658:             hr = hr * 0.5
659:             ht = ht * 0.5
660:             hb = hb * 0.5
661:             hbr = hbr * 0.5
662:             htl = htl * 0.5
663:             hc = hc * 0.5 

665:             k = 0

667:             if (j .gt. 0) then
668:                v(k) = hb
669:                col(k) = row - gxm
670:                k=k+1
671:             endif

673:             if ((j .gt. 0) .and. (i .lt. mx-1)) then
674:                v(k) = hbr
675:                col(k) = row-gxm+1
676:                k=k+1
677:             endif

679:             if (i .gt. 0) then
680:                v(k) = hl
681:                col(k) = row - 1
682:                k = k+1
683:             endif

685:             v(k) = hc
686:             col(k) = row
687:             k=k+1

689:             if (i .lt. mx-1) then
690:                v(k) = hr
691:                col(k) = row + 1
692:                k=k+1
693:             endif

695:             if ((i .gt. 0) .and. (j .lt. my-1)) then
696:                v(k) = htl
697:                col(k) = row + gxm - 1
698:                k=k+1
699:             endif

701:             if (j .lt. my-1) then
702:                v(k) = ht
703:                col(k) = row + gxm
704:                k=k+1
705:             endif

707: ! Set matrix values using local numbering, defined earlier in main routine
708:             call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES,      &
709:      &                              info)

711:             

713:          enddo
714:       enddo
715:       
716: ! restore vectors
717:       call VecRestoreArray(localX,x_v,x_i,info)
718:       call VecRestoreArray(Left,left_v,left_i,info)
719:       call VecRestoreArray(Right,right_v,right_i,info)
720:       call VecRestoreArray(Top,top_v,top_i,info)
721:       call VecRestoreArray(Bottom,bottom_v,bottom_i,info)


724: ! Assemble the matrix
725:       call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,info)
726:       call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,info)

728:       call PetscLogFlops(199*xm*ym,info)

730:       return
731:       end  
732:       
733:       



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

739: ! ----------------------------------------------------------------------------
740: !
741: !/*
742: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
743: !
744: !
745: !*/

747:       subroutine MSA_BoundaryConditions(info)
748:       implicit none

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

753:       integer          info
754:       PetscInt         i,j,k,limit,maxits
755:       PetscInt         xs, xm, gxs, gxm
756:       PetscInt         ys, ym, gys, gym
757:       PetscInt         bsize, lsize
758:       PetscInt         tsize, rsize
759:       PetscScalar      one,two,three,tol
760:       PetscScalar      scl,fnorm,det,xt
761:       PetscScalar      yt,hx,hy,u1,u2,nf1,nf2
762:       PetscScalar      njac11,njac12,njac21,njac22
763:       PetscScalar      b, t, l, r
764:       PetscScalar      boundary_v(0:1)
765:       PetscOffset      boundary_i
766:       logical exitloop
767:       TaoTruth flg

769:       limit=0
770:       maxits = 5
771:       tol=1e-10
772:       b=-0.5d0
773:       t= 0.5d0
774:       l=-0.5d0
775:       r= 0.5d0
776:       xt=0
777:       yt=0
778:       one=1.0d0
779:       two=2.0d0
780:       three=3.0d0


783:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
784:      &                  PETSC_NULL_INTEGER,info)
785:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,              &
786:      &                       gxm,gym,PETSC_NULL_INTEGER,info)

788:       bsize = xm + 2
789:       lsize = ym + 2
790:       rsize = ym + 2
791:       tsize = xm + 2
792:       

794:       call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,info)
795:       call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,info)
796:       call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,info)
797:       call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,info)

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

802:       do j=0,3
803:          
804:          if (j.eq.0) then
805:             yt=b
806:             xt=l+hx*xs
807:             limit=bsize
808:             call VecGetArray(Bottom,boundary_v,boundary_i,info)
809:             

811:          elseif (j.eq.1) then
812:             yt=t
813:             xt=l+hx*xs
814:             limit=tsize
815:             call VecGetArray(Top,boundary_v,boundary_i,info)

817:          elseif (j.eq.2) then
818:             yt=b+hy*ys
819:             xt=l
820:             limit=lsize
821:             call VecGetArray(Left,boundary_v,boundary_i,info)

823:          elseif (j.eq.3) then
824:             yt=b+hy*ys
825:             xt=r
826:             limit=rsize
827:             call VecGetArray(Right,boundary_v,boundary_i,info)
828:          endif
829:          

831:          do i=0,limit-1
832:             
833:             u1=xt
834:             u2=-yt
835:             k = 0
836:             exitloop = .false.
837:             do while (k .lt. maxits .and. (.not. exitloop) )

839:                nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
840:                nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
841:                fnorm=sqrt(nf1*nf1+nf2*nf2)
842:                if (fnorm .gt. tol) then
843:                   njac11=one+u2*u2-u1*u1
844:                   njac12=two*u1*u2
845:                   njac21=-two*u1*u2
846:                   njac22=-one - u1*u1 + u2*u2
847:                   det = njac11*njac22-njac21*njac12
848:                   u1 = u1-(njac22*nf1-njac12*nf2)/det
849:                   u2 = u2-(njac11*nf2-njac21*nf1)/det
850:                else 
851:                   exitloop = .true.
852:                endif
853:                k=k+1
854:             enddo

856:             boundary_v(i + boundary_i) = u1*u1-u2*u2
857:             if ((j .eq. 0) .or. (j .eq. 1)) then
858:                xt = xt + hx
859:             else
860:                yt = yt + hy
861:             endif

863:          enddo
864:                

866:          if (j.eq.0) then
867:             call VecRestoreArray(Bottom,boundary_v,boundary_i,info)
868:          elseif (j.eq.1) then
869:             call VecRestoreArray(Top,boundary_v,boundary_i,info)
870:          elseif (j.eq.2) then
871:             call VecRestoreArray(Left,boundary_v,boundary_i,info)
872:          elseif (j.eq.3) then
873:             call VecRestoreArray(Right,boundary_v,boundary_i,info)
874:          endif
875:          
876:       enddo


879: ! Scale the boundary if desired
880:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom",            &
881:      &                         scl,flg,info)
882:       if (flg .eqv. PETSC_TRUE) then
883:          call VecScale(scl,Bottom,info)
884:       endif

886:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top",               &
887:      &                         scl,flg,info)
888:       if (flg .eqv. PETSC_TRUE) then
889:          call VecScale(scl,Top,info)
890:       endif

892:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right",             &
893:      &                         scl,flg,info)
894:       if (flg .eqv. PETSC_TRUE) then
895:          call VecScale(scl,Right,info)
896:       endif

898:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left",              &
899:      &                         scl,flg,info)
900:       if (flg .eqv. PETSC_TRUE) then
901:          call VecScale(scl,Left,info)
902:       endif
903:          
904:       
905:       return
906:       end

908: ! ----------------------------------------------------------------------------
909: !
910: !/*
911: !     MSA_Plate - Calculates an obstacle for surface to stretch over
912: !
913: !     Output Parameter:
914: !.    xl - lower bound vector
915: !.    xu - upper bound vector
916: !
917: !*/

919:       subroutine MSA_Plate(xl,xu,info)
920:       implicit none

922: ! mx,my,bmx,bmy,da,bheight defined in plate2f.h
923: #include "plate2f.h"
924:       Vec              xl,xu
925:       integer          info
926:       PetscInt         i,j,row
927:       PetscInt         xs, xm, ys, ym
928:       PetscScalar      lb,ub

930: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
931: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
932: ! will return an array of doubles referenced by x_array offset by x_index.
933: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
934: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
935:       PetscScalar      xl_v(0:1)
936:       PetscOffset      xl_i


939:       lb = -1.0d300
940:       ub = 1.0d300

942:       if (bmy .lt. 0) bmy = 0
943:       if (bmy .gt. my) bmy = my
944:       if (bmx .lt. 0) bmx = 0
945:       if (bmx .gt. mx) bmx = mx
946:       

948:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
949:      &             PETSC_NULL_INTEGER,info)

951:       call VecSet(xl,lb,info)
952:       call VecSet(xu,ub,info)

954:       call VecGetArray(xl,xl_v,xl_i,info)
955:       

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

959:          do j=ys,ys+ym-1
960:             
961:             row=(j-ys)*xm + (i-xs)

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

967:             endif

969:          enddo
970:       enddo


973:       call VecRestoreArray(xl,xl_v,xl_i,info)
974:       
975:       return
976:       end



980:       
981:       
982: ! ----------------------------------------------------------------------------
983: !
984: !/*
985: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
986: !
987: !     Output Parameter:
988: !.    X - vector for initial guess
989: !
990: !*/

992:       subroutine MSA_InitialPoint(X, info)
993:       implicit none

995: ! mx,my,localX,da,Top,Left,Bottom,Right defined in plate2f.h
996: #include "plate2f.h"
997:       Vec               X
998:       integer           info
999:       PetscInt          start,i,j
1000:       PetscInt          row
1001:       PetscInt          xs,xm,gxs,gxm
1002:       PetscInt          ys,ym,gys,gym
1003:       PetscScalar       zero, np5

1005: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1006: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (integer) x_index, info)
1007: ! will return an array of doubles referenced by x_array offset by x_index.
1008: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
1009: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1010:       PetscScalar   left_v(0:1),right_v(0:1)
1011:       PetscScalar   bottom_v(0:1),top_v(0:1)
1012:       PetscScalar   x_v(0:1)
1013:       PetscOffset   left_i, right_i, top_i
1014:       PetscOffset   bottom_i,x_i
1015:       PetscTruth    flg
1016:       PetscRandom   rctx
1017:       
1018:       zero = 0.0d0
1019:       np5 = -0.5d0

1021:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start",            &
1022:      &                        start,flg,info)

1024:       if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then  ! the zero vector is reasonable
1025:          call VecSet(X,zero,info)

1027:       elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then  ! random start -0.5 < xi < 0.5 
1028:          call PetscRandomCreate(MPI_COMM_WORLD,rctx,info)
1029:          do i=0,start-1
1030:             call VecSetRandom(X,rctx,info)
1031:          enddo

1033:          call PetscRandomDestroy(rctx,info)
1034:          call VecShift(X,np5,info)

1036:       else   ! average of boundary conditions
1037:          
1038: !        Get Local mesh boundaries
1039:          call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,             &
1040:      &                     PETSC_NULL_INTEGER,info) 
1041:          call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,    &
1042:      &                     PETSC_NULL_INTEGER,info)



1046: !        Get pointers to vector data
1047:          call VecGetArray(Top,top_v,top_i,info)
1048:          call VecGetArray(Bottom,bottom_v,bottom_i,info)
1049:          call VecGetArray(Left,left_v,left_i,info)
1050:          call VecGetArray(Right,right_v,right_i,info)
1051:          
1052:          call VecGetArray(localX,x_v,x_i,info)
1053:       
1054: !        Perform local computations 
1055:          do  j=ys,ys+ym-1
1056:             do i=xs,xs+xm-1
1057:                row = (j-gys)*gxm  + (i-gxs)
1058:                x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my        &
1059:      &             + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) +                  &
1060:      &              (i+1)*left_v(left_i+j-ys+1)/mx       +                  &
1061:      &              (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1062:             enddo
1063:          enddo

1065: !        Restore vectors
1066:          call VecRestoreArray(localX,x_v,x_i,info)

1068:          call VecRestoreArray(Left,left_v,left_i,info)
1069:          call VecRestoreArray(Top,top_v,top_i,info)
1070:          call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
1071:          call VecRestoreArray(Right,right_v,right_i,info)

1073:          call DALocalToGlobal(da,localX,INSERT_VALUES,X,info)

1075:       endif

1077:       return
1078:       end