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: !
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 mymodule
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 mymodule
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: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
65: if (ierr .ne. 0) then
66: print*,'Unable to initialize PETSc'
67: stop
68: endif
70: ! Specify default parameters
71: param = 5.0
72: mx = 10
73: my = 10
74: Nx = PETSC_DECIDE
75: Ny = PETSC_DECIDE
77: ! Check for any command line arguments that might override defaults
78: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
79: & '-mx',mx,flg,ierr)
80: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
81: & '-my',my,flg,ierr)
82: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
83: & '-par',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)
90: call DMSetFromOptions(dm,ierr)
91: call DMSetUp(dm,ierr)
93: ! Create vectors
94: call DMCreateGlobalVector(dm,x,ierr)
95: call DMCreateLocalVector(dm,localX,ierr)
97: ! Create Hessian
98: call DMCreateMatrix(dm,H,ierr)
99: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
101: ! The TAO code begins here
103: ! Create TAO solver
104: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
105: call TaoSetType(tao,TAOCG,ierr)
107: ! Set routines for function and gradient evaluation
109: call TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC, &
110: & FormFunctionGradient,0,ierr)
111: call TaoSetHessian(tao,H,H,ComputeHessian, &
112: & 0,ierr)
114: ! Set initial guess
115: call FormInitialGuess(x,ierr)
116: call TaoSetSolution(tao,x,ierr)
118: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
119: & '-testmonitor',flg,ierr)
120: if (flg) then
121: call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION, &
122: & ierr)
123: endif
125: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
126: & '-testconvergence',flg, ierr)
127: if (flg) then
128: call TaoSetConvergenceTest(tao,ConvergenceTest,dummy, &
129: & ierr)
130: endif
132: ! Check for any TAO command line options
133: call TaoSetFromOptions(tao,ierr)
135: ! SOLVE THE APPLICATION
136: call TaoSolve(tao,ierr)
138: ! Free TAO data structures
139: call TaoDestroy(tao,ierr)
141: ! Free PETSc data structures
142: call VecDestroy(x,ierr)
143: call VecDestroy(localX,ierr)
144: call MatDestroy(H,ierr)
145: call DMDestroy(dm,ierr)
147: ! Finalize TAO and PETSc
148: call PetscFinalize(ierr)
149: end
151: ! ---------------------------------------------------------------------
152: !
153: ! FormInitialGuess - Computes an initial approximation to the solution.
154: !
155: ! Input Parameters:
156: ! X - vector
157: !
158: ! Output Parameters:
159: ! X - vector
160: ! ierr - error code
161: !
162: subroutine FormInitialGuess(X,ierr)
163: use mymodule
164: implicit none
166: ! Input/output variables:
167: Vec X
168: PetscErrorCode ierr
170: ! Local variables:
171: PetscInt i, j, k, xe, ye
172: PetscReal temp, val, hx, hy
173: PetscInt xs, ys, xm, ym
174: PetscInt gxm, gym, gxs, gys
175: PetscInt i1
177: i1 = 1
178: hx = 1.0/real(mx + 1)
179: hy = 1.0/real(my + 1)
181: ! Get corner information
182: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
183: & PETSC_NULL_INTEGER,ierr)
184: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
185: & gxm,gym,PETSC_NULL_INTEGER,ierr)
187: ! Compute initial guess over locally owned part of mesh
188: xe = xs+xm
189: ye = ys+ym
190: do j=ys,ye-1
191: temp = min(j+1,my-j)*hy
192: do i=xs,xe-1
193: k = (j-gys)*gxm + i-gxs
194: val = min((min(i+1,mx-i))*hx,temp)
195: call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
196: end do
197: end do
198: call VecAssemblyBegin(X,ierr)
199: call VecAssemblyEnd(X,ierr)
200: return
201: end
203: ! ---------------------------------------------------------------------
204: !
205: ! FormFunctionGradient - Evaluates gradient G(X).
206: !
207: ! Input Parameters:
208: ! tao - the Tao context
209: ! X - input vector
210: ! dummy - optional user-defined context (not used here)
211: !
212: ! Output Parameters:
213: ! f - the function value at X
214: ! G - vector containing the newly evaluated gradient
215: ! ierr - error code
216: !
217: ! Notes:
218: ! This routine serves as a wrapper for the lower-level routine
219: ! "ApplicationGradient", where the actual computations are
220: ! done using the standard Fortran style of treating the local
221: ! input vector data as an array over the local mesh.
222: !
223: subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
224: use mymodule
225: implicit none
227: ! Input/output variables:
228: Tao tao
229: Vec X, G
230: PetscReal f
231: PetscErrorCode ierr
232: PetscInt dummy
234: ! Declarations for use with local array:
236: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
237: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
238: ! will return an array of doubles referenced by x_array offset by x_index.
239: ! i.e., to reference the kth element of X, use x_array(k + x_index).
240: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
241: PetscReal lx_v(0:1)
242: PetscOffset lx_i
244: ! Local variables:
245: PetscReal zero, p5, area, cdiv3
246: PetscReal val, flin, fquad,floc
247: PetscReal v, vb, vl, vr, vt, dvdx
248: PetscReal dvdy, hx, hy
249: PetscInt xe, ye, xsm, ysm
250: PetscInt xep, yep, i, j, k, ind
251: PetscInt xs, ys, xm, ym
252: PetscInt gxs, gys, gxm, gym
253: PetscInt i1
255: i1 = 1
256: 0
257: cdiv3 = param/3.0
258: zero = 0.0
259: p5 = 0.5
260: hx = 1.0/real(mx + 1)
261: hy = 1.0/real(my + 1)
262: fquad = zero
263: flin = zero
265: ! Initialize gradient to zero
266: call VecSet(G,zero,ierr)
268: ! Scatter ghost points to local vector
269: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
270: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
272: ! Get corner information
273: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
274: & PETSC_NULL_INTEGER,ierr)
275: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
276: & gxm,gym,PETSC_NULL_INTEGER,ierr)
278: ! Get pointer to vector data.
279: call VecGetArray(localX,lx_v,lx_i,ierr)
281: ! Set local loop dimensions
282: xe = xs+xm
283: ye = ys+ym
284: if (xs .eq. 0) then
285: xsm = xs-1
286: else
287: xsm = xs
288: endif
289: if (ys .eq. 0) then
290: ysm = ys-1
291: else
292: ysm = ys
293: endif
294: if (xe .eq. mx) then
295: xep = xe+1
296: else
297: xep = xe
298: endif
299: if (ye .eq. my) then
300: yep = ye+1
301: else
302: yep = ye
303: endif
305: ! Compute local gradient contributions over the lower triangular elements
307: do j = ysm, ye-1
308: do i = xsm, xe-1
309: k = (j-gys)*gxm + i-gxs
310: v = zero
311: vr = zero
312: vt = zero
313: if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
314: if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
315: if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
316: dvdx = (vr-v)/hx
317: dvdy = (vt-v)/hy
318: if (i .ne. -1 .and. j .ne. -1) then
319: ind = k
320: val = - dvdx/hx - dvdy/hy - cdiv3
321: call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
322: endif
323: if (i .ne. mx-1 .and. j .ne. -1) then
324: ind = k+1
325: val = dvdx/hx - cdiv3
326: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
327: endif
328: if (i .ne. -1 .and. j .ne. my-1) then
329: ind = k+gxm
330: val = dvdy/hy - cdiv3
331: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
332: endif
333: fquad = fquad + dvdx*dvdx + dvdy*dvdy
334: flin = flin - cdiv3 * (v+vr+vt)
335: end do
336: end do
338: ! Compute local gradient contributions over the upper triangular elements
340: do j = ys, yep-1
341: do i = xs, xep-1
342: k = (j-gys)*gxm + i-gxs
343: vb = zero
344: vl = zero
345: v = zero
346: if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
347: if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
348: if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
349: dvdx = (v-vl)/hx
350: dvdy = (v-vb)/hy
351: if (i .ne. mx .and. j .ne. 0) then
352: ind = k-gxm
353: val = - dvdy/hy - cdiv3
354: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
355: endif
356: if (i .ne. 0 .and. j .ne. my) then
357: ind = k-1
358: val = - dvdx/hx - cdiv3
359: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
360: endif
361: if (i .ne. mx .and. j .ne. my) then
362: ind = k
363: val = dvdx/hx + dvdy/hy - cdiv3
364: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
365: endif
366: fquad = fquad + dvdx*dvdx + dvdy*dvdy
367: flin = flin - cdiv3*(vb + vl + v)
368: end do
369: end do
371: ! Restore vector
372: call VecRestoreArray(localX,lx_v,lx_i,ierr)
374: ! Assemble gradient vector
375: call VecAssemblyBegin(G,ierr)
376: call VecAssemblyEnd(G,ierr)
378: ! Scale the gradient
379: area = p5*hx*hy
380: floc = area *(p5*fquad+flin)
381: call VecScale(G,area,ierr)
383: ! Sum function contributions from all processes
384: call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM, &
385: & PETSC_COMM_WORLD,ierr)
386: call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+ &
387: & 16.0d0*(xep-xs)*(yep-ys),ierr)
388: return
389: end
391: subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
392: use mymodule
393: implicit none
395: Tao tao
396: Vec X
397: Mat H,Hpre
398: PetscErrorCode ierr
399: PetscInt dummy
401: PetscInt i,j,k
402: PetscInt col(0:4),row
403: PetscInt xs,xm,gxs,gxm
404: PetscInt ys,ym,gys,gym
405: PetscReal v(0:4)
406: PetscInt i1
408: i1 = 1
410: ! Get local grid boundaries
411: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
412: & PETSC_NULL_INTEGER,ierr)
413: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
414: & PETSC_NULL_INTEGER,ierr)
416: do j=ys,ys+ym-1
417: do i=xs,xs+xm-1
418: row = (j-gys)*gxm + (i-gxs)
420: k = 0
421: if (j .gt. gys) then
422: v(k) = -1.0
423: col(k) = row-gxm
424: k = k + 1
425: endif
427: if (i .gt. gxs) then
428: v(k) = -1.0
429: col(k) = row - 1
430: k = k +1
431: endif
433: v(k) = 4.0
434: col(k) = row
435: k = k + 1
437: if (i+1 .lt. gxs + gxm) then
438: v(k) = -1.0
439: col(k) = row + 1
440: k = k + 1
441: endif
443: if (j+1 .lt. gys + gym) then
444: v(k) = -1.0
445: col(k) = row + gxm
446: k = k + 1
447: endif
449: call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
450: enddo
451: enddo
453: ! Assemble matrix
454: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
455: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
457: ! Tell the matrix we will never add a new nonzero location to the
458: ! matrix. If we do it will generate an error.
460: call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
461: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
463: call PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)
465: 0
466: return
467: end
469: subroutine Monitor(tao, dummy, ierr)
470: use mymodule
471: implicit none
473: Tao tao
474: PetscInt dummy
475: PetscErrorCode ierr
477: PetscInt its
478: PetscReal f,gnorm,cnorm,xdiff
479: TaoConvergedReason reason
481: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
482: & reason,ierr)
483: if (mod(its,5) .ne. 0) then
484: call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n', &
485: & ierr)
486: endif
488: 0
490: return
491: end
493: subroutine ConvergenceTest(tao, dummy, ierr)
494: use mymodule
495: implicit none
497: Tao tao
498: PetscInt dummy
499: PetscErrorCode ierr
501: PetscInt its
502: PetscReal f,gnorm,cnorm,xdiff
503: TaoConvergedReason reason
505: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
506: & reason,ierr)
507: if (its .eq. 7) then
508: call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
509: endif
511: 0
513: return
514: end
516: !/*TEST
517: !
518: ! build:
519: ! requires: !complex
520: !
521: ! test:
522: ! args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
523: !
524: ! test:
525: ! suffix: 2
526: ! nsize: 2
527: ! args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
528: !
529: ! test:
530: ! suffix: 3
531: ! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
532: !TEST*/