Actual source code: eptorsion2f.F
1: ! "$Id$";
2: !
3: ! Program usage: mpirun -np <proc> eptorsion2f [all TAO options]
4: !
5: ! Description: This example demonstrates use of the TAO package to solve
6: ! unconstrained minimization problems in parallel. This example is based
7: ! on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.
8: ! The command line options are:
9: ! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
10: ! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
11: ! -par <param>, where <param> = angle of twist per unit length
12: !
13: !/*T
14: ! Concepts: TAO - Solving an unconstrained minimization problem
15: ! Routines: TaoInitialize(); TaoFinalize();
16: ! Routines: TaoCreate(); TaoDestroy();
17: ! Routines: TaoApplicationCreate(); TaoAppDestroy();
18: ! Routines: TaoAppSetObjectiveAndGradientRoutine();
19: ! Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
20: ! Routines: TaoSetOptions();
21: ! Routines: TaoAppSetInitialSolutionVec(); TaoSolveApplication(); TaoDestroy();
22: ! Routines: TaoGetSolutionStatus();
23: ! Processors: n
24: !T*/
25: !
26: ! ----------------------------------------------------------------------
27: !
28: ! Elastic-plastic torsion problem.
29: !
30: ! The elastic plastic torsion problem arises from the determination
31: ! of the stress field on an infinitely long cylindrical bar, which is
32: ! equivalent to the solution of the following problem:
33: ! min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
34: ! where C is the torsion angle per unit length.
35: !
36: ! The C version of this code is eptorsion2.c
37: !
38: ! ----------------------------------------------------------------------
40: implicit none
41: #include "eptorsion2f.h"
43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: ! Variable declarations
45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: !
47: ! See additional variable declarations in the file eptorsion2f.h
48: !
49: integer info ! used to check for functions returning nonzeros
50: Vec x ! solution vector
51: Mat H ! hessian matrix
52: PetscInt Nx, Ny ! number of processes in x- and y- directions
53: TAO_SOLVER tao ! TAO_SOLVER solver context
54: TAO_APPLICATION torsionapp ! TAO application context (PETSc)
55: TaoTerminateReason reason
56: PetscTruth flg
57: PetscInt iter ! iteration information
58: PetscScalar ff,gnorm,cnorm,xdiff
59: PetscInt i1
60:
62: ! Note: Any user-defined Fortran routines (such as FormGradient)
63: ! MUST be declared as external.
65: external FormInitialGuess,FormFunctionGradient,ComputeHessian
67: i1 = 1
69: ! Initialize TAO, PETSc contexts
70: call PetscInitialize(PETSC_NULL_CHARACTER,info)
71: call TaoInitialize(PETSC_NULL_CHARACTER,info)
73: ! Specify default parameters
74: param = 5.0d0
75: mx = 10
76: my = 10
77: Nx = PETSC_DECIDE
78: Ny = PETSC_DECIDE
80: ! Check for any command line arguments that might override defaults
81: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mx",mx,flg,info)
82: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,info)
83: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-par", &
84: & param,flg,info)
86:
87: ! Set up distributed array and vectors
88: call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX, &
89: & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
90: & da,info)
92: ! Create vectors
93: call DACreateGlobalVector(da,x,info)
94: call DACreateLocalVector(da,localX,info)
96: ! Create Hessian
97: call DAGetMatrix(da,MATAIJ,H,info)
98: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)
100: ! The TAO code begins here
102: ! Create TAO solver
103: call TaoCreate(MPI_COMM_WORLD,'tao_cg',tao,info)
104: call TaoApplicationCreate(MPI_COMM_WORLD,torsionapp,info)
106: ! Set routines for function and gradient evaluation
108: ! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
109: call TaoAppSetObjectiveAndGradientRo(torsionapp, &
110: & FormFunctionGradient,TAO_NULL_OBJECT,info)
111: call TaoAppSetHessianMat(torsionapp,H,H,info)
112: call TaoAppSetHessianRoutine(torsionapp,ComputeHessian, &
113: & TAO_NULL_OBJECT,info)
115: ! Set initial guess
116: call FormInitialGuess(x,info)
117: call TaoAppSetInitialSolutionVec(torsionapp,x,info)
119: ! Check for any TAO command line options
120: call TaoSetOptions(torsionapp, tao,info)
123: ! SOLVE THE APPLICATION
124: call TaoSolveApplication(torsionapp,tao,info)
126: ! Get information on termination
127: call TaoGetSolutionStatus(tao,iter,ff,gnorm,cnorm,xdiff, &
128: & reason,info)
129: if (reason .lt. 0) then
130: print *,'TAO did not terminate successfully'
131: endif
133:
134: ! Free TAO data structures
135: call TaoDestroy(tao,info)
136: call TaoAppDestroy(torsionapp,info);
138:
139: ! Free PETSc data structures
140: call VecDestroy(x,info)
141: call VecDestroy(localX,info)
142: call MatDestroy(H,info)
143: call DADestroy(da,info)
146: ! Finalize TAO and PETSc
147: call PetscFinalize(info)
148: call TaoFinalize(info)
150: end
153: ! ---------------------------------------------------------------------
154: !
155: ! FormInitialGuess - Computes an initial approximation to the solution.
156: !
157: ! Input Parameters:
158: ! X - vector
159: !
160: ! Output Parameters:
161: ! X - vector
162: ! info - error code
163: !
164: subroutine FormInitialGuess(X,info)
165: implicit none
167: ! mx, my are defined in eptorsion2f.h
168: #include "eptorsion2f.h"
170: ! Input/output variables:
171: Vec X
172: integer info
174: ! Local variables:
175: PetscInt i, j, k, xe, ye
176: PetscScalar temp, val, hx, hy
177: PetscInt xs, ys, xm, ym
178: PetscInt gxm, gym, gxs, gys
179: PetscInt i1
181: i1 = 1
182: hx = 1.0d0/(mx + 1)
183: hy = 1.0d0/(my + 1)
185: ! Get corner information
186: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
187: & PETSC_NULL_INTEGER,info)
188: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
189: & gxm,gym,PETSC_NULL_INTEGER,info)
193: ! Compute initial guess over locally owned part of mesh
194: xe = xs+xm
195: ye = ys+ym
196: do j=ys,ye-1
197: temp = min(j+1,my-j)*hy
198: do i=xs,xe-1
199: k = (j-gys)*gxm + i-gxs
200: val = min((min(i+1,mx-i))*hx,temp)
201: call VecSetValuesLocal(X,i1,k,val,ADD_VALUES,info)
202: end do
203: end do
204: call VecAssemblyBegin(X,info)
205: call VecAssemblyEnd(X,info)
206: return
207: end
210: ! ---------------------------------------------------------------------
211: !
212: ! FormFunctionGradient - Evaluates gradient G(X).
213: !
214: ! Input Parameters:
215: ! tao - the TAO_SOLVER context
216: ! X - input vector
217: ! dummy - optional user-defined context (not used here)
218: !
219: ! Output Parameters:
220: ! f - the function value at X
221: ! G - vector containing the newly evaluated gradient
222: ! info - error code
223: !
224: ! Notes:
225: ! This routine serves as a wrapper for the lower-level routine
226: ! "ApplicationGradient", where the actual computations are
227: ! done using the standard Fortran style of treating the local
228: ! input vector data as an array over the local mesh.
229: !
230: subroutine FormFunctionGradient(taoapp,X,f,G,dummy,info)
231: implicit none
233: ! da, mx, my, param, localX declared in eptorsion2f.h
234: #include "eptorsion2f.h"
236: ! Input/output variables:
237: TAO_APPLICATION taoapp
238: Vec X, G
239: PetscScalar f
240: integer info
241: PetscInt dummy
243: ! Declarations for use with local array:
246: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
247: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
248: ! will return an array of doubles referenced by x_array offset by x_index.
249: ! i.e., to reference the kth element of X, use x_array(k + x_index).
250: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
251: PetscScalar lx_v(0:1)
252: PetscOffset lx_i
254: ! Local variables:
255: PetscScalar zero, p5, area, cdiv3
256: PetscScalar val, flin, fquad,floc
257: PetscScalar v, vb, vl, vr, vt, dvdx
258: PetscScalar dvdy, hx, hy
259: PetscInt xe, ye, xsm, ysm
260: PetscInt xep, yep, i, j, k, ind
261: PetscInt xs, ys, xm, ym
262: PetscInt gxs, gys, gxm, gym
263: PetscInt i1
265: i1 = 1
266: info = 0
267: cdiv3 = param/3.0d0
268: zero = 0.0d0
269: p5 = 0.5d0
270: hx = 1.0d0/(mx + 1)
271: hy = 1.0d0/(my + 1)
272: fquad = zero
273: flin = zero
275: ! Initialize gradient to zero
276: call VecSet(G,zero,info)
278: ! Scatter ghost points to local vector
279: call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
280: call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
283: ! Get corner information
284: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
285: & PETSC_NULL_INTEGER,info)
286: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
287: & gxm,gym,PETSC_NULL_INTEGER,info)
289: ! Get pointer to vector data.
290: call VecGetArray(localX,lx_v,lx_i,info)
293: ! Set local loop dimensions
294: xe = xs+xm
295: ye = ys+ym
296: if (xs .eq. 0) then
297: xsm = xs-1
298: else
299: xsm = xs
300: endif
301: if (ys .eq. 0) then
302: ysm = ys-1
303: else
304: ysm = ys
305: endif
306: if (xe .eq. mx) then
307: xep = xe+1
308: else
309: xep = xe
310: endif
311: if (ye .eq. my) then
312: yep = ye+1
313: else
314: yep = ye
315: endif
317: ! Compute local gradient contributions over the lower triangular elements
318:
319: do j = ysm, ye-1
320: do i = xsm, xe-1
321: k = (j-gys)*gxm + i-gxs
322: v = zero
323: vr = zero
324: vt = zero
325: if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
326: if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
327: if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
328: dvdx = (vr-v)/hx
329: dvdy = (vt-v)/hy
330: if (i .ne. -1 .and. j .ne. -1) then
331: ind = k
332: val = - dvdx/hx - dvdy/hy - cdiv3
333: call VecSetValuesLocal(G,i1,k,val,ADD_VALUES,info)
334: endif
335: if (i .ne. mx-1 .and. j .ne. -1) then
336: ind = k+1
337: val = dvdx/hx - cdiv3
338: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,info)
339: endif
340: if (i .ne. -1 .and. j .ne. my-1) then
341: ind = k+gxm
342: val = dvdy/hy - cdiv3
343: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,info)
344: endif
345: fquad = fquad + dvdx*dvdx + dvdy*dvdy
346: flin = flin - cdiv3 * (v+vr+vt)
347: end do
348: end do
350: ! Compute local gradient contributions over the upper triangular elements
352: do j = ys, yep-1
353: do i = xs, xep-1
354: k = (j-gys)*gxm + i-gxs
355: vb = zero
356: vl = zero
357: v = zero
358: if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
359: if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
360: if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
361: dvdx = (v-vl)/hx
362: dvdy = (v-vb)/hy
363: if (i .ne. mx .and. j .ne. 0) then
364: ind = k-gxm
365: val = - dvdy/hy - cdiv3
366: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,info)
367: endif
368: if (i .ne. 0 .and. j .ne. my) then
369: ind = k-1
370: val = - dvdx/hx - cdiv3
371: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,info)
372: endif
373: if (i .ne. mx .and. j .ne. my) then
374: ind = k
375: val = dvdx/hx + dvdy/hy - cdiv3
376: call VecSetValuesLocal(G,i1,ind,val,ADD_VALUES,info)
377: endif
378: fquad = fquad + dvdx*dvdx + dvdy*dvdy
379: flin = flin - cdiv3*(vb + vl + v)
380: end do
381: end do
383: ! Restore vector
384: call VecRestoreArray(localX,lx_v,lx_i,info)
385: if (info .ne. 0) print *,'VecRestoreArray'
386: ! Assemble gradient vector
387: call VecAssemblyBegin(G,info)
388: call VecAssemblyEnd(G,info)
390: ! Scale the gradient
391: area = p5*hx*hy
392: floc = area *(p5*fquad+flin)
393: call VecScale(G,area,info)
396: ! Sum function contributions from all processes
397: call MPI_Allreduce(floc,f,1,MPI_DOUBLE_PRECISION,MPI_SUM, &
398: & MPI_COMM_WORLD,info)
399: if (info .ne. 0) print *,'MPI_Allreduce'
400: call PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16, &
401: & info)
405: return
406: end
411: subroutine ComputeHessian(taoapp, X, H, Hpre, flag, dummy, info)
412: implicit none
413: #include "eptorsion2f.h"
414: TAO_APPLICATION taoapp
415: Vec X
416: Mat H,Hpre
417: MatStructure flag
418: integer info
419: PetscInt dummy
421:
422: PetscInt i,j,k
423: PetscInt col(0:4),row
424: PetscInt xs,xm,gxs,gxm
425: PetscInt ys,ym,gys,gym
426: PetscScalar v(0:4)
427: PetscInt i1
429: i1 = 1
431: ! Get local grid boundaries
432: call DAGetCorners(da,xs,ys,TAO_NULL_INTEGER,xm,ym, &
433: & TAO_NULL_INTEGER,info)
434: call DAGetGhostCorners(da,gxs,gys,TAO_NULL_INTEGER,gxm,gym, &
435: & TAO_NULL_INTEGER,info)
437: do j=ys,ys+ym-1
438: do i=xs,xs+xm-1
439: row = (j-gys)*gxm + (i-gxs)
441: k = 0
442: if (j .gt. gys) then
443: v(k) = -1.0d0
444: col(k) = row-gxm
445: k = k + 1
446: endif
448: if (i .gt. gxs) then
449: v(k) = -1.0d0
450: col(k) = row - 1
451: k = k +1
452: endif
454: v(k) = 4.0d0
455: col(k) = row
456: k = k + 1
458: if (i+1 .lt. gxs + gxm) then
459: v(k) = -1.0d0
460: col(k) = row + 1
461: k = k + 1
462: endif
464: if (j+1 .lt. gys + gym) then
465: v(k) = -1.0d0
466: col(k) = row + gxm
467: k = k + 1
468: endif
470: call MatSetValuesLocal(H,i1,row,k,col,v,INSERT_VALUES,info)
471: enddo
472: enddo
474:
475: ! Assemble matrix
476: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
477: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)
480: ! Tell the matrix we will never add a new nonzero location to the
481: ! matrix. If we do it will generate an error.
483: call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,info)
484: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)
487: call PetscLogFlops(9*xm*ym + 49*xm,info)
489: info = 0
490: return
491: end
492:
493:
494: