Actual source code: eptorsion2f.F

petsc-3.7.1 2016-05-15
Report Typos and Errors
  1: !  Program usage: mpiexec -n <proc> eptorsion2f [all TAO options]
  2: !
  3: !  Description:  This example demonstrates use of the TAO package to solve
  4: !  unconstrained minimization problems in parallel.  This example is based
  5: !  on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.
  6: !  The command line options are:
  7: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
  8: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
  9: !    -par <param>, where <param> = angle of twist per unit length
 10: !
 11: !/*T
 12: !   Concepts: TAO^Solving an unconstrained minimization problem
 13: !   Routines: TaoCreate(); TaoSetType();
 14: !   Routines: TaoSetInitialVector();
 15: !   Routines: TaoSetObjectiveAndGradientRoutine();
 16: !   Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 17: !   Routines: TaoSetMonitor(); TaoSetConvergenceTest()
 18: !   Routines: TaoSolve(); TaoGetSolutionStatus()
 19: !   Routines: TaoDestroy();

 21: !   Processors: n
 22: !T*/
 23: !
 24: ! ----------------------------------------------------------------------
 25: !
 26: !  Elastic-plastic torsion problem.
 27: !
 28: !  The elastic plastic torsion problem arises from the deconverged
 29: !  of the stress field on an infinitely long cylindrical bar, which is
 30: !  equivalent to the solution of the following problem:
 31: !     min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
 32: !  where C is the torsion angle per unit length.
 33: !
 34: !  The C version of this code is eptorsion2.c
 35: !
 36: ! ----------------------------------------------------------------------

 38:       implicit none
 39: #include "eptorsion2f.h"

 41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 42: !                   Variable declarations
 43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 44: !
 45: !  See additional variable declarations in the file eptorsion2f.h
 46: !
 47:       PetscErrorCode   ierr           ! used to check for functions returning nonzeros
 48:       Vec              x              ! solution vector
 49:       Mat              H              ! hessian matrix
 50:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
 51:       Tao        tao            ! Tao solver context
 52:       PetscBool        flg
 53:       PetscInt         i1
 54:       PetscInt         dummy


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

 60:       external FormInitialGuess,FormFunctionGradient,ComputeHessian
 61:       external Monitor,ConvergenceTest

 63:       i1 = 1

 65: !     Initialize TAO, PETSc  contexts
 66:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 68: !     Specify default parameters
 69:       param = 5.0
 70:       mx = 10
 71:       my = 10
 72:       Nx = PETSC_DECIDE
 73:       Ny = PETSC_DECIDE

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


 84: !     Set up distributed array and vectors
 85:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,               &
 86:      &     DM_BOUNDARY_NONE,                                             &
 87:      &     DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,        &
 88:      &     PETSC_NULL_INTEGER,dm,ierr)

 90: !     Create vectors
 91:       call DMCreateGlobalVector(dm,x,ierr)
 92:       call DMCreateLocalVector(dm,localX,ierr)

 94: !     Create Hessian
 95:       call DMCreateMatrix(dm,H,ierr)
 96:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

 98: !     The TAO code begins here

100: !     Create TAO solver
101:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
102:       call TaoSetType(tao,TAOCG,ierr)

104: !     Set routines for function and gradient evaluation

106:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
107:      &     FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
108:       call TaoSetHessianRoutine(tao,H,H,ComputeHessian,                 &
109:      &     PETSC_NULL_OBJECT,ierr)

111: !     Set initial guess
112:       call FormInitialGuess(x,ierr)
113:       call TaoSetInitialVector(tao,x,ierr)

115:       call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,   &
116:      &                         '-testmonitor',flg,ierr)
117:       if (flg) then
118:          call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION,      &
119:      &        ierr)
120:       endif

122:       call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,  &
123:      &                         '-testconvergence',flg, ierr)
124:       if (flg) then
125:          call TaoSetConvergenceTest(tao,ConvergenceTest,dummy,          &
126:      &        ierr)
127:       endif

129: !     Check for any TAO command line options
130:       call TaoSetFromOptions(tao,ierr)


133: !     SOLVE THE APPLICATION
134:       call TaoSolve(tao,ierr)

136: !     Free TAO data structures
137:       call TaoDestroy(tao,ierr)


140: !     Free PETSc data structures
141:       call VecDestroy(x,ierr)
142:       call VecDestroy(localX,ierr)
143:       call MatDestroy(H,ierr)
144:       call DMDestroy(dm,ierr)


