Actual source code: eptorsion2f.F

petsc-master 2018-05-25
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:       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: #if defined(PETSC_USING_F90) && !defined(PETSC_USE_FORTRANKIND)
 53:       external PETSC_NULL_FUNCTION
 54: #endif
 55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56: !                   Variable declarations
 57: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 58: !
 59: !  See additional variable declarations in the file eptorsion2f.h
 60: !
 61:       PetscErrorCode   ierr           ! used to check for functions returning nonzeros
 62:       Vec              x              ! solution vector
 63:       Mat              H              ! hessian matrix
 64:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
 65:       Tao        tao            ! Tao solver context
 66:       PetscBool        flg
 67:       PetscInt         i1
 68:       PetscInt         dummy


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

 74:       external FormInitialGuess,FormFunctionGradient,ComputeHessian
 75:       external Monitor,ConvergenceTest

 77:       i1 = 1

 79: !     Initialize TAO, PETSc  contexts
 80:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 81:       if (ierr .ne. 0) then
 82:          print*,'Unable to initialize PETSc'
 83:          stop
 84:       endif

 86: !     Specify default parameters
 87:       param = 5.0
 88:       mx = 10
 89:       my = 10
 90:       Nx = PETSC_DECIDE
 91:       Ny = PETSC_DECIDE

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


102: !     Set up distributed array and vectors
103:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,               &
104:      &     DM_BOUNDARY_NONE,                                             &
105:      &     DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,        &
106:      &     PETSC_NULL_INTEGER,dm,ierr)
107:       call DMSetFromOptions(dm,ierr)
108:       call DMSetUp(dm,ierr)

110: !     Create vectors
111:       call DMCreateGlobalVector(dm,x,ierr)
112:       call DMCreateLocalVector(dm,localX,ierr)

114: !     Create Hessian
115:       call DMCreateMatrix(dm,H,ierr)
116:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

118: !     The TAO code begins here

120: !     Create TAO solver
121:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
122:       call TaoSetType(tao,TAOCG,ierr)

124: !     Set routines for function and gradient evaluation

126:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
127:      &     FormFunctionGradient,0,ierr)
128:       call TaoSetHessianRoutine(tao,H,H,ComputeHessian,                 &
129:      &     0,ierr)

131: !     Set initial guess
132:       call FormInitialGuess(x,ierr)
133:       call TaoSetInitialVector(tao,x,ierr)

135:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
136:      &                         '-testmonitor',flg,ierr)
137:       if (flg) then
138:          call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION,      &
139:      &        ierr)
140:       endif

142:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
143:      &                         '-testconvergence',flg, ierr)
144:       if (flg) then
145:          call TaoSetConvergenceTest(tao,ConvergenceTest,dummy,          &
146:      &        ierr)
147:       endif

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


153: !     SOLVE THE APPLICATION
154:       call TaoSolve(tao,ierr)

156: !     Free TAO data structures
157:       call TaoDestroy(tao,ierr)


160: !     Free PETSc data structures
161:       call VecDestroy(x,ierr)
162:       call VecDestroy(localX,ierr)
163:       call MatDestroy(H,ierr)
164:       call DMDestroy(dm,ierr)


167: !     Finalize TAO and PETSc
168:       call PetscFinalize(ierr)
169:       end


172: ! ---------------------------------------------------------------------
173: !
174: !   FormInitialGuess - Computes an initial approximation to the solution.
175: !
176: !   Input Parameters:
177: !   X    - vector
178: !
179: !   Output Parameters:
180: !   X    - vector
181: !   ierr - error code
182: !
183:       subroutine FormInitialGuess(X,ierr)
184:       use mymodule
185:       implicit none

187: !  Input/output variables:
188:       Vec              X
189:       PetscErrorCode   ierr

191: !  Local variables:
192:       PetscInt         i, j, k, xe, ye
193:       PetscReal      temp, val, hx, hy
194:       PetscInt         xs, ys, xm, ym
195:       PetscInt         gxm, gym, gxs, gys
196:       PetscInt         i1

198:       i1 = 1
199:       hx = 1.0/real(mx + 1)
200:       hy = 1.0/real(my + 1)

202: !  Get corner information
203:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
204:      &                  PETSC_NULL_INTEGER,ierr)
205:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
206:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)



