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