Actual source code: eptorsion2f.F

  1: ! "$Id$";
  2: !
  3: !  Program usage: mpirun -np <proc> eptorsion2f [all TAO options] 
  4: !
  5: !  Description:  This example demonstrates use of the TAO package to solve
  6: !  unconstrained minimization problems in parallel.  This example is based 
  7: !  on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.
  8: !  The command line options are:
  9: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
 10: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
 11: !    -par <param>, where <param> = angle of twist per unit length
 12: !
 13: !/*T
 14: !   Concepts: TAO - Solving an unconstrained minimization problem
 15: !   Routines: TaoInitialize(); TaoFinalize(); 
 16: !   Routines: TaoCreate(); TaoDestroy();
 17: !   Routines: TaoApplicationCreate(); TaoAppDestroy();
 18: !   Routines: TaoAppSetObjectiveAndGradientRoutine();
 19: !   Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
 20: !   Routines: TaoSetOptions();
 21: !   Routines: TaoAppSetInitialSolutionVec(); TaoSolveApplication(); TaoDestroy();
 22: !   Routines: TaoGetSolutionStatus();
 23: !   Processors: n
 24: !T*/
 25: !
 26: ! ---------------------------------------------------------------------- 
 27: !
 28: !  Elastic-plastic torsion problem.  
 29: !
 30: !  The elastic plastic torsion problem arises from the determination 
 31: !  of the stress field on an infinitely long cylindrical bar, which is
 32: !  equivalent to the solution of the following problem:
 33: !     min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
 34: !  where C is the torsion angle per unit length.
 35: !
 36: !  The C version of this code is eptorsion2.c
 37: !
 38: ! ---------------------------------------------------------------------- 

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

 43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 44: !                   Variable declarations
 45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 46: !
 47: !  See additional variable declarations in the file eptorsion2f.h
 48: !
 49:       integer          info           ! used to check for functions returning nonzeros
 50:       Vec              x              ! solution vector
 51:       Mat              H              ! hessian matrix
 52:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
 53:       TAO_SOLVER       tao            ! TAO_SOLVER solver context
 54:       TAO_APPLICATION  torsionapp     ! TAO application context (PETSc)
 55:       TaoTerminateReason reason
 56:       PetscTruth       flg
 57:       PetscInt         iter           ! iteration information
 58:       PetscScalar      ff,gnorm,cnorm,xdiff
 59:       PetscInt         i1
 60:       

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

 65:       external FormInitialGuess,FormFunctionGradient,ComputeHessian

 67:       i1 = 1

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

 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,info)
 82:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,info)
 83:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-par",              &
 84:      &                         param,flg,info)

 86:       
 87: !     Set up distributed array and vectors
 88:       call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,      &
 89:      &     mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,      &
 90:      &     da,info)

 92: !     Create vectors
 93:       call DACreateGlobalVector(da,x,info)
 94:       call DACreateLocalVector(da,localX,info)

 96: !     Create Hessian
 97:       call DAGetMatrix(da,MATAIJ,H,info)
 98:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)

100: !     The TAO code begins here

102: !     Create TAO solver
103:       call TaoCreate(MPI_COMM_WORLD,'tao_cg',tao,info)
104:       call TaoApplicationCreate(MPI_COMM_WORLD,torsionapp,info)

106: !     Set routines for function and gradient evaluation

108: !     TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
109:       call TaoAppSetObjectiveAndGradientRo(torsionapp,                       &
110:      &     FormFunctionGradient,TAO_NULL_OBJECT,info)
111:       call TaoAppSetHessianMat(torsionapp,H,H,info)
112:       call TaoAppSetHessianRoutine(torsionapp,ComputeHessian,                &
113:      &    TAO_NULL_OBJECT,info)

115: !     Set initial guess
116:       call FormInitialGuess(x,info)
117:       call TaoAppSetInitialSolutionVec(torsionapp,x,info)

119: !     Check for any TAO command line options
120:       call TaoSetOptions(torsionapp, tao,info)


123: !     SOLVE THE APPLICATION
124:       call TaoSolveApplication(torsionapp,tao,info)

126: !     Get information on termination
127:       call TaoGetSolutionStatus(tao,iter,ff,gnorm,cnorm,xdiff,           &
128:      &                          reason,info)
129:       if (reason .lt. 0) then
130:          print *,'TAO did not terminate successfully'
131:       endif

