Actual source code: eptorsion2f.F
petsc-3.12.5 2020-03-29
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: #if defined(PETSC_USING_F90) && !defined(PETSC_USE_FORTRANKIND)
53: external PETSC_NULL_FUNCTION
54: #endif
55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: ! Variable declarations
57: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: !
59: ! See additional variable declarations in the file eptorsion2f.h
60: !
61: PetscErrorCode ierr ! used to check for functions returning nonzeros
62: Vec x ! solution vector
63: Mat H ! hessian matrix
64: PetscInt Nx, Ny ! number of processes in x- and y- directions
65: Tao tao ! Tao solver context
66: PetscBool flg
67: PetscInt i1
68: PetscInt dummy
71: ! Note: Any user-defined Fortran routines (such as FormGradient)
72: ! MUST be declared as external.
74: external FormInitialGuess,FormFunctionGradient,ComputeHessian
75: external Monitor,ConvergenceTest
77: i1 = 1
79: ! Initialize TAO, PETSc contexts
80: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
81: if (ierr .ne. 0) then
82: print*,'Unable to initialize PETSc'
83: stop
84: endif
86: ! Specify default parameters
87: param = 5.0
88: mx = 10
89: my = 10
90: Nx = PETSC_DECIDE
91: Ny = PETSC_DECIDE
93: ! Check for any command line arguments that might override defaults
94: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
95: & '-mx',mx,flg,ierr)
96: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
97: & '-my',my,flg,ierr)
98: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
99: & '-par',param,flg,ierr)
102: ! Set up distributed array and vectors
103: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, &
104: & DM_BOUNDARY_NONE, &
105: & DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER, &
106: & PETSC_NULL_INTEGER,dm,ierr)
107: call DMSetFromOptions(dm,ierr)
108: call DMSetUp(dm,ierr)
110: ! Create vectors
111: call DMCreateGlobalVector(dm,x,ierr)
112: call DMCreateLocalVector(dm,localX,ierr)
114: ! Create Hessian
115: call DMCreateMatrix(dm,H,ierr)
116: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
118: ! The TAO code begins here
120: ! Create TAO solver
121: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
122: call TaoSetType(tao,TAOCG,ierr)
124: ! Set routines for function and gradient evaluation
126: call TaoSetObjectiveAndGradientRoutine(tao, &
127: & FormFunctionGradient,0,ierr)
128: call TaoSetHessianRoutine(tao,H,H,ComputeHessian, &
129: & 0,ierr)
131: ! Set initial guess
132: call FormInitialGuess(x,ierr)
133: call TaoSetInitialVector(tao,x,ierr)
135: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
136: & '-testmonitor',flg,ierr)
137: if (flg) then
138: call TaoSetMonitor(tao,Monitor,dummy,PETSC_NULL_FUNCTION, &
139: & ierr)
140: endif
142: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
143: & '-testconvergence',flg, ierr)
144: if (flg) then
145: call TaoSetConvergenceTest(tao,ConvergenceTest,dummy, &
146: & ierr)
147: endif
149: ! Check for any TAO command line options
150: call TaoSetFromOptions(tao,ierr)
153: ! SOLVE THE APPLICATION
154: call TaoSolve(tao,ierr)
156: ! Free TAO data structures
157: call TaoDestroy(tao,ierr)
160: ! Free PETSc data structures
161: call VecDestroy(x,ierr)
162: call VecDestroy(localX,ierr)
163: call MatDestroy(H,ierr)
164: call DMDestroy(dm,ierr)
167: ! Finalize TAO and PETSc
168: call PetscFinalize(ierr)
169: end
172: ! ---------------------------------------------------------------------
173: !
174: ! FormInitialGuess - Computes an initial approximation to the solution.
175: !
176: ! Input Parameters:
177: ! X - vector
178: !
179: ! Output Parameters:
180: ! X - vector
181: ! ierr - error code
182: !
183: subroutine FormInitialGuess(X,ierr)
184: use mymodule
185: implicit none
187: ! Input/output variables:
188: Vec X
189: PetscErrorCode ierr
191: ! Local variables:
192: PetscInt i, j, k, xe, ye
193: PetscReal temp, val, hx, hy
194: PetscInt xs, ys, xm, ym
195: PetscInt gxm, gym, gxs, gys
196: PetscInt i1
198: i1 = 1
199: hx = 1.0/real(mx + 1)
200: hy = 1.0/real(my + 1)
202: ! Get corner information
203: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
204: & PETSC_NULL_INTEGER,ierr)
205: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
206: & gxm,gym,PETSC_NULL_INTEGER,ierr)
210: ! Compute initial guess over locally owned part of mesh
211: xe = xs+xm
212: ye = ys+ym
213: do j=ys,ye-1
214: temp = min(j+1,my-j)*hy
215: do i=xs,xe-1
216: k = (j-gys)*gxm + i-gxs
217: val = min((min(i+1,mx-i))*hx,temp)
218: call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,ierr)
219: end do
220: end do
221: call VecAssemblyBegin(X,ierr)
222: call VecAssemblyEnd(X,ierr)
223: return
224: end
227: ! ---------------------------------------------------------------------
228: !
229: ! FormFunctionGradient - Evaluates gradient G(X).
230: !
231: ! Input Parameters:
232: ! tao - the Tao context
233: ! X - input vector
234: ! dummy - optional user-defined context (not used here)
235: !
236: ! Output Parameters:
237: ! f - the function value at X
238: ! G - vector containing the newly evaluated gradient
239: ! ierr - error code
240: !
241: ! Notes:
242: ! This routine serves as a wrapper for the lower-level routine
243: ! "ApplicationGradient", where the actual computations are
244: ! done using the standard Fortran style of treating the local
245: ! input vector data as an array over the local mesh.
246: !
247: subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
248: use mymodule
249: implicit none
251: ! Input/output variables:
252: Tao tao
253: Vec X, G
254: PetscReal f
255: PetscErrorCode ierr
256: PetscInt dummy
258: ! Declarations for use with local array:
261: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
262: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
263: ! will return an array of doubles referenced by x_array offset by x_index.
264: ! i.e., to reference the kth element of X, use x_array(k + x_index).
265: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
266: PetscReal lx_v(0:1)
267: PetscOffset lx_i
269: ! Local variables:
270: PetscReal zero, p5, area, cdiv3
271: PetscReal val, flin, fquad,floc
272: PetscReal v, vb, vl, vr, vt, dvdx
273: PetscReal dvdy, hx, hy
274: PetscInt xe, ye, xsm, ysm
275: PetscInt xep, yep, i, j, k, ind
276: PetscInt xs, ys, xm, ym
277: PetscInt gxs, gys, gxm, gym
278: PetscInt i1
280: i1 = 1
281: 0
282: cdiv3 = param/3.0
283: zero = 0.0
284: p5 = 0.5
285: hx = 1.0/real(mx + 1)
286: hy = 1.0/real(my + 1)
287: fquad = zero
288: flin = zero
290: ! Initialize gradient to zero
291: call VecSet(G,zero,ierr)
293: ! Scatter ghost points to local vector
294: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
295: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
298: ! Get corner information
299: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
300: & PETSC_NULL_INTEGER,ierr)
301: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
302: & gxm,gym,PETSC_NULL_INTEGER,ierr)
304: ! Get pointer to vector data.
305: call VecGetArray(localX,lx_v,lx_i,ierr)
308: ! Set local loop dimensions
309: xe = xs+xm
310: ye = ys+ym
311: if (xs .eq. 0) then
312: xsm = xs-1
313: else
314: xsm = xs
315: endif
316: if (ys .eq. 0) then
317: ysm = ys-1
318: else
319: ysm = ys
320: endif
321: if (xe .eq. mx) then
322: xep = xe+1
323: else
324: xep = xe
325: endif
326: if (ye .eq. my) then
327: yep = ye+1
328: else
329: yep = ye
330: endif
332: ! Compute local gradient contributions over the lower triangular elements
334: do j = ysm, ye-1
335: do i = xsm, xe-1
336: k = (j-gys)*gxm + i-gxs
337: v = zero
338: vr = zero
339: vt = zero
340: if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
341: if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
342: if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
343: dvdx = (vr-v)/hx
344: dvdy = (vt-v)/hy
345: if (i .ne. -1 .and. j .ne. -1) then
346: ind = k
347: val = - dvdx/hx - dvdy/hy - cdiv3
348: call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,ierr)
349: endif
350: if (i .ne. mx-1 .and. j .ne. -1) then
351: ind = k+1
352: val = dvdx/hx - cdiv3
353: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
354: endif
355: if (i .ne. -1 .and. j .ne. my-1) then
356: ind = k+gxm
357: val = dvdy/hy - cdiv3
358: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
359: endif
360: fquad = fquad + dvdx*dvdx + dvdy*dvdy
361: flin = flin - cdiv3 * (v+vr+vt)
362: end do
363: end do
365: ! Compute local gradient contributions over the upper triangular elements
367: do j = ys, yep-1
368: do i = xs, xep-1
369: k = (j-gys)*gxm + i-gxs
370: vb = zero
371: vl = zero
372: v = zero
373: if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
374: if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
375: if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
376: dvdx = (v-vl)/hx
377: dvdy = (v-vb)/hy
378: if (i .ne. mx .and. j .ne. 0) then
379: ind = k-gxm
380: val = - dvdy/hy - cdiv3
381: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
382: endif
383: if (i .ne. 0 .and. j .ne. my) then
384: ind = k-1
385: val = - dvdx/hx - cdiv3
386: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
387: endif
388: if (i .ne. mx .and. j .ne. my) then
389: ind = k
390: val = dvdx/hx + dvdy/hy - cdiv3
391: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,ierr)
392: endif
393: fquad = fquad + dvdx*dvdx + dvdy*dvdy
394: flin = flin - cdiv3*(vb + vl + v)
395: end do
396: end do
398: ! Restore vector
399: call VecRestoreArray(localX,lx_v,lx_i,ierr)
401: ! Assemble gradient vector
402: call VecAssemblyBegin(G,ierr)
403: call VecAssemblyEnd(G,ierr)
405: ! Scale the gradient
406: area = p5*hx*hy
407: floc = area *(p5*fquad+flin)
408: call VecScale(G,area,ierr)
411: ! Sum function contributions from all processes
412: call MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM, &
413: & PETSC_COMM_WORLD,ierr)
414: call PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+ &
415: & 16.0d0*(xep-xs)*(yep-ys),ierr)
416: return
417: end
419: subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
420: use mymodule
421: implicit none
423: Tao tao
424: Vec X
425: Mat H,Hpre
426: PetscErrorCode ierr
427: PetscInt dummy
430: PetscInt i,j,k
431: PetscInt col(0:4),row
432: PetscInt xs,xm,gxs,gxm
433: PetscInt ys,ym,gys,gym
434: PetscReal v(0:4)
435: PetscInt i1
437: i1 = 1
439: ! Get local grid boundaries
440: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
441: & PETSC_NULL_INTEGER,ierr)
442: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
443: & PETSC_NULL_INTEGER,ierr)
445: do j=ys,ys+ym-1
446: do i=xs,xs+xm-1
447: row = (j-gys)*gxm + (i-gxs)
449: k = 0
450: if (j .gt. gys) then
451: v(k) = -1.0
452: col(k) = row-gxm
453: k = k + 1
454: endif
456: if (i .gt. gxs) then
457: v(k) = -1.0
458: col(k) = row - 1
459: k = k +1
460: endif
462: v(k) = 4.0
463: col(k) = row
464: k = k + 1
466: if (i+1 .lt. gxs + gxm) then
467: v(k) = -1.0
468: col(k) = row + 1
469: k = k + 1
470: endif
472: if (j+1 .lt. gys + gym) then
473: v(k) = -1.0
474: col(k) = row + gxm
475: k = k + 1
476: endif
478: call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,ierr)
479: enddo
480: enddo
483: ! Assemble matrix
484: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
485: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
488: ! Tell the matrix we will never add a new nonzero location to the
489: ! matrix. If we do it will generate an error.
491: call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)
492: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
495: call PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr)
497: 0
498: return
499: end
503: subroutine Monitor(tao, dummy, ierr)
504: use mymodule
505: implicit none
507: Tao tao
508: PetscInt dummy
509: PetscErrorCode ierr
511: PetscInt its
512: PetscReal f,gnorm,cnorm,xdiff
513: TaoConvergedReason reason
515: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
516: & reason,ierr)
517: if (mod(its,5) .ne. 0) then
518: call PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n', &
519: & ierr)
520: endif
522: 0
524: return
525: end
527: subroutine ConvergenceTest(tao, dummy, ierr)
528: use mymodule
529: implicit none
531: Tao tao
532: PetscInt dummy
533: PetscErrorCode ierr
535: PetscInt its
536: PetscReal f,gnorm,cnorm,xdiff
537: TaoConvergedReason reason
539: call TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff, &
540: & reason,ierr)
541: if (its .eq. 7) then
542: call TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr)
543: endif
545: 0
547: return
548: end
550: !/*TEST
551: !
552: ! build:
553: ! requires: !complex
554: !
555: ! test:
556: ! args: -tao_smonitor -tao_type nls -tao_gttol 1.e-2
557: !
558: ! test:
559: ! suffix: 2
560: ! nsize: 2
561: ! args: -tao_smonitor -tao_type lmvm -tao_gttol 1.e-2
562: !
563: ! test:
564: ! suffix: 3
565: ! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
566: !TEST*/