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*/