210: !  Compute initial guess over locally owned part of mesh
211:       xe = xs+xm
212:       ye = ys+ym
213:       do j=ys,ye-1
214:          temp = min(j+1,my-j)*hy
215:          do i=xs,xe-1
216:             k   = (j-gys)*gxm + i-gxs
217:             val = min((min(i+1,mx-i))*hx,temp)
218:             call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
219:          end do
220:       end do
221:       call VecAssemblyBegin(X,ierr)
222:       call VecAssemblyEnd(X,ierr)
223:       return
224:       end


227: ! ---------------------------------------------------------------------
228: !
229: !  FormFunctionGradient - Evaluates gradient G(X).
230: !
231: !  Input Parameters:
232: !  tao   - the Tao context
233: !  X     - input vector
234: !  dummy - optional user-defined context (not used here)
235: !
236: !  Output Parameters:
237: !  f     - the function value at X
238: !  G     - vector containing the newly evaluated gradient
239: !  ierr  - error code
240: !
241: !  Notes:
242: !  This routine serves as a wrapper for the lower-level routine
243: !  "ApplicationGradient", where the actual computations are
244: !  done using the standard Fortran style of treating the local
245: !  input vector data as an array over the local mesh.
246: !
247:       subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
248:       use mymodule
249:       implicit none

251: !  Input/output variables:
252:       Tao        tao
253:       Vec              X, G
254:       PetscReal      f
255:       PetscErrorCode   ierr
256:       PetscInt         dummy

258: !  Declarations for use with local array:


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

269: !  Local variables:
270:       PetscReal      zero, p5, area, cdiv3
271:       PetscReal      val, flin, fquad,floc
272:       PetscReal      v, vb, vl, vr, vt, dvdx
273:       PetscReal      dvdy, hx, hy
274:       PetscInt         xe, ye, xsm, ysm
275:       PetscInt         xep, yep, i, j, k, ind
276:       PetscInt         xs, ys, xm, ym
277:       PetscInt         gxs, gys, gxm, gym
278:       PetscInt         i1

280:       i1 = 1
281:       0
282:       cdiv3 = param/3.0
283:       zero = 0.0
284:       p5   = 0.5
285:       hx = 1.0/real(mx + 1)
286:       hy = 1.0/real(my + 1)
287:       fquad = zero
288:       flin = zero

290: !  Initialize gradient to zero
291:       call VecSet(G,zero,ierr)

293: !  Scatter ghost points to local vector
294:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
295:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)


298: !  Get corner information
299:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
300:      &                  PETSC_NULL_INTEGER,ierr)
301:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
302:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)

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


308: !  Set local loop dimensions
309:       xe = xs+xm
310:       ye = ys+ym
311:       if (xs .eq. 0) then
312:          xsm = xs-1
313:       else
314:          xsm = xs
315:       endif
316:       if (ys .eq. 0) then
317:          ysm = ys-1
318:       else
319:          ysm = ys
320:       endif
321:       if (xe .eq. mx) then
322:          xep = xe+1
323:       else
324:          xep = xe
325:       endif
326:       if (ye .eq. my) then
327:          yep = ye+1
328:       else
329:          yep = ye
330:       endif

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

334:       do j = ysm, ye-1
335:          do i = xsm, xe-1
336:             k  = (j-gys)*gxm + i-gxs
337:             v  = zero
338:             vr = zero
339:             vt = zero
340:             if (i .ge. 0 .and. j .ge. 0)      v = lx_v(lx_i+k)
341:             if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
342:             if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
343:             dvdx = (vr-v)/hx
344:             dvdy = (vt-v)/hy
345:             if (i .ne. -1 .and. j .ne. -1) then
346:                ind = k
347:                val = - dvdx/hx - dvdy/hy - cdiv3
348:                call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
349:             endif
350:             if (i .ne. mx-1 .and. j .ne. -1) then
351:                ind = k+1
352:                val =  dvdx/hx - cdiv3
353:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
354:             endif
355:             if (i .ne. -1 .and. j .ne. my-1) then
356:               ind = k+gxm
357:               val = dvdy/hy - cdiv3
358:               call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
359:             endif
360:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
361:             flin = flin - cdiv3 * (v+vr+vt)
362:          end do
363:       end do

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

