Actual source code: eptorsion2f.F
petsc-3.6.1 2015-08-06
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: TaoCreate(); TaoSetType();
14: ! Routines: TaoSetInitialVector();
15: ! Routines: TaoSetObjectiveAndGradientRoutine();
16: ! Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
17: ! Routines: TaoSetMonitor(); TaoSetConvergenceTest()
18: ! Routines: TaoSolve(); TaoGetSolutionStatus()
19: ! Routines: TaoGetConvergedReason(); 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: TaoConvergedReason reason
53: PetscBool flg
54: PetscInt iter ! iteration information
55: PetscReal ff,gnorm,cnorm,xdiff
56: PetscInt i1
57: PetscInt dummy
60: ! Note: Any user-defined Fortran routines (such as FormGradient)
61: ! MUST be declared as external.
63: external FormInitialGuess,FormFunctionGradient,ComputeHessian
64: external Monitor,ConvergenceTest
66: i1 = 1
68: ! Initialize TAO, PETSc contexts
69: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
71: ! Specify default parameters
72: param = 5.0
73: mx = 10
74: my = 10
75: Nx = PETSC_DECIDE
76: Ny = PETSC_DECIDE
78: ! Check for any command line arguments that might override defaults
79: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mx",mx,flg,ierr)
80: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,ierr)
81: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-par", &
82: & param,flg,ierr)
85: ! Set up distributed array and vectors
86: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, &
87: & DM_BOUNDARY_NONE, &
88: & DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER, &
89: & PETSC_NULL_INTEGER,dm,ierr)
91: ! Create vectors
92: call DMCreateGlobalVector(dm,x,ierr)
93: call DMCreateLocalVector(dm,localX,ierr)
95: ! Create Hessian
96: call DMCreateMatrix(dm,H,ierr)
97: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
99: ! The TAO code begins here
101: ! Create TAO solver
102: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
103: call TaoSetType(tao,TAOCG,ierr)
105: ! Set routines for function and gradient evaluation
107: call TaoSetObjectiveAndGradientRoutine(tao, &
108: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
109: call TaoSetHessianRoutine(tao,H,H,ComputeHessian, &
110: & PETSC_NULL_OBJECT,ierr)
112: ! Set initial guess
113: call FormInitialGuess(x,ierr)
114: call TaoSetInitialVector(tao,x,ierr)
116: call PetscOptionsHasName(PETSC_NULL_CHARACTER,"-testmonitor",flg, &
117: & ierr)
118: if (flg) then
119: call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION, &
120: & ierr)
121: endif
123: call PetscOptionsHasName(PETSC_NULL_CHARACTER,"-testconvergence", &
124: & flg, ierr)
125: if (flg) then
126: call TaoSetConvergenceTest(tao,ConvergenceTest,dummy, &
127: & ierr)
128: endif
130: ! Check for any TAO command line options
131: call TaoSetFromOptions(tao,ierr)
134: ! SOLVE THE APPLICATION
135: call TaoSolve(tao,ierr)
137: ! Get information on termination
138: call TaoGetConvergedReason(tao,reason,ierr)
139: if (reason .lt. 0) then
140: call PetscPrintF(PETSC_COMM_WORLD, &
141: & 'TAO did not terminate successfully\n',ierr)
142: endif
145: ! Free TAO data structures
146: call TaoDestroy(tao,ierr)
149: ! Free PETSc data structures
150: call VecDestroy(x,ierr)
151: call VecDestroy(localX,ierr)
152: call MatDestroy(H,ierr)
153: call DMDestroy(dm,ierr)
156: ! Finalize TAO and PETSc
157: call PetscFinalize(ierr)
158: end
161: ! ---------------------------------------------------------------------
162: !
163: ! FormInitialGuess - Computes an initial approximation to the solution.
164: !
165: ! Input Parameters:
166: ! X - vector
167: !
168: ! Output Parameters:
169: ! X - vector
170: ! ierr - error code
171: !
172: subroutine FormInitialGuess(X,ierr)
173: implicit none
175: ! mx, my are defined in eptorsion2f.h
176: #include "eptorsion2f.h"
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/(mx + 1)
191: hy = 1.0/(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)
201: ! Compute initial guess over locally owned part of mesh
202: xe = xs+xm
203: ye = ys+ym
204: do j=ys,ye-1
205: temp = min(j+1,my-j)*hy
206: do i=xs,xe-1
207: k = (j-gys)*gxm + i-gxs
208: val = min((min(i+1,mx-i))*hx,temp)
209: call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
210: end do
211: end do
212: call VecAssemblyBegin(X,ierr)
213: call VecAssemblyEnd(X,ierr)
214: return
215: end
218: ! ---------------------------------------------------------------------
219: !
220: ! FormFunctionGradient - Evaluates gradient G(X).
221: !
222: ! Input Parameters:
223: ! tao - the Tao context
224: ! X - input vector
225: ! dummy - optional user-defined context (not used here)
226: !
227: ! Output Parameters:
228: ! f - the function value at X
229: ! G - vector containing the newly evaluated gradient
230: ! ierr - error code
231: !
232: ! Notes:
233: ! This routine serves as a wrapper for the lower-level routine
234: ! "ApplicationGradient", where the actual computations are
235: ! done using the standard Fortran style of treating the local
236: ! input vector data as an array over the local mesh.
237: !
238: subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
239: implicit none
241: ! dm, mx, my, param, localX declared in eptorsion2f.h
242: #include "eptorsion2f.h"
244: ! Input/output variables:
245: Tao tao
246: Vec X, G
247: PetscReal f
248: PetscErrorCode ierr
249: PetscInt dummy
251: ! Declarations for use with local array:
254: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
255: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
256: ! will return an array of doubles referenced by x_array offset by x_index.
257: ! i.e., to reference the kth element of X, use x_array(k + x_index).
258: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
259: PetscReal lx_v(0:1)
260: PetscOffset lx_i
262: ! Local variables:
263: PetscReal zero, p5, area, cdiv3
264: PetscReal val, flin, fquad,floc
265: PetscReal v, vb, vl, vr, vt, dvdx
266: PetscReal dvdy, hx, hy
267: PetscInt xe, ye, xsm, ysm
268: PetscInt xep, yep, i, j, k, ind
269: PetscInt xs, ys, xm, ym
270: PetscInt gxs, gys, gxm, gym
271: PetscInt i1
273: i1 = 1
274: 0
275: cdiv3 = param/3.0
276: zero = 0.0
277: p5 = 0.5
278: hx = 1.0/(mx + 1)
279: hy = 1.0/(my + 1)
280: fquad = zero
281: flin = zero
283: ! Initialize gradient to zero
284: call VecSet(G,zero,ierr)
286: ! Scatter ghost points to local vector
287: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
288: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
291: ! Get corner information
292: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
293: & PETSC_NULL_INTEGER,ierr)
294: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
295: & gxm,gym,PETSC_NULL_INTEGER,ierr)
297: ! Get pointer to vector data.
298: call VecGetArray(localX,lx_v,lx_i,ierr)
301: ! Set local loop dimensions
302: xe = xs+xm
303: ye = ys+ym
304: if (xs .eq. 0) then
305: xsm = xs-1
306: else
307: xsm = xs
308: endif
309: if (ys .eq. 0) then
310: ysm = ys-1
311: else
312: ysm = ys
313: endif
314: if (xe .eq. mx) then
315: xep = xe+1
316: else
317: xep = xe
318: endif
319: if (ye .eq. my) then
320: yep = ye+1
321: else
322: yep = ye
323: endif
325: ! Compute local gradient contributions over the lower triangular elements
327: do j = ysm, ye-1
328: do i = xsm, xe-1
329: k = (j-gys)*gxm + i-gxs
330: v = zero
331: vr = zero
332: vt = zero
333: if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
334: if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
335: if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
336: dvdx = (vr-v)/hx
337: dvdy = (vt-v)/hy
338: if (i .ne. -1 .and. j .ne. -1) then
339: ind = k
340: val = - dvdx/hx - dvdy/hy - cdiv3
341: call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
342: endif
343: if (i .ne. mx-1 .and. j .ne. -1) then
344: ind = k+1
345: val = dvdx/hx - cdiv3
346: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
347: endif
348: if (i .ne. -1 .and. j .ne. my-1) then
349: ind = k+gxm
350: val = dvdy/hy - cdiv3
351: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
352: endif
353: fquad = fquad + dvdx*dvdx + dvdy*dvdy
354: flin = flin - cdiv3 * (v+vr+vt)
355: end do
356: end do
358: ! Compute local gradient contributions over the upper triangular elements
360: do j = ys, yep-1
361: do i = xs, xep-1
362: k = (j-gys)*gxm + i-gxs
363: vb = zero
364: vl = zero
365: v = zero
366: if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
367: if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
368: if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
369: dvdx = (v-vl)/hx
370: dvdy = (v-vb)/hy
371: if (i .ne. mx .and. j .ne. 0) then
372: ind = k-gxm
373: val = - dvdy/hy - cdiv3
374: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
375: endif
376: if (i .ne. 0 .and. j .ne. my) then
377: ind = k-1
378: val = - dvdx/hx - cdiv3
379: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
380: endif
381: if (i .ne. mx .and. j .ne. my) then
382: ind = k
383: val = dvdx/hx + dvdy/hy - cdiv3
384: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
385: endif
386: fquad = fquad + dvdx*dvdx + dvdy*dvdy
387: flin = flin - cdiv3*(vb + vl + v)
388: end do
389: end do
391: ! Restore vector
392: call VecRestoreArray(localX,lx_v,lx_i,ierr)
394: ! Assemble gradient vector
395: call VecAssemblyBegin(G,ierr)
396: call VecAssemblyEnd(G,ierr)
398: ! Scale the gradient
399: area = p5*hx*hy
400: floc = area *(p5*fquad+flin)
401: call VecScale(G,area,ierr)
404: ! Sum function contributions from all processes
405: call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM, &
406: & PETSC_COMM_WORLD,ierr)
407: call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+ &
408: & 16.0d0*(xep-xs)*(yep-ys),ierr)
409: return
410: end
415: subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
416: implicit none
417: #include "eptorsion2f.h"
418: Tao tao
419: Vec X
420: Mat H,Hpre
421: PetscErrorCode ierr
422: PetscInt dummy
425: PetscInt i,j,k
426: PetscInt col(0:4),row
427: PetscInt xs,xm,gxs,gxm
428: PetscInt ys,ym,gys,gym
429: PetscReal v(0:4)
430: PetscInt i1
432: i1 = 1
434: ! Get local grid boundaries
435: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
436: & PETSC_NULL_INTEGER,ierr)
437: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
438: & PETSC_NULL_INTEGER,ierr)
440: do j=ys,ys+ym-1
441: do i=xs,xs+xm-1
442: row = (j-gys)*gxm + (i-gxs)
444: k = 0
445: if (j .gt. gys) then
446: v(k) = -1.0
447: col(k) = row-gxm
448: k = k + 1
449: endif
451: if (i .gt. gxs) then
452: v(k) = -1.0
453: col(k) = row - 1
454: k = k +1
455: endif
457: v(k) = 4.0
458: col(k) = row
459: k = k + 1
461: if (i+1 .lt. gxs + gxm) then
462: v(k) = -1.0
463: col(k) = row + 1
464: k = k + 1
465: endif
467: if (j+1 .lt. gys + gym) then
468: v(k) = -1.0
469: col(k) = row + gxm
470: k = k + 1
471: endif
473: call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
474: enddo
475: enddo
478: ! Assemble matrix
479: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
480: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
483: ! Tell the matrix we will never add a new nonzero location to the
484: ! matrix. If we do it will generate an error.
486: call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
487: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
490: call PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)
492: 0
493: return
494: end
498: subroutine Monitor(tao, dummy, ierr)
499: implicit none
500: #include "eptorsion2f.h"
501: Tao tao
502: PetscInt dummy
503: PetscErrorCode ierr
505: PetscInt its
506: PetscReal f,gnorm,cnorm,xdiff
507: TaoConvergedReason reason
509: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
510: & reason,ierr)
511: if (mod(its,5) .ne. 0) then
512: call PetscPrintf(PETSC_COMM_WORLD,"iteration multiple of 5\n", &
513: & ierr)
514: endif
516: 0
518: return
519: end
521: subroutine ConvergenceTest(tao, dummy, ierr)
522: implicit none
523: #include "eptorsion2f.h"
524: Tao tao
525: PetscInt dummy
526: PetscErrorCode ierr
528: PetscInt its
529: PetscReal f,gnorm,cnorm,xdiff
530: TaoConvergedReason reason
532: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
533: & reason,ierr)
534: if (its .eq. 7) then
535: call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
536: endif
538: 0
540: return
541: end