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
68: ! Note: Any user-defined Fortran routines (such as FormGradient)
69: ! MUST be declared as external.
71: external FormInitialGuess,FormFunctionGradient,ComputeHessian
72: external Monitor,ConvergenceTest
74: i1 = 1
76: ! Initialize TAO, PETSc contexts
77: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
78: if (ierr .ne. 0) then
79: print*,'Unable to initialize PETSc'
80: stop
81: endif
83: ! Specify default parameters
84: param = 5.0
85: mx = 10
86: my = 10
87: Nx = PETSC_DECIDE
88: Ny = PETSC_DECIDE
90: ! Check for any command line arguments that might override defaults
91: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
92: & '-mx',mx,flg,ierr)
93: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
94: & '-my',my,flg,ierr)
95: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
96: & '-par',param,flg,ierr)
99: ! Set up distributed array and vectors
100: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, &
101: & DM_BOUNDARY_NONE, &
102: & DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER, &
103: & PETSC_NULL_INTEGER,dm,ierr)
104: call DMSetFromOptions(dm,ierr)
105: call DMSetUp(dm,ierr)
107: ! Create vectors
108: call DMCreateGlobalVector(dm,x,ierr)
109: call DMCreateLocalVector(dm,localX,ierr)
111: ! Create Hessian
112: call DMCreateMatrix(dm,H,ierr)
113: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
115: ! The TAO code begins here
117: ! Create TAO solver
118: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
119: call TaoSetType(tao,TAOCG,ierr)
121: ! Set routines for function and gradient evaluation
123: call TaoSetObjectiveAndGradientRoutine(tao, &
124: & FormFunctionGradient,0,ierr)
125: call TaoSetHessianRoutine(tao,H,H,ComputeHessian, &
126: & 0,ierr)
128: ! Set initial guess
129: call FormInitialGuess(x,ierr)
130: call TaoSetInitialVector(tao,x,ierr)
132: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
133: & '-testmonitor',flg,ierr)
134: if (flg) then
135: call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION, &
136: & ierr)
137: endif
139: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
140: & '-testconvergence',flg, ierr)
141: if (flg) then
142: call TaoSetConvergenceTest(tao,ConvergenceTest,dummy, &
143: & ierr)
144: endif
146: ! Check for any TAO command line options
147: call TaoSetFromOptions(tao,ierr)
150: ! SOLVE THE APPLICATION
151: call TaoSolve(tao,ierr)
153: ! Free TAO data structures
154: call TaoDestroy(tao,ierr)
157: ! Free PETSc data structures
158: call VecDestroy(x,ierr)
159: call VecDestroy(localX,ierr)
160: call MatDestroy(H,ierr)
161: call DMDestroy(dm,ierr)
164: ! Finalize TAO and PETSc
165: call PetscFinalize(ierr)
166: end
169: ! ---------------------------------------------------------------------
170: !
171: ! FormInitialGuess - Computes an initial approximation to the solution.
172: !
173: ! Input Parameters:
174: ! X - vector
175: !
176: ! Output Parameters:
177: ! X - vector
178: ! ierr - error code
179: !
180: subroutine FormInitialGuess(X,ierr)
181: use mymodule
182: implicit none
184: ! Input/output variables:
185: Vec X
186: PetscErrorCode ierr
188: ! Local variables:
189: PetscInt i, j, k, xe, ye
190: PetscReal temp, val, hx, hy
191: PetscInt xs, ys, xm, ym
192: PetscInt gxm, gym, gxs, gys
193: PetscInt i1
195: i1 = 1
196: hx = 1.0/real(mx + 1)
197: hy = 1.0/real(my + 1)
199: ! Get corner information
200: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
201: & PETSC_NULL_INTEGER,ierr)
202: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
203: & gxm,gym,PETSC_NULL_INTEGER,ierr)
207: ! Compute initial guess over locally owned part of mesh
208: xe = xs+xm
209: ye = ys+ym
210: do j=ys,ye-1
211: temp = min(j+1,my-j)*hy
212: do i=xs,xe-1
213: k = (j-gys)*gxm + i-gxs
214: val = min((min(i+1,mx-i))*hx,temp)
215: call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
216: end do
217: end do
218: call VecAssemblyBegin(X,ierr)
219: call VecAssemblyEnd(X,ierr)
220: return
221: end
224: ! ---------------------------------------------------------------------
225: !
226: ! FormFunctionGradient - Evaluates gradient G(X).
227: !
228: ! Input Parameters:
229: ! tao - the Tao context
230: ! X - input vector
231: ! dummy - optional user-defined context (not used here)
232: !
233: ! Output Parameters:
234: ! f - the function value at X
235: ! G - vector containing the newly evaluated gradient
236: ! ierr - error code
237: !
238: ! Notes:
239: ! This routine serves as a wrapper for the lower-level routine
240: ! "ApplicationGradient", where the actual computations are
241: ! done using the standard Fortran style of treating the local
242: ! input vector data as an array over the local mesh.
243: !
244: subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
245: use mymodule
246: implicit none
248: ! Input/output variables:
249: Tao tao
250: Vec X, G
251: PetscReal f
252: PetscErrorCode ierr
253: PetscInt dummy
255: ! Declarations for use with local array:
258: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
259: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
260: ! will return an array of doubles referenced by x_array offset by x_index.
261: ! i.e., to reference the kth element of X, use x_array(k + x_index).
262: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
263: PetscReal lx_v(0:1)
264: PetscOffset lx_i
266: ! Local variables:
267: PetscReal zero, p5, area, cdiv3
268: PetscReal val, flin, fquad,floc
269: PetscReal v, vb, vl, vr, vt, dvdx
270: PetscReal dvdy, hx, hy
271: PetscInt xe, ye, xsm, ysm
272: PetscInt xep, yep, i, j, k, ind
273: PetscInt xs, ys, xm, ym
274: PetscInt gxs, gys, gxm, gym
275: PetscInt i1
277: i1 = 1
278: 0
279: cdiv3 = param/3.0
280: zero = 0.0
281: p5 = 0.5
282: hx = 1.0/real(mx + 1)
283: hy = 1.0/real(my + 1)
284: fquad = zero
285: flin = zero
287: ! Initialize gradient to zero
288: call VecSet(G,zero,ierr)
290: ! Scatter ghost points to local vector
291: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
292: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
295: ! Get corner information
296: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
297: & PETSC_NULL_INTEGER,ierr)
298: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
299: & gxm,gym,PETSC_NULL_INTEGER,ierr)
301: ! Get pointer to vector data.
302: call VecGetArray(localX,lx_v,lx_i,ierr)
305: ! Set local loop dimensions
306: xe = xs+xm
307: ye = ys+ym
308: if (xs .eq. 0) then
309: xsm = xs-1
310: else
311: xsm = xs
312: endif
313: if (ys .eq. 0) then
314: ysm = ys-1
315: else
316: ysm = ys
317: endif
318: if (xe .eq. mx) then
319: xep = xe+1
320: else
321: xep = xe
322: endif
323: if (ye .eq. my) then
324: yep = ye+1
325: else
326: yep = ye
327: endif
329: ! Compute local gradient contributions over the lower triangular elements
331: do j = ysm, ye-1
332: do i = xsm, xe-1
333: k = (j-gys)*gxm + i-gxs
334: v = zero
335: vr = zero
336: vt = zero
337: if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
338: if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
339: if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
340: dvdx = (vr-v)/hx
341: dvdy = (vt-v)/hy
342: if (i .ne. -1 .and. j .ne. -1) then
343: ind = k
344: val = - dvdx/hx - dvdy/hy - cdiv3
345: call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
346: endif
347: if (i .ne. mx-1 .and. j .ne. -1) then
348: ind = k+1
349: val = dvdx/hx - cdiv3
350: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
351: endif
352: if (i .ne. -1 .and. j .ne. my-1) then
353: ind = k+gxm
354: val = dvdy/hy - cdiv3
355: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
356: endif
357: fquad = fquad + dvdx*dvdx + dvdy*dvdy
358: flin = flin - cdiv3 * (v+vr+vt)
359: end do
360: end do
362: ! Compute local gradient contributions over the upper triangular elements
364: do j = ys, yep-1
365: do i = xs, xep-1
366: k = (j-gys)*gxm + i-gxs
367: vb = zero
368: vl = zero
369: v = zero
370: if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
371: if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
372: if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
373: dvdx = (v-vl)/hx
374: dvdy = (v-vb)/hy
375: if (i .ne. mx .and. j .ne. 0) then
376: ind = k-gxm
377: val = - dvdy/hy - cdiv3
378: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
379: endif
380: if (i .ne. 0 .and. j .ne. my) then
381: ind = k-1
382: val = - dvdx/hx - cdiv3
383: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
384: endif
385: if (i .ne. mx .and. j .ne. my) then
386: ind = k
387: val = dvdx/hx + dvdy/hy - cdiv3
388: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
389: endif
390: fquad = fquad + dvdx*dvdx + dvdy*dvdy
391: flin = flin - cdiv3*(vb + vl + v)
392: end do
393: end do
395: ! Restore vector
396: call VecRestoreArray(localX,lx_v,lx_i,ierr)
398: ! Assemble gradient vector
399: call VecAssemblyBegin(G,ierr)
400: call VecAssemblyEnd(G,ierr)
402: ! Scale the gradient
403: area = p5*hx*hy
404: floc = area *(p5*fquad+flin)
405: call VecScale(G,area,ierr)
408: ! Sum function contributions from all processes
409: call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM, &
410: & PETSC_COMM_WORLD,ierr)
411: call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+ &
412: & 16.0d0*(xep-xs)*(yep-ys),ierr)
413: return
414: end
416: subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
417: use mymodule
418: implicit none
420: Tao tao
421: Vec X
422: Mat H,Hpre
423: PetscErrorCode ierr
424: PetscInt dummy
427: PetscInt i,j,k
428: PetscInt col(0:4),row
429: PetscInt xs,xm,gxs,gxm
430: PetscInt ys,ym,gys,gym
431: PetscReal v(0:4)
432: PetscInt i1
434: i1 = 1
436: ! Get local grid boundaries
437: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
438: & PETSC_NULL_INTEGER,ierr)
439: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
440: & PETSC_NULL_INTEGER,ierr)
442: do j=ys,ys+ym-1
443: do i=xs,xs+xm-1
444: row = (j-gys)*gxm + (i-gxs)
446: k = 0
447: if (j .gt. gys) then
448: v(k) = -1.0
449: col(k) = row-gxm
450: k = k + 1
451: endif
453: if (i .gt. gxs) then
454: v(k) = -1.0
455: col(k) = row - 1
456: k = k +1
457: endif
459: v(k) = 4.0
460: col(k) = row
461: k = k + 1
463: if (i+1 .lt. gxs + gxm) then
464: v(k) = -1.0
465: col(k) = row + 1
466: k = k + 1
467: endif
469: if (j+1 .lt. gys + gym) then
470: v(k) = -1.0
471: col(k) = row + gxm
472: k = k + 1
473: endif
475: call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
476: enddo
477: enddo
480: ! Assemble matrix
481: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
482: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
485: ! Tell the matrix we will never add a new nonzero location to the
486: ! matrix. If we do it will generate an error.
488: call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
489: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
492: call PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)
494: 0
495: return
496: end
500: subroutine Monitor(tao, dummy, ierr)
501: use mymodule
502: implicit none
504: Tao tao
505: PetscInt dummy
506: PetscErrorCode ierr
508: PetscInt its
509: PetscReal f,gnorm,cnorm,xdiff
510: TaoConvergedReason reason
512: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
513: & reason,ierr)
514: if (mod(its,5) .ne. 0) then
515: call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n', &
516: & ierr)
517: endif
519: 0
521: return
522: end
524: subroutine ConvergenceTest(tao, dummy, ierr)
525: use mymodule
526: implicit none
528: Tao tao
529: PetscInt dummy
530: PetscErrorCode ierr
532: PetscInt its
533: PetscReal f,gnorm,cnorm,xdiff
534: TaoConvergedReason reason
536: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
537: & reason,ierr)
538: if (its .eq. 7) then
539: call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
540: endif
542: 0
544: return
545: end
547: !/*TEST
548: !
549: ! build:
550: ! requires: !complex
551: !
552: ! test:
553: ! args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
554: !
555: ! test:
556: ! suffix: 2
557: ! nsize: 2
558: ! args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
559: !
560: ! test:
561: ! suffix: 3
562: ! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
563: !TEST*/