Actual source code: eptorsion2f.F

  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:       module mymodule
 39: #include "petsc/finclude/petsctao.h"
 40:       use petscdmda
 41:       use petsctao
 42:       implicit none

 44:       Vec              localX
 45:       DM               dm
 46:       PetscReal      param
 47:       PetscInt         mx, my
 48:       end module

 50:       use mymodule
 51:       implicit none
 52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53: !                   Variable declarations
 54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55: !
 56: !  See additional variable declarations in the file eptorsion2f.h
 57: !
 58:       PetscErrorCode   ierr           ! used to check for functions returning nonzeros
 59:       Vec              x              ! solution vector
 60:       Mat              H              ! hessian matrix
 61:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
 62:       Tao        tao            ! Tao solver context
 63:       PetscBool        flg
 64:       PetscInt         i1
 65:       PetscInt         dummy

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

 70:       external FormInitialGuess,FormFunctionGradient,ComputeHessian
 71:       external Monitor,ConvergenceTest

 73:       i1 = 1

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

 82: !     Specify default parameters
 83:       param = 5.0
 84:       mx = 10
 85:       my = 10
 86:       Nx = PETSC_DECIDE
 87:       Ny = PETSC_DECIDE

 89: !     Check for any command line arguments that might override defaults
 90:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,    &
 91:      &                        '-mx',mx,flg,ierr)
 92:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,    &
 93:      &                        '-my',my,flg,ierr)
 94:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
 95:      &                         '-par',param,flg,ierr)

 97: !     Set up distributed array and vectors
 98:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,               &
 99:      &     DM_BOUNDARY_NONE,                                             &
100:      &     DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,        &
101:      &     PETSC_NULL_INTEGER,dm,ierr)
102:       call DMSetFromOptions(dm,ierr)
103:       call DMSetUp(dm,ierr)

105: !     Create vectors
106:       call DMCreateGlobalVector(dm,x,ierr)
107:       call DMCreateLocalVector(dm,localX,ierr)

109: !     Create Hessian
110:       call DMCreateMatrix(dm,H,ierr)
111:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

113: !     The TAO code begins here

115: !     Create TAO solver
116:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
117:       call TaoSetType(tao,TAOCG,ierr)

119: !     Set routines for function and gradient evaluation

121:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
122:      &     FormFunctionGradient,0,ierr)
123:       call TaoSetHessianRoutine(tao,H,H,ComputeHessian,                 &
124:      &     0,ierr)

126: !     Set initial guess
127:       call FormInitialGuess(x,ierr)
128:       call TaoSetInitialVector(tao,x,ierr)

130:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
131:      &                         '-testmonitor',flg,ierr)
132:       if (flg) then
133:          call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION,      &
134:      &        ierr)
135:       endif

137:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
138:      &                         '-testconvergence',flg, ierr)
139:       if (flg) then
140:          call TaoSetConvergenceTest(tao,ConvergenceTest,dummy,          &
141:      &        ierr)
142:       endif

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

147: !     SOLVE THE APPLICATION
148:       call TaoSolve(tao,ierr)

150: !     Free TAO data structures
151:       call TaoDestroy(tao,ierr)

153: !     Free PETSc data structures
154:       call VecDestroy(x,ierr)
155:       call VecDestroy(localX,ierr)
156:       call MatDestroy(H,ierr)
157:       call DMDestroy(dm,ierr)

159: !     Finalize TAO and PETSc
160:       call PetscFinalize(ierr)
161:       end

163: ! ---------------------------------------------------------------------
164: !
165: !   FormInitialGuess - Computes an initial approximation to the solution.
166: !
167: !   Input Parameters:
168: !   X    - vector
169: !
170: !   Output Parameters:
171: !   X    - vector
172: !   ierr - error code
173: !
174:       subroutine FormInitialGuess(X,ierr)
175:       use mymodule
176:       implicit none

178: !  Input/output variables:
179:       Vec              X
180:       PetscErrorCode   ierr

182: !  Local variables:
183:       PetscInt         i, j, k, xe, ye
184:       PetscReal      temp, val, hx, hy
185:       PetscInt         xs, ys, xm, ym
186:       PetscInt         gxm, gym, gxs, gys
187:       PetscInt         i1

189:       i1 = 1
190:       hx = 1.0/real(mx + 1)
191:       hy = 1.0/real(my + 1)

193: !  Get corner information
194:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
195:      &                  PETSC_NULL_INTEGER,ierr)
196:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
197:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)

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

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

239: !  Input/output variables:
240:       Tao        tao
241:       Vec              X, G
242:       PetscReal      f
243:       PetscErrorCode   ierr
244:       PetscInt         dummy

246: !  Declarations for use with local array:

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

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

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

277: !  Initialize gradient to zero
278:       call VecSet(G,zero,ierr)

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

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

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

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

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

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

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

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

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

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

390: ! Scale the gradient
391:       area = p5*hx*hy
392:       floc = area *(p5*fquad+flin)
393:       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

403:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
404:       use mymodule
405:       implicit none

407:       Tao       tao
408:       Vec             X
409:       Mat             H,Hpre
410:       PetscErrorCode  ierr
411:       PetscInt        dummy

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

420:       i1 = 1

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

428:       do j=ys,ys+ym-1
429:          do i=xs,xs+xm-1
430:             row = (j-gys)*gxm + (i-gxs)

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

439:             if (i .gt. gxs) then
440:                v(k) = -1.0
441:                col(k) = row - 1
442:                k = k +1
443:             endif

445:             v(k) = 4.0
446:             col(k) = row
447:             k = k + 1

449:             if (i+1 .lt. gxs + gxm) then
450:                v(k) = -1.0
451:                col(k) = row + 1
452:                k = k + 1
453:             endif

455:             if (j+1 .lt. gys + gym) then
456:                v(k) = -1.0
457:                col(k) = row + gxm
458:                k = k + 1
459:             endif

461:             call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
462:          enddo
463:       enddo

465: !     Assemble matrix
466:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
467:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)

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

472:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
473:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

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

477:       0
478:       return
479:       end

481:       subroutine Monitor(tao, dummy, ierr)
482:       use mymodule
483:       implicit none

485:       Tao tao
486:       PetscInt dummy
487:       PetscErrorCode ierr

489:       PetscInt its
490:       PetscReal f,gnorm,cnorm,xdiff
491:       TaoConvergedReason reason

493:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,             &
494:      &     reason,ierr)
495:       if (mod(its,5) .ne. 0) then
496:          call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',  &
497:      &        ierr)
498:       endif

500:       0

502:       return
503:       end

505:       subroutine ConvergenceTest(tao, dummy, ierr)
506:       use mymodule
507:       implicit none

509:       Tao tao
510:       PetscInt dummy
511:       PetscErrorCode ierr

513:       PetscInt its
514:       PetscReal f,gnorm,cnorm,xdiff
515:       TaoConvergedReason reason

517:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,            &
518:      &     reason,ierr)
519:       if (its .eq. 7) then
520:        call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
521:       endif

523:       0

525:       return
526:       end

528: !/*TEST
529: !
530: !   build:
531: !      requires: !complex
532: !
533: !   test:
534: !      args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
535: !
536: !   test:
537: !      suffix: 2
538: !      nsize: 2
539: !      args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
540: !
541: !   test:
542: !      suffix: 3
543: !      args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
544: !TEST*/