133:       
134: !     Free TAO data structures
135:       call TaoDestroy(tao,info)
136:       call TaoAppDestroy(torsionapp,info);  

138:   
139: !     Free PETSc data structures 
140:       call VecDestroy(x,info)
141:       call VecDestroy(localX,info)
142:       call MatDestroy(H,info)
143:       call DADestroy(da,info)


146: !     Finalize TAO and PETSc
147:       call PetscFinalize(info)
148:       call TaoFinalize(info)

150:       end


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

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

170: !  Input/output variables:
171:       Vec              X
172:       integer          info

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

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

185: !  Get corner information
186:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
187:      &                  PETSC_NULL_INTEGER,info) 
188:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,              &
189:      &                   gxm,gym,PETSC_NULL_INTEGER,info)



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


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

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

236: !  Input/output variables:
237:       TAO_APPLICATION  taoapp
238:       Vec              X, G
239:       PetscScalar      f
240:       integer          info
241:       PetscInt         dummy

243: !  Declarations for use with local array:


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

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

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

275: !  Initialize gradient to zero
276:       call VecSet(G,zero,info)

278: !  Scatter ghost points to local vector
279:       call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
280:       call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)


283: !  Get corner information
284:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
285:      &                  PETSC_NULL_INTEGER,info)
286:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,              &
287:      &                   gxm,gym,PETSC_NULL_INTEGER,info)

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


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
318:      
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,info)
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,info)
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,info)
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,info)
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,info)
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,info)
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,info)
385:       if (info .ne. 0) print *,'VecRestoreArray'
386: !  Assemble gradient vector
387:       call VecAssemblyBegin(G,info)
388:       call VecAssemblyEnd(G,info)

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


396: !  Sum function contributions from all processes
397:       call MPI_Allreduce(floc,f,1,MPI_DOUBLE_PRECISION,MPI_SUM,           &
398:      &                   MPI_COMM_WORLD,info)
399:       if (info .ne. 0) print *,'MPI_Allreduce'
400:       call PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16,       &
401:      &                    info)



405:       return
406:       end




411:       subroutine ComputeHessian(taoapp, X, H, Hpre, flag, dummy, info)
412:       implicit none
413: #include "eptorsion2f.h"      
414:       TAO_APPLICATION taoapp
415:       Vec             X
416:       Mat             H,Hpre
417:       MatStructure    flag
418:       integer         info
419:       PetscInt        dummy

421:       
422:       PetscInt i,j,k
423:       PetscInt col(0:4),row
424:       PetscInt xs,xm,gxs,gxm
425:       PetscInt ys,ym,gys,gym
426:       PetscScalar v(0:4)
427:       PetscInt i1

429:       i1 = 1

431: !     Get local grid boundaries
432:       call DAGetCorners(da,xs,ys,TAO_NULL_INTEGER,xm,ym,                &
433:      &                TAO_NULL_INTEGER,info)
434:       call DAGetGhostCorners(da,gxs,gys,TAO_NULL_INTEGER,gxm,gym,        &
435:      &                TAO_NULL_INTEGER,info)

437:       do j=ys,ys+ym-1
438:          do i=xs,xs+xm-1
439:             row = (j-gys)*gxm + (i-gxs)

441:             k = 0
442:             if (j .gt. gys) then
443:                v(k) = -1.0d0
444:                col(k) = row-gxm
445:                k = k + 1
446:             endif

448:             if (i .gt. gxs) then
449:                v(k) = -1.0d0
450:                col(k) = row - 1
451:                k = k +1
452:             endif

454:             v(k) = 4.0d0
455:             col(k) = row
456:             k = k + 1

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

464:             if (j+1 .lt. gys + gym) then
465:                v(k) = -1.0d0
466:                col(k) = row + gxm
467:                k = k + 1
468:             endif

470:             call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,info)
471:          enddo
472:       enddo

474:       
475: !     Assemble matrix
476:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
477:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)


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

483:       call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,info)
484:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)


487:       call PetscLogFlops(9*xm*ym + 49*xm,info)

489:       info = 0
490:       return 
491:       end
492:       
493:                
494: