Actual source code: eptorsion2f.F
petsc-3.7.5 2017-01-01
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: 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: PetscBool flg
53: PetscInt i1
54: PetscInt dummy
57: ! Note: Any user-defined Fortran routines (such as FormGradient)
58: ! MUST be declared as external.
60: external FormInitialGuess,FormFunctionGradient,ComputeHessian
61: external Monitor,ConvergenceTest
63: i1 = 1
65: ! Initialize TAO, PETSc contexts
66: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
68: ! Specify default parameters
69: param = 5.0
70: mx = 10
71: my = 10
72: Nx = PETSC_DECIDE
73: Ny = PETSC_DECIDE
75: ! Check for any command line arguments that might override defaults
76: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
77: & '-mx',mx,flg,ierr)
78: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
79: & '-my',my,flg,ierr)
80: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
81: & '-par',param,flg,ierr)
84: ! Set up distributed array and vectors
85: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, &
86: & DM_BOUNDARY_NONE, &
87: & DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER, &
88: & PETSC_NULL_INTEGER,dm,ierr)
90: ! Create vectors
91: call DMCreateGlobalVector(dm,x,ierr)
92: call DMCreateLocalVector(dm,localX,ierr)
94: ! Create Hessian
95: call DMCreateMatrix(dm,H,ierr)
96: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
98: ! The TAO code begins here
100: ! Create TAO solver
101: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
102: call TaoSetType(tao,TAOCG,ierr)
104: ! Set routines for function and gradient evaluation
106: call TaoSetObjectiveAndGradientRoutine(tao, &
107: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
108: call TaoSetHessianRoutine(tao,H,H,ComputeHessian, &
109: & PETSC_NULL_OBJECT,ierr)
111: ! Set initial guess
112: call FormInitialGuess(x,ierr)
113: call TaoSetInitialVector(tao,x,ierr)
115: call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
116: & '-testmonitor',flg,ierr)
117: if (flg) then
118: call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION, &
119: & ierr)
120: endif
122: call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
123: & '-testconvergence',flg, ierr)
124: if (flg) then
125: call TaoSetConvergenceTest(tao,ConvergenceTest,dummy, &
126: & ierr)
127: endif
129: ! Check for any TAO command line options
130: call TaoSetFromOptions(tao,ierr)
133: ! SOLVE THE APPLICATION
134: call TaoSolve(tao,ierr)
136: ! Free TAO data structures
137: call TaoDestroy(tao,ierr)
140: ! Free PETSc data structures
141: call VecDestroy(x,ierr)
142: call VecDestroy(localX,ierr)
143: call MatDestroy(H,ierr)
144: call DMDestroy(dm,ierr)
147: ! Finalize TAO and PETSc
148: call PetscFinalize(ierr)
149: end
152: ! ---------------------------------------------------------------------
153: !
154: ! FormInitialGuess - Computes an initial approximation to the solution.
155: !
156: ! Input Parameters:
157: ! X - vector
158: !
159: ! Output Parameters:
160: ! X - vector
161: ! ierr - error code
162: !
163: subroutine FormInitialGuess(X,ierr)
164: implicit none
166: ! mx, my are defined in eptorsion2f.h
167: #include "eptorsion2f.h"
169: ! Input/output variables:
170: Vec X
171: PetscErrorCode ierr
173: ! Local variables:
174: PetscInt i, j, k, xe, ye
175: PetscReal temp, val, hx, hy
176: PetscInt xs, ys, xm, ym
177: PetscInt gxm, gym, gxs, gys
178: PetscInt i1
180: i1 = 1
181: hx = 1.0/(mx + 1)
182: hy = 1.0/(my + 1)
184: ! Get corner information
185: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
186: & PETSC_NULL_INTEGER,ierr)
187: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
188: & gxm,gym,PETSC_NULL_INTEGER,ierr)
192: ! Compute initial guess over locally owned part of mesh
193: xe = xs+xm
194: ye = ys+ym
195: do j=ys,ye-1
196: temp = min(j+1,my-j)*hy
197: do i=xs,xe-1
198: k = (j-gys)*gxm + i-gxs
199: val = min((min(i+1,mx-i))*hx,temp)
200: call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
201: end do
202: end do
203: call VecAssemblyBegin(X,ierr)
204: call VecAssemblyEnd(X,ierr)
205: return
206: end
209: ! ---------------------------------------------------------------------
210: !
211: ! FormFunctionGradient - Evaluates gradient G(X).
212: !
213: ! Input Parameters:
214: ! tao - the Tao context
215: ! X - input vector
216: ! dummy - optional user-defined context (not used here)
217: !
218: ! Output Parameters:
219: ! f - the function value at X
220: ! G - vector containing the newly evaluated gradient
221: ! ierr - error code
222: !
223: ! Notes:
224: ! This routine serves as a wrapper for the lower-level routine
225: ! "ApplicationGradient", where the actual computations are
226: ! done using the standard Fortran style of treating the local
227: ! input vector data as an array over the local mesh.
228: !
229: subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
230: implicit none
232: ! dm, mx, my, param, localX declared in eptorsion2f.h
233: #include "eptorsion2f.h"
235: ! Input/output variables:
236: Tao tao
237: Vec X, G
238: PetscReal f
239: PetscErrorCode ierr
240: PetscInt dummy
242: ! Declarations for use with local array:
245: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
246: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
247: ! will return an array of doubles referenced by x_array offset by x_index.
248: ! i.e., to reference the kth element of X, use x_array(k + x_index).
249: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
250: PetscReal lx_v(0:1)
251: PetscOffset lx_i
253: ! Local variables:
254: PetscReal zero, p5, area, cdiv3
255: PetscReal val, flin, fquad,floc
256: PetscReal v, vb, vl, vr, vt, dvdx
257: PetscReal dvdy, hx, hy
258: PetscInt xe, ye, xsm, ysm
259: PetscInt xep, yep, i, j, k, ind
260: PetscInt xs, ys, xm, ym
261: PetscInt gxs, gys, gxm, gym
262: PetscInt i1
264: i1 = 1
265: 0
266: cdiv3 = param/3.0
267: zero = 0.0
268: p5 = 0.5
269: hx = 1.0/(mx + 1)
270: hy = 1.0/(my + 1)
271: fquad = zero
272: flin = zero
274: ! Initialize gradient to zero
275: call VecSet(G,zero,ierr)
277: ! Scatter ghost points to local vector
278: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
279: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
282: ! Get corner information
283: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
284: & PETSC_NULL_INTEGER,ierr)
285: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
286: & gxm,gym,PETSC_NULL_INTEGER,ierr)
288: ! Get pointer to vector data.
289: call VecGetArray(localX,lx_v,lx_i,ierr)
292: ! Set local loop dimensions
293: xe = xs+xm
294: ye = ys+ym
295: if (xs .eq. 0) then
296: xsm = xs-1
297: else
298: xsm = xs
299: endif
300: if (ys .eq. 0) then
301: ysm = ys-1
302: else
303: ysm = ys
304: endif
305: if (xe .eq. mx) then
306: xep = xe+1
307: else
308: xep = xe
309: endif
310: if (ye .eq. my) then
311: yep = ye+1
312: else
313: yep = ye
314: endif
316: ! Compute local gradient contributions over the lower triangular elements
318: do j = ysm, ye-1
319: do i = xsm, xe-1
320: k = (j-gys)*gxm + i-gxs
321: v = zero
322: vr = zero
323: vt = zero
324: if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
325: if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
326: if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
327: dvdx = (vr-v)/hx
328: dvdy = (vt-v)/hy
329: if (i .ne. -1 .and. j .ne. -1) then
330: ind = k
331: val = - dvdx/hx - dvdy/hy - cdiv3
332: call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
333: endif
334: if (i .ne. mx-1 .and. j .ne. -1) then
335: ind = k+1
336: val = dvdx/hx - cdiv3
337: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
338: endif
339: if (i .ne. -1 .and. j .ne. my-1) then
340: ind = k+gxm
341: val = dvdy/hy - cdiv3
342: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
343: endif
344: fquad = fquad + dvdx*dvdx + dvdy*dvdy
345: flin = flin - cdiv3 * (v+vr+vt)
346: end do
347: end do
349: ! Compute local gradient contributions over the upper triangular elements
351: do j = ys, yep-1
352: do i = xs, xep-1
353: k = (j-gys)*gxm + i-gxs
354: vb = zero
355: vl = zero
356: v = zero
357: if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
358: if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
359: if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
360: dvdx = (v-vl)/hx
361: dvdy = (v-vb)/hy
362: if (i .ne. mx .and. j .ne. 0) then
363: ind = k-gxm
364: val = - dvdy/hy - cdiv3
365: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
366: endif
367: if (i .ne. 0 .and. j .ne. my) then
368: ind = k-1
369: val = - dvdx/hx - cdiv3
370: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
371: endif
372: if (i .ne. mx .and. j .ne. my) then
373: ind = k
374: val = dvdx/hx + dvdy/hy - cdiv3
375: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
376: endif
377: fquad = fquad + dvdx*dvdx + dvdy*dvdy
378: flin = flin - cdiv3*(vb + vl + v)
379: end do
380: end do
382: ! Restore vector
383: call VecRestoreArray(localX,lx_v,lx_i,ierr)
385: ! Assemble gradient vector
386: call VecAssemblyBegin(G,ierr)
387: call VecAssemblyEnd(G,ierr)
389: ! Scale the gradient
390: area = p5*hx*hy
391: floc = area *(p5*fquad+flin)
392: 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
406: subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
407: implicit none
408: #include "eptorsion2f.h"
409: Tao tao
410: Vec X
411: Mat H,Hpre
412: PetscErrorCode ierr
413: PetscInt dummy
416: PetscInt i,j,k
417: PetscInt col(0:4),row
418: PetscInt xs,xm,gxs,gxm
419: PetscInt ys,ym,gys,gym
420: PetscReal v(0:4)
421: PetscInt i1
423: i1 = 1
425: ! Get local grid boundaries
426: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
427: & PETSC_NULL_INTEGER,ierr)
428: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
429: & PETSC_NULL_INTEGER,ierr)
431: do j=ys,ys+ym-1
432: do i=xs,xs+xm-1
433: row = (j-gys)*gxm + (i-gxs)
435: k = 0
436: if (j .gt. gys) then
437: v(k) = -1.0
438: col(k) = row-gxm
439: k = k + 1
440: endif
442: if (i .gt. gxs) then
443: v(k) = -1.0
444: col(k) = row - 1
445: k = k +1
446: endif
448: v(k) = 4.0
449: col(k) = row
450: k = k + 1
452: if (i+1 .lt. gxs + gxm) then
453: v(k) = -1.0
454: col(k) = row + 1
455: k = k + 1
456: endif
458: if (j+1 .lt. gys + gym) then
459: v(k) = -1.0
460: col(k) = row + gxm
461: k = k + 1
462: endif
464: call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
465: enddo
466: enddo
469: ! Assemble matrix
470: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
471: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
474: ! Tell the matrix we will never add a new nonzero location to the
475: ! matrix. If we do it will generate an error.
477: call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
478: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
481: call PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)
483: 0
484: return
485: end
489: subroutine Monitor(tao, dummy, ierr)
490: implicit none
491: #include "eptorsion2f.h"
492: Tao tao
493: PetscInt dummy
494: PetscErrorCode ierr
496: PetscInt its
497: PetscReal f,gnorm,cnorm,xdiff
498: TaoConvergedReason reason
500: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
501: & reason,ierr)
502: if (mod(its,5) .ne. 0) then
503: call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n', &
504: & ierr)
505: endif
507: 0
509: return
510: end
512: subroutine ConvergenceTest(tao, dummy, ierr)
513: implicit none
514: #include "eptorsion2f.h"
515: Tao tao
516: PetscInt dummy
517: PetscErrorCode ierr
519: PetscInt its
520: PetscReal f,gnorm,cnorm,xdiff
521: TaoConvergedReason reason
523: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
524: & reason,ierr)
525: if (its .eq. 7) then
526: call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
527: endif
529: 0
531: return
532: end