Actual source code: eptorsion2f.F

tao-2.1-p0 2012-07-24
  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: TaoInitialize(); TaoFinalize(); 
 14: !   Routines: TaoCreate(); TaoSetType();
 15: !   Routines: TaoSetInitialVector(); 
 16: !   Routines: TaoSetObjectiveAndGradientRoutine();
 17: !   Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 18: !   Routines: TaoSetMonitor(); TaoSetConvergenceTest()
 19: !   Routines: TaoSolve(); TaoGetSolutionStatus()
 20: !   Routines: TaoGetTerminationReason(); TaoDestroy();

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

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

 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 43: !                   Variable declarations
 44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 45: !
 46: !  See additional variable declarations in the file eptorsion2f.h
 47: !
 48:       PetscErrorCode   ierr           ! used to check for functions returning nonzeros
 49:       Vec              x              ! solution vector
 50:       Mat              H              ! hessian matrix
 51:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
 52:       TaoSolver        tao            ! TaoSolver solver context
 53:       TaoSolverTerminationReason reason
 54:       PetscBool        flg
 55:       PetscInt         iter           ! iteration information
 56:       PetscReal      ff,gnorm,cnorm,xdiff
 57:       PetscInt         i1
 58:       PetscInt         dummy
 59:       

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

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

 67:       i1 = 1

 69: !     Initialize TAO, PETSc  contexts
 70:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 71:       call TaoInitialize(PETSC_NULL_CHARACTER,ierr)

 73: !     Specify default parameters 
 74:       param = 5.0d0
 75:       mx = 10
 76:       my = 10
 77:       Nx = PETSC_DECIDE
 78:       Ny = PETSC_DECIDE

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

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

 93: !     Create vectors
 94:       call DMCreateGlobalVector(dm,x,ierr)
 95:       call DMCreateLocalVector(dm,localX,ierr)

 97: !     Create Hessian
 98:       call DMCreateMatrix(dm,MATAIJ,H,ierr)
 99:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

101: !     The TAO code begins here

103: !     Create TAO solver
104:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
105:       call TaoSetType(tao,"tao_cg",ierr)

107: !     Set routines for function and gradient evaluation

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

114: !     Set initial guess
115:       call FormInitialGuess(x,ierr)
116:       call TaoSetInitialVector(tao,x,ierr)

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

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

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


136: !     SOLVE THE APPLICATION
137:       call TaoSolve(tao,ierr)

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

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

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


158: !     Finalize TAO and PETSc
159:       call PetscFinalize(ierr)
160:       call TaoFinalize(ierr)

162:       end


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

179: ! mx, my are defined in eptorsion2f.h
180: #include "eptorsion2f.h"

182: !  Input/output variables:
183:       Vec              X
184:       PetscErrorCode   ierr

186: !  Local variables:
187:       PetscInt         i, j, k, xe, ye
188:       PetscReal      temp, val, hx, hy
189:       PetscInt         xs, ys, xm, ym
190:       PetscInt         gxm, gym, gxs, gys
191:       PetscInt         i1

193:       i1 = 1
194:       hx = 1.0d0/(mx + 1)
195:       hy = 1.0d0/(my + 1)

197: !  Get corner information
198:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
199:      &                  PETSC_NULL_INTEGER,ierr) 
200:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
201:      &                   gxm,gym,PETSC_NULL_INTEGER,ierr)



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


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

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

248: !  Input/output variables:
249:       TaoSolver        tao
250:       Vec              X, G
251:       PetscReal      f
252:       PetscErrorCode   ierr
253:       PetscInt         dummy

255: !  Declarations for use with local array:


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

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

277:       i1 = 1
278:       0
279:       cdiv3 = param/3.0d0
280:       zero = 0.0d0
281:       p5   = 0.5d0
282:       hx = 1.0d0/(mx + 1)
283:       hy = 1.0d0/(my + 1)
284:       fquad = zero
285:       flin = zero

287: !  Initialize gradient to zero
288:       call VecSet(G,zero,ierr)

290: !  Scatter ghost points to local vector
291:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
292:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)


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

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


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

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

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

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

395: !  Restore vector
396:       call VecRestoreArray(localX,lx_v,lx_i,ierr)

398: !  Assemble gradient vector
399:       call VecAssemblyBegin(G,ierr)
400:       call VecAssemblyEnd(G,ierr)

402: ! Scale the gradient      
403:       area = p5*hx*hy
404:       floc = area *(p5*fquad+flin)
405:       call VecScale(G,area,ierr)


408: !  Sum function contributions from all processes
409:       call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPI_SUM,                   &
410:      &                   PETSC_COMM_WORLD,ierr)
411:       call PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16.0d0,     &
412:      &                    ierr)



416:       return
417:       end




422:       subroutine ComputeHessian(tao, X, H, Hpre, flag, dummy, ierr)
423:       implicit none
424: #include "eptorsion2f.h"      
425:       TaoSolver       tao
426:       Vec             X
427:       Mat             H,Hpre
428:       MatStructure    flag
429:       PetscErrorCode  ierr
430:       PetscInt        dummy

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

440:       i1 = 1

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

448:       do j=ys,ys+ym-1
449:          do i=xs,xs+xm-1
450:             row = (j-gys)*gxm + (i-gxs)

452:             k = 0
453:             if (j .gt. gys) then
454:                v(k) = -1.0d0
455:                col(k) = row-gxm
456:                k = k + 1
457:             endif

459:             if (i .gt. gxs) then
460:                v(k) = -1.0d0
461:                col(k) = row - 1
462:                k = k +1
463:             endif

465:             v(k) = 4.0d0
466:             col(k) = row
467:             k = k + 1

469:             if (i+1 .lt. gxs + gxm) then
470:                v(k) = -1.0d0
471:                col(k) = row + 1
472:                k = k + 1
473:             endif

475:             if (j+1 .lt. gys + gym) then
476:                v(k) = -1.0d0
477:                col(k) = row + gxm
478:                k = k + 1
479:             endif

481:             call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
482:          enddo
483:       enddo

485:       
486: !     Assemble matrix
487:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
488:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)


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

494:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
495:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)


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

500:       0
501:       return 
502:       end
503:       
504:                
505:       
506:       subroutine Monitor(tao, dummy, ierr)
507:       implicit none
508: #include "eptorsion2f.h"      
509:       TaoSolver tao
510:       PetscInt dummy
511:       PetscErrorCode ierr

513:       PetscInt its
514:       PetscReal f,gnorm,cnorm,xdiff
515:       TaoSolverTerminationReason reason
516:       
517:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,             &
518:      &     reason,ierr)
519:       if (mod(its,5) .ne. 0) then
520:          call PetscPrintf(PETSC_COMM_WORLD,"iteration multiple of 5\n",  &
521:      &        ierr)
522:       endif
523:       
524:       0

526:       return
527:       end

529:       subroutine ConvergenceTest(tao, dummy, ierr)
530:       implicit none
531: #include "eptorsion2f.h"      
532:       TaoSolver tao
533:       PetscInt dummy
534:       PetscErrorCode ierr

536:       PetscInt its
537:       PetscReal f,gnorm,cnorm,xdiff
538:       TaoSolverTerminationReason reason
539:       
540:       call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,            &
541:      &     reason,ierr)
542:       if (its .eq. 7) then
543:          call TaoSetTerminationReason(tao,TAO_DIVERGED_MAXITS,          &
544:      &        ierr)
545:       endif
546:       
547:       0

549:       return
550:       end