Actual source code: eptorsion2f.F90
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: !
12: ! ----------------------------------------------------------------------
13: !
14: ! Elastic-plastic torsion problem.
15: !
16: ! The elastic plastic torsion problem arises from the deconverged
17: ! of the stress field on an infinitely long cylindrical bar, which is
18: ! equivalent to the solution of the following problem:
19: ! min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
20: ! where C is the torsion angle per unit length.
21: !
22: ! The C version of this code is eptorsion2.c
23: !
24: ! ----------------------------------------------------------------------
26: module eptorsion2fmodule
27: #include "petsc/finclude/petsctao.h"
28: use petscdmda
29: use petsctao
30: implicit none
32: Vec localX
33: DM dm
34: PetscReal param
35: PetscInt mx, my
36: end module
38: use eptorsion2fmodule
39: implicit none
40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: ! Variable declarations
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: !
44: ! See additional variable declarations in the file eptorsion2f.h
45: !
46: PetscErrorCode ierr ! used to check for functions returning nonzeros
47: Vec x ! solution vector
48: Mat H ! hessian matrix
49: PetscInt Nx, Ny ! number of processes in x- and y- directions
50: Tao tao ! Tao solver context
51: PetscBool flg
52: PetscInt i1
53: PetscInt dummy
55: ! Note: Any user-defined Fortran routines (such as FormGradient)
56: ! MUST be declared as external.
58: external FormInitialGuess,FormFunctionGradient,ComputeHessian
59: external Monitor,ConvergenceTest
61: i1 = 1
63: ! Initialize TAO, PETSc contexts
64: PetscCallA(PetscInitialize(ierr))
66: ! Specify default parameters
67: param = 5.0
68: mx = 10
69: my = 10
70: Nx = PETSC_DECIDE
71: Ny = PETSC_DECIDE
73: ! Check for any command line arguments that might override defaults
74: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
75: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
76: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',param,flg,ierr))
78: ! Set up distributed array and vectors
79: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,dm,ierr))
80: PetscCallA(DMSetFromOptions(dm,ierr))
81: PetscCallA(DMSetUp(dm,ierr))
83: ! Create vectors
84: PetscCallA(DMCreateGlobalVector(dm,x,ierr))
85: PetscCallA(DMCreateLocalVector(dm,localX,ierr))
87: ! Create Hessian
88: PetscCallA(DMCreateMatrix(dm,H,ierr))
89: PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))
91: ! The TAO code begins here
93: ! Create TAO solver
94: PetscCallA(TaoCreate(PETSC_COMM_WORLD,tao,ierr))
95: PetscCallA(TaoSetType(tao,TAOCG,ierr))
97: ! Set routines for function and gradient evaluation
99: PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr))
100: PetscCallA(TaoSetHessian(tao,H,H,ComputeHessian,0,ierr))
102: ! Set initial guess
103: PetscCallA(FormInitialGuess(x,ierr))
104: PetscCallA(TaoSetSolution(tao,x,ierr))
106: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testmonitor',flg,ierr))
107: if (flg) then
108: PetscCallA(TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION,ierr))
109: endif
111: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testconvergence',flg, ierr))
112: if (flg) then
113: PetscCallA(TaoSetConvergenceTest(tao,ConvergenceTest,dummy,ierr))
114: endif
116: ! Check for any TAO command line options
117: PetscCallA(TaoSetFromOptions(tao,ierr))
119: ! SOLVE THE APPLICATION
120: PetscCallA(TaoSolve(tao,ierr))
122: ! Free TAO data structures
123: PetscCallA(TaoDestroy(tao,ierr))
125: ! Free PETSc data structures
126: PetscCallA(VecDestroy(x,ierr))
127: PetscCallA(VecDestroy(localX,ierr))
128: PetscCallA(MatDestroy(H,ierr))
129: PetscCallA(DMDestroy(dm,ierr))
131: ! Finalize TAO and PETSc
132: PetscCallA(PetscFinalize(ierr))
133: end
135: ! ---------------------------------------------------------------------
136: !
137: ! FormInitialGuess - Computes an initial approximation to the solution.
138: !
139: ! Input Parameters:
140: ! X - vector
141: !
142: ! Output Parameters:
143: ! X - vector
144: ! ierr - error code
145: !
146: subroutine FormInitialGuess(X,ierr)
147: use eptorsion2fmodule
148: implicit none
150: ! Input/output variables:
151: Vec X
152: PetscErrorCode ierr
154: ! Local variables:
155: PetscInt i, j, k, xe, ye
156: PetscReal temp, val, hx, hy
157: PetscInt xs, ys, xm, ym
158: PetscInt gxm, gym, gxs, gys
159: PetscInt i1
161: i1 = 1
162: hx = 1.0/real(mx + 1)
163: hy = 1.0/real(my + 1)
165: ! Get corner information
166: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
167: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
169: ! Compute initial guess over locally owned part of mesh
170: xe = xs+xm
171: ye = ys+ym
172: do j=ys,ye-1
173: temp = min(j+1,my-j)*hy
174: do i=xs,xe-1
175: k = (j-gys)*gxm + i-gxs
176: val = min((min(i+1,mx-i))*hx,temp)
177: PetscCall(VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr))
178: end do
179: end do
180: PetscCall(VecAssemblyBegin(X,ierr))
181: PetscCall(VecAssemblyEnd(X,ierr))
182: return
183: end
185: ! ---------------------------------------------------------------------
186: !
187: ! FormFunctionGradient - Evaluates gradient G(X).
188: !
189: ! Input Parameters:
190: ! tao - the Tao context
191: ! X - input vector
192: ! dummy - optional user-defined context (not used here)
193: !
194: ! Output Parameters:
195: ! f - the function value at X
196: ! G - vector containing the newly evaluated gradient
197: ! ierr - error code
198: !
199: ! Notes:
200: ! This routine serves as a wrapper for the lower-level routine
201: ! "ApplicationGradient", where the actual computations are
202: ! done using the standard Fortran style of treating the local
203: ! input vector data as an array over the local mesh.
204: !
205: subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
206: use eptorsion2fmodule
207: implicit none
209: ! Input/output variables:
210: Tao tao
211: Vec X, G
212: PetscReal f
213: PetscErrorCode ierr
214: PetscInt dummy
216: ! Declarations for use with local array:
218: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
219: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
220: ! will return an array of doubles referenced by x_array offset by x_index.
221: ! i.e., to reference the kth element of X, use x_array(k + x_index).
222: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
223: PetscReal lx_v(0:1)
224: PetscOffset lx_i
226: ! Local variables:
227: PetscReal zero, p5, area, cdiv3
228: PetscReal val, flin, fquad,floc
229: PetscReal v, vb, vl, vr, vt, dvdx
230: PetscReal dvdy, hx, hy
231: PetscInt xe, ye, xsm, ysm
232: PetscInt xep, yep, i, j, k, ind
233: PetscInt xs, ys, xm, ym
234: PetscInt gxs, gys, gxm, gym
235: PetscInt i1
237: i1 = 1
238: 0
239: cdiv3 = param/3.0
240: zero = 0.0
241: p5 = 0.5
242: hx = 1.0/real(mx + 1)
243: hy = 1.0/real(my + 1)
244: fquad = zero
245: flin = zero
247: ! Initialize gradient to zero
248: PetscCall(VecSet(G,zero,ierr))
250: ! Scatter ghost points to local vector
251: PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
252: PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))
254: ! Get corner information
255: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
256: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
258: ! Get pointer to vector data.
259: PetscCall(VecGetArray(localX,lx_v,lx_i,ierr))
261: ! Set local loop dimensions
262: xe = xs+xm
263: ye = ys+ym
264: if (xs .eq. 0) then
265: xsm = xs-1
266: else
267: xsm = xs
268: endif
269: if (ys .eq. 0) then
270: ysm = ys-1
271: else
272: ysm = ys
273: endif
274: if (xe .eq. mx) then
275: xep = xe+1
276: else
277: xep = xe
278: endif
279: if (ye .eq. my) then
280: yep = ye+1
281: else
282: yep = ye
283: endif
285: ! Compute local gradient contributions over the lower triangular elements
287: do j = ysm, ye-1
288: do i = xsm, xe-1
289: k = (j-gys)*gxm + i-gxs
290: v = zero
291: vr = zero
292: vt = zero
293: if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
294: if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
295: if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
296: dvdx = (vr-v)/hx
297: dvdy = (vt-v)/hy
298: if (i .ne. -1 .and. j .ne. -1) then
299: ind = k
300: val = - dvdx/hx - dvdy/hy - cdiv3
301: PetscCall(VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr))
302: endif
303: if (i .ne. mx-1 .and. j .ne. -1) then
304: ind = k+1
305: val = dvdx/hx - cdiv3
306: PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
307: endif
308: if (i .ne. -1 .and. j .ne. my-1) then
309: ind = k+gxm
310: val = dvdy/hy - cdiv3
311: PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
312: endif
313: fquad = fquad + dvdx*dvdx + dvdy*dvdy
314: flin = flin - cdiv3 * (v+vr+vt)
315: end do
316: end do
318: ! Compute local gradient contributions over the upper triangular elements
320: do j = ys, yep-1
321: do i = xs, xep-1
322: k = (j-gys)*gxm + i-gxs
323: vb = zero
324: vl = zero
325: v = zero
326: if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
327: if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
328: if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
329: dvdx = (v-vl)/hx
330: dvdy = (v-vb)/hy
331: if (i .ne. mx .and. j .ne. 0) then
332: ind = k-gxm
333: val = - dvdy/hy - cdiv3
334: PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
335: endif
336: if (i .ne. 0 .and. j .ne. my) then
337: ind = k-1
338: val = - dvdx/hx - cdiv3
339: PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
340: endif
341: if (i .ne. mx .and. j .ne. my) then
342: ind = k
343: val = dvdx/hx + dvdy/hy - cdiv3
344: PetscCall(VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr))
345: endif
346: fquad = fquad + dvdx*dvdx + dvdy*dvdy
347: flin = flin - cdiv3*(vb + vl + v)
348: end do
349: end do
351: ! Restore vector
352: PetscCall(VecRestoreArray(localX,lx_v,lx_i,ierr))
354: ! Assemble gradient vector
355: PetscCall(VecAssemblyBegin(G,ierr))
356: PetscCall(VecAssemblyEnd(G,ierr))
358: ! Scale the gradient
359: area = p5*hx*hy
360: floc = area *(p5*fquad+flin)
361: PetscCall(VecScale(G,area,ierr))
363: ! Sum function contributions from all processes
364: PetscCallMPI(MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr))
365: PetscCall(PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+16.0d0*(xep-xs)*(yep-ys),ierr))
366: return
367: end
369: subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
370: use eptorsion2fmodule
371: implicit none
373: Tao tao
374: Vec X
375: Mat H,Hpre
376: PetscErrorCode ierr
377: PetscInt dummy
379: PetscInt i,j,k
380: PetscInt col(0:4),row
381: PetscInt xs,xm,gxs,gxm
382: PetscInt ys,ym,gys,gym
383: PetscReal v(0:4)
384: PetscInt i1
386: i1 = 1
388: ! Get local grid boundaries
389: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
390: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
392: do j=ys,ys+ym-1
393: do i=xs,xs+xm-1
394: row = (j-gys)*gxm + (i-gxs)
396: k = 0
397: if (j .gt. gys) then
398: v(k) = -1.0
399: col(k) = row-gxm
400: k = k + 1
401: endif
403: if (i .gt. gxs) then
404: v(k) = -1.0
405: col(k) = row - 1
406: k = k +1
407: endif
409: v(k) = 4.0
410: col(k) = row
411: k = k + 1
413: if (i+1 .lt. gxs + gxm) then
414: v(k) = -1.0
415: col(k) = row + 1
416: k = k + 1
417: endif
419: if (j+1 .lt. gys + gym) then
420: v(k) = -1.0
421: col(k) = row + gxm
422: k = k + 1
423: endif
425: PetscCall(MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr))
426: enddo
427: enddo
429: ! Assemble matrix
430: PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
431: PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))
433: ! Tell the matrix we will never add a new nonzero location to the
434: ! matrix. If we do it will generate an error.
436: PetscCall(MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
437: PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))
439: PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr))
441: 0
442: return
443: end
445: subroutine Monitor(tao, dummy, ierr)
446: use eptorsion2fmodule
447: implicit none
449: Tao tao
450: PetscInt dummy
451: PetscErrorCode ierr
453: PetscInt its
454: PetscReal f,gnorm,cnorm,xdiff
455: TaoConvergedReason reason
457: PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
458: if (mod(its,5) .ne. 0) then
459: PetscCall(PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',ierr))
460: endif
462: 0
464: return
465: end
467: subroutine ConvergenceTest(tao, dummy, ierr)
468: use eptorsion2fmodule
469: implicit none
471: Tao tao
472: PetscInt dummy
473: PetscErrorCode ierr
475: PetscInt its
476: PetscReal f,gnorm,cnorm,xdiff
477: TaoConvergedReason reason
479: PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
480: if (its .eq. 7) then
481: PetscCall(TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr))
482: endif
484: 0
486: return
487: end
489: !/*TEST
490: !
491: ! build:
492: ! requires: !complex
493: !
494: ! test:
495: ! args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
496: !
497: ! test:
498: ! suffix: 2
499: ! nsize: 2
500: ! args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
501: !
502: ! test:
503: ! suffix: 3
504: ! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
505: !TEST*/