Actual source code: eptorsion2f.F

petsc-master 2014-11-25
Report Typos and Errors
  1: !  Program usage: mpirun -np <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: TaoGetConvergedReason(); 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:       TaoConvergedReason reason
 53:       PetscBool        flg
 54:       PetscInt         iter           ! iteration information
 55:       PetscReal      ff,gnorm,cnorm,xdiff
 56:       PetscInt         i1
 57:       PetscInt         dummy


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

 63:       external FormInitialGuess,FormFunctionGradient,ComputeHessian
 64:       external Monitor,ConvergenceTest

 66:       i1 = 1

 68: !     Initialize TAO, PETSc  contexts
 69:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 71: !     Specify default parameters
 72:       param = 5.0
 73:       mx = 10
 74:       my = 10
 75:       Nx = PETSC_DECIDE
 76:       Ny = PETSC_DECIDE

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


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

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

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

 99: !     The TAO code begins here

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

105: !     Set routines for function and gradient evaluation

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

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

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

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

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


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

137: !     Get information on termination
138:       call TaoGetConvergedReason(tao,reason,ierr)
139:       if (reason .lt. 0) then
140:           call PetscPrintF(PETSC_COMM_WORLD,                              &
141:      &        'TAO did not terminate successfully\n',ierr)
142:       endif


145: !     Free TAO data structures
146:       call TaoDestroy(tao,ierr)


149: !     Free PETSc data structures
150:       call VecDestroy(x,ierr)
151:       call VecDestroy(localX,ierr)
152:       call MatDestroy(H,ierr)
153:       call DMDestroy(dm,ierr)


156: !     Finalize TAO and PETSc
157:       call PetscFinalize(ierr)
158:       end


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

175: ! mx, my are defined in eptorsion2f.h
176: #include "eptorsion2f.h"

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/(mx + 1)
191:       hy = 1.0/(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)



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


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

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

244: !  Input/output variables:
245:       Tao        tao
246:       Vec              X, G
247:       PetscReal      f
248:       PetscErrorCode   ierr
249:       PetscInt         dummy

251: !  Declarations for use with local array:


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

262: !  Local variables:
263:       PetscReal      zero, p5, area, cdiv3
264:       PetscReal      val, flin, fquad,floc
265:       PetscReal      v, vb, vl, vr, vt, dvdx
266:       PetscReal      dvdy, hx, hy
267:       PetscInt         xe, ye, xsm, ysm
268:       PetscInt         xep, yep, i, j, k, ind
269:       PetscInt         xs, ys, xm, ym
270:       PetscInt         gxs, gys, gxm, gym
271:       PetscInt         i1

273:       i1 = 1
274:       0
275:       cdiv3 = param/3.0
276:       zero = 0.0
277:       p5   = 0.5
278:       hx = 1.0/(mx + 1)
279:       hy = 1.0/(my + 1)
280:       fquad = zero
281:       flin = zero

283: !  Initialize gradient to zero
284:       call VecSet(G,zero,ierr)

286: !  Scatter ghost points to local vector
287:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
288:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)


291: !  Get corner information
292:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
293:      &                  PETSC_NULL_INTEGER,ierr)
294:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
295:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)

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


301: !  Set local loop dimensions
302:       xe = xs+xm
303:       ye = ys+ym
304:       if (xs .eq. 0) then
305:          xsm = xs-1
306:       else
307:          xsm = xs
308:       endif
309:       if (ys .eq. 0) then
310:          ysm = ys-1
311:       else
312:          ysm = ys
313:       endif
314:       if (xe .eq. mx) then
315:          xep = xe+1
316:       else
317:          xep = xe
318:       endif
319:       if (ye .eq. my) then
320:          yep = ye+1
321:       else
322:          yep = ye
323:       endif

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

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

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

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

391: !  Restore vector
392:       call VecRestoreArray(localX,lx_v,lx_i,ierr)

394: !  Assemble gradient vector
395:       call VecAssemblyBegin(G,ierr)
396:       call VecAssemblyEnd(G,ierr)

398: ! Scale the gradient
399:       area = p5*hx*hy
400:       floc = area *(p5*fquad+flin)
401:       call VecScale(G,area,ierr)


404: !  Sum function contributions from all processes
405:       call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,                   &
406:      &                   PETSC_COMM_WORLD,ierr)
407:       call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+                        &
408:      &                   16.0d0*(xep-xs)*(yep-ys),ierr)
409:       return
410:       end




415:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
416:       implicit none
417: #include "eptorsion2f.h"
418:       Tao       tao
419:       Vec             X
420:       Mat             H,Hpre
421:       PetscErrorCode  ierr
422:       PetscInt        dummy


425:       PetscInt i,j,k
426:       PetscInt col(0:4),row
427:       PetscInt xs,xm,gxs,gxm
428:       PetscInt ys,ym,gys,gym
429:       PetscReal v(0:4)
430:       PetscInt i1

432:       i1 = 1

434: !     Get local grid boundaries
435:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
436:      &                PETSC_NULL_INTEGER,ierr)
437:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,      &
438:      &                PETSC_NULL_INTEGER,ierr)

440:       do j=ys,ys+ym-1
441:          do i=xs,xs+xm-1
442:             row = (j-gys)*gxm + (i-gxs)

444:             k = 0
445:             if (j .gt. gys) then
446:                v(k) = -1.0
447:                col(k) = row-gxm
448:                k = k + 1
449:             endif

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

457:             v(k) = 4.0
458:             col(k) = row
459:             k = k + 1

461:             if (i+1 .lt. gxs + gxm) then
462:                v(k) = -1.0
463:                col(k) = row + 1
464:                k = k + 1
465:             endif

467:             if (j+1 .lt. gys + gym) then
468:                v(k) = -1.0
469:                col(k) = row + gxm
470:                k = k + 1
471:             endif

473:             call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
474:          enddo
475:       enddo


478: !     Assemble matrix
479:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
480:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)


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

486:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
487:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)


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

492:       0
493:       return
494:       end



498:       subroutine Monitor(tao, dummy, ierr)
499:       implicit none
500: #include "eptorsion2f.h"
501:       Tao tao
502:       PetscInt dummy
503:       PetscErrorCode ierr

505:       PetscInt its
506:       PetscReal f,gnorm,cnorm,xdiff
507:       TaoConvergedReason reason

509:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,             &
510:      &     reason,ierr)
511:       if (mod(its,5) .ne. 0) then
512:          call PetscPrintf(PETSC_COMM_WORLD,"iteration multiple of 5\n",  &
513:      &        ierr)
514:       endif

516:       0

518:       return
519:       end

521:       subroutine ConvergenceTest(tao, dummy, ierr)
522:       implicit none
523: #include "eptorsion2f.h"
524:       Tao tao
525:       PetscInt dummy
526:       PetscErrorCode ierr

528:       PetscInt its
529:       PetscReal f,gnorm,cnorm,xdiff
530:       TaoConvergedReason reason

532:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,            &
533:      &     reason,ierr)
534:       if (its .eq. 7) then
535:        call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
536:       endif

538:       0

540:       return
541:       end