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, intent(out) :: 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
60: hx = 1.0/real(mx + 1)
61: hy = 1.0/real(my + 1)
63: ! Get corner information
64: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
65: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
67: ! Compute initial guess over locally owned part of mesh
68: xe = xs + xm
69: ye = ys + ym
70: do j = ys, ye - 1
71: temp = min(j + 1, my - j)*hy
72: do i = xs, xe - 1
73: k = (j - gys)*gxm + i - gxs
74: val = min((min(i + 1, mx - i))*hx, temp)
75: PetscCall(VecSetValuesLocal(X, 1_PETSC_INT_KIND, [k], [val], ADD_VALUES, ierr))
76: end do
77: end do
78: PetscCall(VecAssemblyBegin(X, ierr))
79: PetscCall(VecAssemblyEnd(X, ierr))
80: end
82: ! ---------------------------------------------------------------------
83: !
84: ! FormFunctionGradient - Evaluates gradient G(X).
85: !
86: ! Input Parameters:
87: ! tao - the Tao context
88: ! X - input vector
89: ! dummy - optional user-defined context (not used here)
90: !
91: ! Output Parameters:
92: ! f - the function value at X
93: ! G - vector containing the newly evaluated gradient
94: ! ierr - error code
95: !
96: ! Notes:
97: ! This routine serves as a wrapper for the lower-level routine
98: ! "ApplicationGradient", where the actual computations are
99: ! done using the standard Fortran style of treating the local
100: ! input vector data as an array over the local mesh.
101: !
102: subroutine FormFunctionGradient(ta, X, f, G, dummy, ierr)
103: ! Input/output variables:
104: type(tTao) ta
105: type(tVec) X, G
106: PetscReal f
107: PetscErrorCode, intent(out) :: ierr
108: PetscInt dummy
110: ! Declarations for use with local array:
112: PetscReal, pointer :: lx_v(:)
114: ! Local variables:
115: PetscReal area, cdiv3
116: PetscReal val, flin, fquad, floc
117: PetscReal v, vb, vl, vr, vt, dvdx
118: PetscReal dvdy, hx, hy
119: PetscInt xe, ye, xsm, ysm
120: PetscInt xep, yep, i, j, k, ind
121: PetscInt xs, ys, xm, ym
122: PetscInt gxs, gys, gxm, gym
124: ierr = 0
125: cdiv3 = param/3.0
126: hx = 1.0/real(mx + 1)
127: hy = 1.0/real(my + 1)
128: fquad = 0.0_PETSC_REAL_KIND
129: flin = 0.0_PETSC_REAL_KIND
131: ! Initialize gradient to zero
132: PetscCall(VecSet(G, 0.0_PETSC_REAL_KIND, ierr))
134: ! Scatter ghost points to local vector
135: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
136: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))
138: ! Get corner information
139: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
140: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
142: ! Get pointer to vector data.
143: PetscCall(VecGetArrayRead(localX, lx_v, ierr))
145: ! Set local loop dimensions
146: xe = xs + xm
147: ye = ys + ym
148: if (xs == 0) then
149: xsm = xs - 1
150: else
151: xsm = xs
152: end if
153: if (ys == 0) then
154: ysm = ys - 1
155: else
156: ysm = ys
157: end if
158: if (xe == mx) then
159: xep = xe + 1
160: else
161: xep = xe
162: end if
163: if (ye == my) then
164: yep = ye + 1
165: else
166: yep = ye
167: end if
169: ! Compute local gradient contributions over the lower triangular elements
171: do j = ysm, ye - 1
172: do i = xsm, xe - 1
173: k = (j - gys)*gxm + i - gxs
174: v = 0.0
175: vr = 0.0
176: vt = 0.0
177: if (i >= 0 .and. j >= 0) v = lx_v(k + 1)
178: if (i < mx - 1 .and. j > -1) vr = lx_v(k + 2)
179: if (i > -1 .and. j < my - 1) vt = lx_v(k + 1 + gxm)
180: dvdx = (vr - v)/hx
181: dvdy = (vt - v)/hy
182: if (i /= -1 .and. j /= -1) then
183: ind = k
184: val = -dvdx/hx - dvdy/hy - cdiv3
185: PetscCall(VecSetValuesLocal(G, 1_PETSC_INT_KIND, [k], [val], ADD_VALUES, ierr))
186: end if
187: if (i /= mx - 1 .and. j /= -1) then
188: ind = k + 1
189: val = dvdx/hx - cdiv3
190: PetscCall(VecSetValuesLocal(G, 1_PETSC_INT_KIND, [ind], [val], ADD_VALUES, ierr))
191: end if
192: if (i /= -1 .and. j /= my - 1) then
193: ind = k + gxm
194: val = dvdy/hy - cdiv3
195: PetscCall(VecSetValuesLocal(G, 1_PETSC_INT_KIND, [ind], [val], ADD_VALUES, ierr))
196: end if
197: fquad = fquad + dvdx*dvdx + dvdy*dvdy
198: flin = flin - cdiv3*(v + vr + vt)
199: end do
200: end do
202: ! Compute local gradient contributions over the upper triangular elements
204: do j = ys, yep - 1
205: do i = xs, xep - 1
206: k = (j - gys)*gxm + i - gxs
207: vb = 0.0
208: vl = 0.0
209: v = 0.0
210: if (i < mx .and. j > 0) vb = lx_v(k + 1 - gxm)
211: if (i > 0 .and. j < my) vl = lx_v(k)
212: if (i < mx .and. j < my) v = lx_v(1 + k)
213: dvdx = (v - vl)/hx
214: dvdy = (v - vb)/hy
215: if (i /= mx .and. j /= 0) then
216: ind = k - gxm
217: val = -dvdy/hy - cdiv3
218: PetscCall(VecSetValuesLocal(G, 1_PETSC_INT_KIND, [ind], [val], ADD_VALUES, ierr))
219: end if
220: if (i /= 0 .and. j /= my) then
221: ind = k - 1
222: val = -dvdx/hx - cdiv3
223: PetscCall(VecSetValuesLocal(G, 1_PETSC_INT_KIND, [ind], [val], ADD_VALUES, ierr))
224: end if
225: if (i /= mx .and. j /= my) then
226: ind = k
227: val = dvdx/hx + dvdy/hy - cdiv3
228: PetscCall(VecSetValuesLocal(G, 1_PETSC_INT_KIND, [ind], [val], ADD_VALUES, ierr))
229: end if
230: fquad = fquad + dvdx*dvdx + dvdy*dvdy
231: flin = flin - cdiv3*(vb + vl + v)
232: end do
233: end do
235: ! Restore vector
236: PetscCall(VecRestoreArrayRead(localX, lx_v, ierr))
238: ! Assemble gradient vector
239: PetscCall(VecAssemblyBegin(G, ierr))
240: PetscCall(VecAssemblyEnd(G, ierr))
242: ! Scale the gradient
243: area = .5*hx*hy
244: floc = area*(.5*fquad + flin)
245: PetscCall(VecScale(G, area, ierr))
247: ! Sum function contributions from all processes
248: PetscCallMPI(MPI_Allreduce(floc, f, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD, ierr))
249: PetscCall(PetscLogFlops(20.0d0*(ye - ysm)*(xe - xsm) + 16.0d0*(xep - xs)*(yep - ys), ierr))
250: end
252: subroutine ComputeHessian(ta, X, H, Hpre, dummy, ierr)
253: type(tTao) ta
254: type(tVec) X
255: type(tMat) H, Hpre
256: PetscErrorCode, intent(out) :: ierr
257: PetscInt dummy
259: PetscInt i, j, k
260: PetscInt col(0:4), row
261: PetscInt xs, xm, gxs, gxm
262: PetscInt ys, ym, gys, gym
263: PetscReal v(0:4)
265: ! Get local grid boundaries
266: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
267: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
269: do j = ys, ys + ym - 1
270: do i = xs, xs + xm - 1
271: row = (j - gys)*gxm + (i - gxs)
273: k = 0
274: if (j > gys) then
275: v(k) = -1.0
276: col(k) = row - gxm
277: k = k + 1
278: end if
280: if (i > gxs) then
281: v(k) = -1.0
282: col(k) = row - 1
283: k = k + 1
284: end if
286: v(k) = 4.0
287: col(k) = row
288: k = k + 1
290: if (i + 1 < gxs + gxm) then
291: v(k) = -1.0
292: col(k) = row + 1
293: k = k + 1
294: end if
296: if (j + 1 < gys + gym) then
297: v(k) = -1.0
298: col(k) = row + gxm
299: k = k + 1
300: end if
302: PetscCall(MatSetValuesLocal(H, 1_PETSC_INT_KIND, [row], k, col, v, INSERT_VALUES, ierr))
303: end do
304: end do
306: ! Assemble matrix
307: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr))
308: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr))
310: ! Tell the matrix we will never add a new nonzero location to the
311: ! matrix. If we do it will generate an error.
313: PetscCall(MatSetOption(H, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))
314: PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))
316: PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm, ierr))
318: ierr = 0
319: end
321: subroutine Monitor(ta, dummy, ierr)
322: type(tTao) ta
323: PetscInt dummy
324: PetscErrorCode, intent(out) :: ierr
326: PetscInt its
327: PetscReal f, gnorm, cnorm, xdiff
328: TaoConvergedReason reason
330: PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr))
331: if (mod(its, 5) /= 0) then
332: PetscCall(PetscPrintf(PETSC_COMM_WORLD, 'iteration multiple of 5\n', ierr))
333: end if
335: ierr = 0
337: end
339: subroutine ConvergenceTest(ta, dummy, ierr)
340: type(tTao) ta
341: PetscInt dummy
342: PetscErrorCode, intent(out) :: ierr
344: PetscInt its
345: PetscReal f, gnorm, cnorm, xdiff
346: TaoConvergedReason reason
348: PetscCall(TaoGetSolutionStatus(ta, its, f, gnorm, cnorm, xdiff, reason, ierr))
349: if (its == 7) then
350: PetscCall(TaoSetConvergedReason(ta, TAO_DIVERGED_MAXITS, ierr))
351: end if
353: ierr = 0
355: end
356: end module
358: program eptorsion2f
359: use eptorsion2fmodule
360: implicit none
361: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
362: ! Variable declarations
363: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
364: !
365: ! See additional variable declarations in the file eptorsion2f.h
366: !
367: PetscErrorCode ierr ! used to check for functions returning nonzeros
368: type(tVec) x ! solution vector
369: type(tMat) H ! hessian matrix
370: PetscInt Nx, Ny ! number of processes in x- and y- directions
371: type(tTao) ta ! Tao solver context
372: PetscBool flg
373: PetscInt dummy
375: ! Initialize TAO, PETSc contexts
376: PetscCallA(PetscInitialize(ierr))
377: Nx = PETSC_DECIDE
378: Ny = PETSC_DECIDE
380: ! Specify default parameters and
381: ! check for any command line arguments that might override defaults
382: mx = 10
383: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
384: my = 10
385: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))
386: param = 5.0
387: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', param, flg, ierr))
389: ! Set up distributed array and vectors
390: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, Nx, Ny, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
391: PetscCallA(DMSetFromOptions(dm, ierr))
392: PetscCallA(DMSetUp(dm, ierr))
394: ! Create vectors
395: PetscCallA(DMCreateGlobalVector(dm, x, ierr))
396: PetscCallA(DMCreateLocalVector(dm, localX, ierr))
398: ! Create Hessian
399: PetscCallA(DMCreateMatrix(dm, H, ierr))
400: PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))
402: ! The TAO code begins here
404: ! Create TAO solver
405: PetscCallA(TaoCreate(PETSC_COMM_WORLD, ta, ierr))
406: PetscCallA(TaoSetType(ta, TAOCG, ierr))
408: ! Set routines for function and gradient evaluation
409: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
410: PetscCallA(TaoSetHessian(ta, H, H, ComputeHessian, 0, ierr))
412: ! Set initial guess
413: PetscCallA(FormInitialGuess(x, ierr))
414: PetscCallA(TaoSetSolution(ta, x, ierr))
416: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testmonitor', flg, ierr))
417: if (flg) then
418: PetscCallA(TaoMonitorSet(ta, Monitor, dummy, PETSC_NULL_FUNCTION, ierr))
419: end if
421: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-testconvergence', flg, ierr))
422: if (flg) then
423: PetscCallA(TaoSetConvergenceTest(ta, ConvergenceTest, dummy, ierr))
424: end if
426: ! Check for any TAO command line options
427: PetscCallA(TaoSetFromOptions(ta, ierr))
429: ! SOLVE THE APPLICATION
430: PetscCallA(TaoSolve(ta, ierr))
432: ! Free TAO data structures
433: PetscCallA(TaoDestroy(ta, ierr))
435: ! Free PETSc data structures
436: PetscCallA(VecDestroy(x, ierr))
437: PetscCallA(VecDestroy(localX, ierr))
438: PetscCallA(MatDestroy(H, ierr))
439: PetscCallA(DMDestroy(dm, ierr))
441: ! Finalize TAO and PETSc
442: PetscCallA(PetscFinalize(ierr))
443: end
445: !/*TEST
446: !
447: ! build:
448: ! requires: !complex
449: !
450: ! test:
451: ! args: -tao_monitor_short -tao_type nls -tao_gttol 1.e-2
452: !
453: ! test:
454: ! suffix: 2
455: ! nsize: 2
456: ! args: -tao_monitor_short -tao_type lmvm -tao_gttol 1.e-2
457: !
458: ! test:
459: ! suffix: 3
460: ! args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
461: !TEST*/