147: !     Finalize TAO and PETSc
148:       call PetscFinalize(ierr)
149:       end


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

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

169: !  Input/output variables:
170:       Vec              X
171:       PetscErrorCode   ierr

173: !  Local variables:
174:       PetscInt         i, j, k, xe, ye
175:       PetscReal      temp, val, hx, hy
176:       PetscInt         xs, ys, xm, ym
177:       PetscInt         gxm, gym, gxs, gys
178:       PetscInt         i1

180:       i1 = 1
181:       hx = 1.0/(mx + 1)
182:       hy = 1.0/(my + 1)

184: !  Get corner information
185:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
186:      &                  PETSC_NULL_INTEGER,ierr)
187:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
188:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)



192: !  Compute initial guess over locally owned part of mesh
193:       xe = xs+xm
194:       ye = ys+ym
195:       do j=ys,ye-1
196:          temp = min(j+1,my-j)*hy
197:          do i=xs,xe-1
198:             k   = (j-gys)*gxm + i-gxs
199:             val = min((min(i+1,mx-i))*hx,temp)
200:             call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
201:          end do
202:       end do
203:       call VecAssemblyBegin(X,ierr)
204:       call VecAssemblyEnd(X,ierr)
205:       return
206:       end


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

232: ! dm, mx, my, param, localX declared in eptorsion2f.h
233: #include "eptorsion2f.h"

235: !  Input/output variables:
236:       Tao        tao
237:       Vec              X, G
238:       PetscReal      f
239:       PetscErrorCode   ierr
240:       PetscInt         dummy

242: !  Declarations for use with local array:


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

253: !  Local variables:
254:       PetscReal      zero, p5, area, cdiv3
255:       PetscReal      val, flin, fquad,floc
256:       PetscReal      v, vb, vl, vr, vt, dvdx
257:       PetscReal      dvdy, hx, hy
258:       PetscInt         xe, ye, xsm, ysm
259:       PetscInt         xep, yep, i, j, k, ind
260:       PetscInt         xs, ys, xm, ym
261:       PetscInt         gxs, gys, gxm, gym
262:       PetscInt         i1

264:       i1 = 1
265:       0
266:       cdiv3 = param/3.0
267:       zero = 0.0
268:       p5   = 0.5
269:       hx = 1.0/(mx + 1)
270:       hy = 1.0/(my + 1)
271:       fquad = zero
272:       flin = zero

274: !  Initialize gradient to zero
275:       call VecSet(G,zero,ierr)

277: !  Scatter ghost points to local vector
278:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
279:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)


282: !  Get corner information
283:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
284:      &                  PETSC_NULL_INTEGER,ierr)
285:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
286:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)

288: !  Get pointer to vector data.
289:       call VecGetArray(localX,lx_v,lx_i,ierr)


292: !  Set local loop dimensions
293:       xe = xs+xm
294:       ye = ys+ym
295:       if (xs .eq. 0) then
296:          xsm = xs-1
297:       else
298:          xsm = xs
299:       endif
300:       if (ys .eq. 0) then
301:          ysm = ys-1
302:       else
303:          ysm = ys
304:       endif
305:       if (xe .eq. mx) then
306:          xep = xe+1
307:       else
308:          xep = xe
309:       endif
310:       if (ye .eq. my) then
311:          yep = ye+1
312:       else
313:          yep = ye
314:       endif

316: !     Compute local gradient contributions over the lower triangular elements

318:       do j = ysm, ye-1
319:          do i = xsm, xe-1
320:             k  = (j-gys)*gxm + i-gxs
321:             v  = zero
322:             vr = zero
323:             vt = zero
324:             if (i .ge. 0 .and. j .ge. 0)      v = lx_v(lx_i+k)
325:             if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
326:             if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
327:             dvdx = (vr-v)/hx
328:             dvdy = (vt-v)/hy
329:             if (i .ne. -1 .and. j .ne. -1) then
330:                ind = k
331:                val = - dvdx/hx - dvdy/hy - cdiv3
332:                call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
333:             endif
334:             if (i .ne. mx-1 .and. j .ne. -1) then
335:                ind = k+1
336:                val =  dvdx/hx - cdiv3
337:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
338:             endif
339:             if (i .ne. -1 .and. j .ne. my-1) then
340:               ind = k+gxm
341:               val = dvdy/hy - cdiv3
342:               call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
343:             endif
344:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
345:             flin = flin - cdiv3 * (v+vr+vt)
346:          end do
347:       end do

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