367:       do j = ys, yep-1
368:          do i = xs, xep-1
369:             k  = (j-gys)*gxm + i-gxs
370:             vb = zero
371:             vl = zero
372:             v  = zero
373:             if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
374:             if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
375:             if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
376:             dvdx = (v-vl)/hx
377:             dvdy = (v-vb)/hy
378:             if (i .ne. mx .and. j .ne. 0) then
379:                ind = k-gxm
380:                val = - dvdy/hy - cdiv3
381:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
382:             endif
383:             if (i .ne. 0 .and. j .ne. my) then
384:                ind = k-1
385:                val =  - dvdx/hx - cdiv3
386:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
387:             endif
388:             if (i .ne. mx .and. j .ne. my) then
389:                ind = k
390:                val =  dvdx/hx + dvdy/hy - cdiv3
391:                call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
392:             endif
393:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
394:             flin = flin - cdiv3*(vb + vl + v)
395:          end do
396:       end do

398: !  Restore vector
399:       call VecRestoreArray(localX,lx_v,lx_i,ierr)

401: !  Assemble gradient vector
402:       call VecAssemblyBegin(G,ierr)
403:       call VecAssemblyEnd(G,ierr)

405: ! Scale the gradient
406:       area = p5*hx*hy
407:       floc = area *(p5*fquad+flin)
408:       call VecScale(G,area,ierr)


411: !  Sum function contributions from all processes
412:       call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,                   &
413:      &                   PETSC_COMM_WORLD,ierr)
414:       call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+                        &
415:      &                   16.0d0*(xep-xs)*(yep-ys),ierr)
416:       return
417:       end

419:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
420:       use mymodule
421:       implicit none

423:       Tao       tao
424:       Vec             X
425:       Mat             H,Hpre
426:       PetscErrorCode  ierr
427:       PetscInt        dummy


430:       PetscInt i,j,k
431:       PetscInt col(0:4),row
432:       PetscInt xs,xm,gxs,gxm
433:       PetscInt ys,ym,gys,gym
434:       PetscReal v(0:4)
435:       PetscInt i1

437:       i1 = 1

439: !     Get local grid boundaries
440:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
441:      &                PETSC_NULL_INTEGER,ierr)
442:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,      &
443:      &                PETSC_NULL_INTEGER,ierr)

445:       do j=ys,ys+ym-1
446:          do i=xs,xs+xm-1
447:             row = (j-gys)*gxm + (i-gxs)

449:             k = 0
450:             if (j .gt. gys) then
451:                v(k) = -1.0
452:                col(k) = row-gxm
453:                k = k + 1
454:             endif

456:             if (i .gt. gxs) then
457:                v(k) = -1.0
458:                col(k) = row - 1
459:                k = k +1
460:             endif

462:             v(k) = 4.0
463:             col(k) = row
464:             k = k + 1

466:             if (i+1 .lt. gxs + gxm) then
467:                v(k) = -1.0
468:                col(k) = row + 1
469:                k = k + 1
470:             endif

472:             if (j+1 .lt. gys + gym) then
473:                v(k) = -1.0
474:                col(k) = row + gxm
475:                k = k + 1
476:             endif

478:             call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
479:          enddo
480:       enddo


483: !     Assemble matrix
484:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
485:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)


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

491:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
492:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)


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

497:       0
498:       return
499:       end



503:       subroutine Monitor(tao, dummy, ierr)
504:       use mymodule
505:       implicit none

507:       Tao tao
508:       PetscInt dummy
509:       PetscErrorCode ierr

511:       PetscInt its
512:       PetscReal f,gnorm,cnorm,xdiff
513:       TaoConvergedReason reason

515:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,             &
516:      &     reason,ierr)
517:       if (mod(its,5) .ne. 0) then
518:          call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',  &
519:      &        ierr)
520:       endif

522:       0

524:       return
525:       end

527:       subroutine ConvergenceTest(tao, dummy, ierr)
528:       use mymodule
529:       implicit none

531:       Tao tao
532:       PetscInt dummy
533:       PetscErrorCode ierr

535:       PetscInt its
536:       PetscReal f,gnorm,cnorm,xdiff
537:       TaoConvergedReason reason

539:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,            &
540:      &     reason,ierr)
541:       if (its .eq. 7) then
542:        call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
543:       endif

545:       0

547:       return
548:       end

550: !/*TEST
551: !
552: !   build:
553: !      requires: !complex
554: !
555: !   test:
556: !      args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
557: !
558: !   test:
559: !      suffix: 2
560: !      nsize: 2
561: !      args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
562: !
563: !TEST*/