Actual source code: eptorsion2f.F90
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: ! ----------------------------------------------------------------------
25: #include "petsc/finclude/petscdmda.h"
26: #include "petsc/finclude/petsctao.h"
27: module eptorsion2fmodule
28: use petscdmda
29: use petsctao
30: implicit none
32: type(tVec) localX
33: type(tDM) dm
34: PetscReal param
35: PetscInt mx, my
37: contains
38: ! ---------------------------------------------------------------------
39: !
40: ! FormInitialGuess - Computes an initial approximation to the solution.
41: !
42: ! Input Parameters:
43: ! X - vector
44: !
45: ! Output Parameters:
46: ! X - vector
47: ! ierr - error code
48: !
49: subroutine FormInitialGuess(X, ierr)
50: ! Input/output variables:
51: Vec X
52: PetscErrorCode ierr
54: ! Local variables:
55: PetscInt i, j, k, xe, ye
56: PetscReal temp, val, hx, hy
57: PetscInt xs, ys, xm, ym
58: PetscInt gxm, gym, gxs, gys
59: PetscInt i1
61: i1 = 1
62: hx = 1.0/real(mx + 1)
63: hy = 1.0/real(my + 1)
65: ! Get corner information
66: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
67: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
69: ! Compute initial guess over locally owned part of mesh
70: xe = xs + xm
71: ye = ys + ym
72: do j = ys, ye - 1
73: temp = min(j + 1, my - j)*hy
74: do i = xs, xe - 1
75: k = (j - gys)*gxm + i - gxs
76: val = min((min(i + 1, mx - i))*hx, temp)
77: PetscCall(VecSetValuesLocal(X, i1, [k], [val], ADD_VALUES, ierr))
78: end do
79: end do
80: PetscCall(VecAssemblyBegin(X, ierr))
81: PetscCall(VecAssemblyEnd(X, ierr))
82: end
84: ! ---------------------------------------------------------------------
85: !
86: ! FormFunctionGradient - Evaluates gradient G(X).
87: !
88: ! Input Parameters:
89: ! tao - the Tao context
90: ! X - input vector
91: ! dummy - optional user-defined context (not used here)
92: !
93: ! Output Parameters:
94: ! f - the function value at X
95: ! G - vector containing the newly evaluated gradient
96: ! ierr - error code
97: !
98: ! Notes:
99: ! This routine serves as a wrapper for the lower-level routine
100: ! "ApplicationGradient", where the actual computations are
101: ! done using the standard Fortran style of treating the local
102: ! input vector data as an array over the local mesh.
103: !
104: subroutine FormFunctionGradient(ta, X, f, G, dummy, ierr)
105: ! Input/output variables:
106: type(tTao) ta
107: type(tVec) X, G
108: PetscReal f
109: PetscErrorCode ierr
110: PetscInt dummy
112: ! Declarations for use with local array:
114: PetscReal, pointer :: lx_v(:)
116: ! Local variables:
117: PetscReal zero, p5, area, cdiv3
118: PetscReal val, flin, fquad, floc
119: PetscReal v, vb, vl, vr, vt, dvdx
120: PetscReal dvdy, hx, hy
121: PetscInt xe, ye, xsm, ysm
122: PetscInt xep, yep, i, j, k, ind
123: PetscInt xs, ys, xm, ym
124: PetscInt gxs, gys, gxm, gym
125: PetscInt i1
127: i1 = 1
128: ierr = 0
129: cdiv3 = param/3.0
130: zero = 0.0
131: p5 = 0.5
132: hx = 1.0/real(mx + 1)
133: hy = 1.0/real(my + 1)
134: fquad = zero
135: flin = zero
137: ! Initialize gradient to zero
138: PetscCall(VecSet(G, zero, ierr))
140: ! Scatter ghost points to local vector
141: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
142: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))
144: ! Get corner information
145: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
146: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
148: ! Get pointer to vector data.
149: PetscCall(VecGetArrayRead(localX, lx_v, ierr))
151: ! Set local loop dimensions
152: xe = xs + xm
153: ye = ys + ym
154: if (xs == 0) then
155: xsm = xs - 1
156: else
157: xsm = xs
158: end if
159: if (ys == 0) then
160: ysm = ys - 1
161: else
162: ysm = ys
163: end if
164: if (xe == mx) then
165: xep = xe + 1
166: else
167: xep = xe
168: end if
169: if (ye == my) then
170: yep = ye + 1
171: else
172: yep = ye
173: end if
175: ! Compute local gradient contributions over the lower triangular elements
177: do j = ysm, ye - 1
178: do i = xsm, xe - 1
179: k = (j - gys)*gxm + i - gxs
180: v = zero
181: vr = zero
182: vt = zero
183: if (i >= 0 .and. j >= 0) v = lx_v(k + 1)
184: if (i < mx - 1 .and. j > -1) vr = lx_v(k + 2)
185: if (i > -1 .and. j < my - 1) vt = lx_v(k + 1 + gxm)
186: dvdx = (vr - v)/hx
187: dvdy = (vt - v)/hy
188: if (i /= -1 .and. j /= -1) then
189: ind = k
190: val = -dvdx/hx - dvdy/hy - cdiv3
191: PetscCall(VecSetValuesLocal(G, i1, [k], [val], ADD_VALUES, ierr))
192: end if
193: if (i /= mx - 1 .and. j /= -1) then
194: ind = k + 1
195: val = dvdx/hx - cdiv3
196: PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr))
197: end if
198: if (i /= -1 .and. j /= my - 1) then
199: ind = k + gxm
200: val = dvdy/hy - cdiv3
201: PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr))
202: end if
203: fquad = fquad + dvdx*dvdx + dvdy*dvdy
204: flin = flin - cdiv3*(v + vr + vt)
205: end do
206: end do
208: ! Compute local gradient contributions over the upper triangular elements
210: do j = ys, yep - 1
211: do i = xs, xep - 1
212: k = (j - gys)*gxm + i - gxs
213: vb = zero
214: vl = zero
215: v = zero
216: if (i < mx .and. j > 0) vb = lx_v(k + 1 - gxm)
217: if (i > 0 .and. j < my) vl = lx_v(k)
218: if (i < mx .and. j < my) v = lx_v(1 + k)
219: dvdx = (v - vl)/hx
220: dvdy = (v - vb)/hy
221: if (i /= mx .and. j /= 0) then
222: ind = k - gxm
223: val = -dvdy/hy - cdiv3
224: PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr))
225: end if
226: if (i /= 0 .and. j /= my) then
227: ind = k - 1
228: val = -dvdx/hx - cdiv3
229: PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr))
230: end if
231: if (i /= mx .and. j /= my) then
232: ind = k
233: val = dvdx/hx + dvdy/hy - cdiv3
234: PetscCall(VecSetValuesLocal(G, i1, [ind], [val], ADD_VALUES, ierr))
235: end if
236: fquad = fquad + dvdx*dvdx + dvdy*dvdy
237: flin = flin - cdiv3*(vb + vl + v)
238: end do
239: end do
241: ! Restore vector
242: PetscCall(VecRestoreArrayRead(localX, lx_v, ierr))
244: ! Assemble gradient vector
245: PetscCall(VecAssemblyBegin(G, ierr))
246: PetscCall(VecAssemblyEnd(G, ierr))
248: ! Scale the gradient
249: area = p5*hx*hy
250: floc = area*(p5*fquad + flin)
251: PetscCall(VecScale(G, area, ierr))
253: ! Sum function contributions from all processes
254: PetscCallMPI(MPI_Allreduce(floc, f, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD, ierr))
255: PetscCall(PetscLogFlops(20.0d0*(ye - ysm)*(xe - xsm) + 16.0d0*(xep - xs)*(yep - ys), ierr))
256: end
258: subroutine ComputeHessian(ta, X, H, Hpre, dummy, ierr)
259: type(tTao) ta
260: type(tVec) X
261: type(tMat) H, Hpre
262: PetscErrorCode ierr
263: PetscInt dummy
265: PetscInt i, j, k
266: PetscInt col(0:4), row
267: PetscInt xs, xm, gxs, gxm
268: PetscInt ys, ym, gys, gym
269: PetscReal v(0:4)
270: PetscInt i1
272: i1 = 1
274: ! Get local grid boundaries
275: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
276: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
278: do j = ys, ys + ym - 1
279: do i = xs, xs + xm - 1
280: row = (j - gys)*gxm + (i - gxs)
282: k = 0
283: if (j > gys) then
284: v(k) = -1.0
285: col(k) = row - gxm
286: k = k + 1
287: end if
289: if (i > gxs) then
290: v(k) = -1.0
291: col(k) = row - 1
292: k = k + 1
293: end if
295: v(k) = 4.0
296: col(k) = row
297: k = k + 1
299: if (i + 1 < gxs + gxm) then
300: v(k) = -1.0
301: col(k) = row + 1
302: k = k + 1
303: end if
305: if (j + 1 < gys + gym) then
306: v(k) = -1.0
307: col(k) = row + gxm
308: k = k + 1
309: end if
311: PetscCall(MatSetValuesLocal(H, i1, [row], k, col, v, INSERT_VALUES, ierr))
312: end do
313: end do
315: ! Assemble matrix
316: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr))
317: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr))
319: ! Tell the matrix we will never add a new nonzero location to the
320: ! matrix. If we do it will generate an error.
322: PetscCall(MatSetOption(H, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))
323: PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))
325: PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm, ierr))
327: ierr = 0
328: end
330: subroutine Monitor(ta, dummy, ierr)
331: type(tTao) ta
332: PetscInt dummy
333: PetscErrorCode ierr
335: PetscInt its
336: PetscReal f, gnorm, cnorm, xdiff
337: TaoConvergedReason reason
339: PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr))
340: if (mod(its, 5) /= 0) then
341: PetscCall(PetscPrintf(PETSC_COMM_WORLD, 'iteration multiple of 5\n', ierr))
342: end if
344: ierr = 0
346: end
348: subroutine ConvergenceTest(ta, dummy, ierr)
349: type(tTao) ta
350: PetscInt dummy
351: PetscErrorCode ierr
353: PetscInt its
354: PetscReal f, gnorm, cnorm, xdiff
355: TaoConvergedReason reason
357: PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr))
358: if (its == 7) then
359: PetscCall(TaoSetConvergedReason(ta, TAO_DIVERGED_MAXITS, ierr))
360: end if
362: ierr = 0
364: end
365: end module
367: program eptorsion2f
368: use eptorsion2fmodule
369: implicit none
370: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371: ! Variable declarations
372: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
373: !
374: ! See additional variable declarations in the file eptorsion2f.h
375: !
376: PetscErrorCode ierr ! used to check for functions returning nonzeros
377: type(tVec) x ! solution vector
378: type(tMat) H ! hessian matrix
379: PetscInt Nx, Ny ! number of processes in x- and y- directions
380: type(tTao) ta ! Tao solver context
381: PetscBool flg
382: PetscInt i1
383: PetscInt dummy
385: i1 = 1
387: ! Initialize TAO, PETSc contexts
388: PetscCallA(PetscInitialize(ierr))
390: ! Specify default parameters
391: param = 5.0
392: mx = 10
393: my = 10
394: Nx = PETSC_DECIDE
395: Ny = PETSC_DECIDE
397: ! Check for any command line arguments that might override defaults
398: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
399: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))
400: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', param, flg, ierr))
402: ! Set up distributed array and vectors
403: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, Nx, Ny, i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
404: PetscCallA(DMSetFromOptions(dm, ierr))
405: PetscCallA(DMSetUp(dm, ierr))
407: ! Create vectors
408: PetscCallA(DMCreateGlobalVector(dm, x, ierr))
409: PetscCallA(DMCreateLocalVector(dm, localX, ierr))
411: ! Create Hessian
412: PetscCallA(DMCreateMatrix(dm, H, ierr))
413: PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))
415: ! The TAO code begins here
417: ! Create TAO solver
418: PetscCallA(TaoCreate(PETSC_COMM_WORLD, ta, ierr))
419: PetscCallA(TaoSetType(ta, TAOCG, ierr))
421: ! Set routines for function and gradient evaluation
423: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
424: PetscCallA(TaoSetHessian(ta, H, H, ComputeHessian, 0, ierr))
426: ! Set initial guess
427: PetscCallA(FormInitialGuess(x, ierr))
428: PetscCallA(TaoSetSolution(ta, x, ierr))
430: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testmonitor', flg, ierr))
431: if (flg) then
432: PetscCallA(TaoMonitorSet(ta, Monitor, dummy, PETSC_NULL_FUNCTION, ierr))
433: end if
435: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testconvergence', flg, ierr))
436: if (flg) then
437: PetscCallA(TaoSetConvergenceTest(ta, ConvergenceTest, dummy, ierr))
438: end if
440: ! Check for any TAO command line options
441: PetscCallA(TaoSetFromOptions(ta, ierr))
443: ! SOLVE THE APPLICATION
444: PetscCallA(TaoSolve(ta, ierr))
446: ! Free TAO data structures
447: PetscCallA(TaoDestroy(ta, ierr))
449: ! Free PETSc data structures
450: PetscCallA(VecDestroy(x, ierr))
451: PetscCallA(VecDestroy(localX, ierr))
452: PetscCallA(MatDestroy(H, ierr))
453: PetscCallA(DMDestroy(dm, ierr))
455: ! Finalize TAO and PETSc
456: PetscCallA(PetscFinalize(ierr))
457: end
459: !/*TEST
460: !
461: ! build:
462: ! requires: !complex
463: !
464: ! test:
465: ! args: -tao_monitor_short -tao_type nls -tao_gttol 1.e-2
466: !
467: ! test:
468: ! suffix: 2
469: ! nsize: 2
470: ! args: -tao_monitor_short -tao_type lmvm -tao_gttol 1.e-2
471: !
472: ! test:
473: ! suffix: 3
474: ! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
475: !TEST*/