351:       do j = ys, yep-1
352:          do i = xs, xep-1
353:             k  = (j-gys)*gxm + i-gxs
354:             vb = zero
355:             vl = zero
356:             v  = zero
357:             if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
358:             if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
359:             if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
360:             dvdx = (v-vl)/hx
361:             dvdy = (v-vb)/hy
362:             if (i .ne. mx .and. j .ne. 0) then
363:                ind = k-gxm
364:                val = - dvdy/hy - cdiv3
365:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
366:             endif
367:             if (i .ne. 0 .and. j .ne. my) then
368:                ind = k-1
369:                val =  - dvdx/hx - cdiv3
370:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
371:             endif
372:             if (i .ne. mx .and. j .ne. my) then
373:                ind = k
374:                val =  dvdx/hx + dvdy/hy - cdiv3
375:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
376:             endif
377:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
378:             flin = flin - cdiv3*(vb + vl + v)
379:          end do
380:       end do

382: !  Restore vector
383:       call VecRestoreArray(localX,lx_v,lx_i,ierr)

385: !  Assemble gradient vector
386:       call VecAssemblyBegin(G,ierr)
387:       call VecAssemblyEnd(G,ierr)

389: ! Scale the gradient
390:       area = p5*hx*hy
391:       floc = area *(p5*fquad+flin)
392:       call VecScale(G,area,ierr)


395: !  Sum function contributions from all processes
396:       call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,                   &
397:      &                   PETSC_COMM_WORLD,ierr)
398:       call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+                        &
399:      &                   16.0d0*(xep-xs)*(yep-ys),ierr)
400:       return
401:       end




406:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
407:       implicit none
408: #include "eptorsion2f.h"
409:       Tao       tao
410:       Vec             X
411:       Mat             H,Hpre
412:       PetscErrorCode  ierr
413:       PetscInt        dummy


416:       PetscInt i,j,k
417:       PetscInt col(0:4),row
418:       PetscInt xs,xm,gxs,gxm
419:       PetscInt ys,ym,gys,gym
420:       PetscReal v(0:4)
421:       PetscInt i1

423:       i1 = 1

425: !     Get local grid boundaries
426:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
427:      &                PETSC_NULL_INTEGER,ierr)
428:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,      &
429:      &                PETSC_NULL_INTEGER,ierr)

431:       do j=ys,ys+ym-1
432:          do i=xs,xs+xm-1
433:             row = (j-gys)*gxm + (i-gxs)

435:             k = 0
436:             if (j .gt. gys) then
437:                v(k) = -1.0
438:                col(k) = row-gxm
439:                k = k + 1
440:             endif

442:             if (i .gt. gxs) then
443:                v(k) = -1.0
444:                col(k) = row - 1
445:                k = k +1
446:             endif

448:             v(k) = 4.0
449:             col(k) = row
450:             k = k + 1

452:             if (i+1 .lt. gxs + gxm) then
453:                v(k) = -1.0
454:                col(k) = row + 1
455:                k = k + 1
456:             endif

458:             if (j+1 .lt. gys + gym) then
459:                v(k) = -1.0
460:                col(k) = row + gxm
461:                k = k + 1
462:             endif

464:             call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
465:          enddo
466:       enddo


469: !     Assemble matrix
470:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
471:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)


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

477:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
478:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)


481:       call PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)

483:       0
484:       return
485:       end



489:       subroutine Monitor(tao, dummy, ierr)
490:       implicit none
491: #include "eptorsion2f.h"
492:       Tao tao
493:       PetscInt dummy
494:       PetscErrorCode ierr

496:       PetscInt its
497:       PetscReal f,gnorm,cnorm,xdiff
498:       TaoConvergedReason reason

500:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,             &
501:      &     reason,ierr)
502:       if (mod(its,5) .ne. 0) then
503:          call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',  &
504:      &        ierr)
505:       endif

507:       0

509:       return
510:       end

512:       subroutine ConvergenceTest(tao, dummy, ierr)
513:       implicit none
514: #include "eptorsion2f.h"
515:       Tao tao
516:       PetscInt dummy
517:       PetscErrorCode ierr

519:       PetscInt its
520:       PetscReal f,gnorm,cnorm,xdiff
521:       TaoConvergedReason reason

523:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,            &
524:      &     reason,ierr)
525:       if (its .eq. 7) then
526:        call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
527:       endif

529:       0

531:       return
532:       end