Actual source code: eptorsion3.c
1: /*$Id: eptorsion.c, v 1.1 2002/08/08 10:30 lopezca@mauddib.mcs.anl.gov $*/
3: /* Program usage: mpirun -np <proc> eptorsion [-help] [all TAO options] */
5: /*
6: Include "tao.h" so we can use TAO solvers.
7: petscda.h for distributed array
8: ad_deriv.h for AD gradient
9: */
11: #include "petscda.h"
12: #include "tao.h"
13: #include "taodaapplication.h"
15: static char help[] = "This example is based on the Elastic-Plastic Torsion (dept)\n\
16: problem from the MINPACK-2 test suite.\n\
17: The command line options are:\n\
18: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
19: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
20: -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
21: -byelement, if computation is made by functions on rectangular elements\n\
22: -adic, if AD is used (AD is not used by default)\n\
23: -u1 <u1>, where <u1> = upper limit in the 1st coordinate direction\n\
24: -u2 <u2>, where <u2> = upper limit in the 2nd coordinate direction\n\
25: -par <param>, where <param> = angle of twist per unit length\n\n";
27: /*T
28: Concepts: TAO - Solving a bounded minimization problem
29: Routines: TaoInitialize(); TaoFinalize();
30: Routines: TaoCreate(); TaoDestroy();
31: Routines: DAApplicationCreate(); DAApplicationDestroy();
32: Routines: DAAppSetVariableBoundsRoutine();
33: Routines: DAAppSetElementObjectiveAndGradientRoutine();
34: Routines: DAAppSetElementHessianRoutine();
35: Routines: DAAppSetObjectiveAndGradientRoutine();
36: Routines: DAAppSetADElementFunctionGradient();
37: Routines: DAAppSetHessianRoutine();
38: Routines: TaoSetOptions();
39: Routines: TaoGetSolutionStatus(); TaoDAAppSolve();
40: Routines: DAAppSetMonitor(); TaoView();
41: Routines: DAAppGetSolution();
42: Routines: DAAppGetInterpolationMatrix();
43: Processors: n
44: T*/
46: /*
47: User-defined application context - contains data needed by the
48: application-provided call-back routines.
49: */
50: typedef struct {
51: InactiveDouble param;
52: InactiveDouble hx, hy; /* increment size in both directions */
53: InactiveDouble area; /* area of the triangles */
54: } ADFGCtx;
56: typedef struct {
57: PetscReal param; /* 'c' parameter */
58: PetscReal u1, u2; /* domain upper limits (lower limits = 0) */
59: double hx, hy; /* increment size in both directions */
60: double area; /* area of the triangles */
61: ADFGCtx fgctx; /* Used only when an ADIC generated gradient is used */
62: } AppCtx;
63: int ad_EPTorsLocalFunction(PetscInt[2], DERIV_TYPE[4], DERIV_TYPE*, void*);
65: /* User-defined routines found in this file */
66: static int AppCtxInitialize(void *ptr);
67: static int FormInitialGuess(DA, Vec);
69: static int EPTorsLocalFunctionGradient(PetscInt[2], double x[4], double *f, double g[4], void *ptr);
70: static int EPTorsLocalHessian(PetscInt[2], double x[4], double H[4][4], void *ptr);
72: static int WholeEPTorsFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);
73: static int WholeEPTorsHessian(TAO_APPLICATION,DA,Vec,Mat,void*);
75: static int DASetBounds(TAO_APPLICATION, DA, Vec, Vec, void*);
77: static int MyGridMonitorBefore(TAO_APPLICATION, DA, PetscInt, void *);
82: int main( int argc, char **argv ) {
84: int info; /* used to check for functions returning nonzeros */
85: PetscInt mx,my,Nx,Ny;
86: double ff,gnorm;
87: PetscInt iter, nlevels; /* multigrid levels */
88: DA DAarray[20];
89: Vec X;
90: PetscTruth flg, PreLoad = PETSC_TRUE; /* flags */
91: TaoMethod method = "tao_gpcg"; /* minimization method */
92: AppCtx user; /* user-defined work context */
93: TAO_SOLVER tao; /* TAO_SOLVER solver context */
94: TAO_APPLICATION EPTorsApp; /* The PETSc application */
95: TaoTerminateReason reason;
96: KSP ksp;
97: PC pc;
99: /* Initialize TAO */
100: PetscInitialize(&argc, &argv, (char *)0, help);
101: TaoInitialize(&argc, &argv, (char *)0, help);
103: PreLoadBegin(PreLoad,"Solve");
104:
105: info = AppCtxInitialize((void*)&user); CHKERRQ(info);
107: nlevels=5;
108: info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
109: mx = my = 11; /* these correspond to 10 segments on each dimension */
110: info = PetscOptionsGetInt(TAO_NULL, "-mx", &mx, &flg); CHKERRQ(info);
111: info = PetscOptionsGetInt(TAO_NULL, "-my", &my, &flg); CHKERRQ(info);
112: if (PreLoadIt == 0) {
113: nlevels = 1; mx = 11; my = 11; }
115: PetscPrintf(MPI_COMM_WORLD,"\n---- Elastic-Plastic Torsion Problem -----\n\n");
117: /* Let PETSc determine the vector distribution */
118: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
120: /* Create distributed array (DA) to manage parallel grid and vectors */
121: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,mx,
122: my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
123: for (iter=1;iter<nlevels;iter++){
124: info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
125: }
127: /* Create TAO solver and set desired solution method */
128: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
129: info = TaoApplicationCreate(PETSC_COMM_WORLD,&EPTorsApp); CHKERRQ(info);
130: info = TaoAppSetDAApp(EPTorsApp, DAarray, nlevels ); CHKERRQ(info);
131: /* Sets routines for function, gradient and bounds evaluation */
132: info = DAAppSetVariableBoundsRoutine(EPTorsApp,DASetBounds,(void *)&user); CHKERRQ(info);
134: info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
135: if (flg) {
137: /* Sets routines for function and gradient evaluation, element by element */
138: info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
139: if (flg) {
140: info = DAAppSetADElementFunctionGradient(EPTorsApp,ad_EPTorsLocalFunction,192,(void *)&user.fgctx); CHKERRQ(info);
141: } else {
142: info = DAAppSetElementObjectiveAndGradientRoutine(EPTorsApp,EPTorsLocalFunctionGradient,42,(void *)&user); CHKERRQ(info);
143: }
144: /* Sets routines for Hessian evaluation, element by element */
145: info = DAAppSetElementHessianRoutine(EPTorsApp,EPTorsLocalHessian,6,(void*)&user); CHKERRQ(info);
147: } else {
149: /* Sets routines for function and gradient evaluation, all in one routine */
150: info = DAAppSetObjectiveAndGradientRoutine(EPTorsApp,WholeEPTorsFunctionGradient,(void *)&user); CHKERRQ(info);
152: /* Sets routines for Hessian evaluation, all in one routine */
153: info = DAAppSetHessianRoutine(EPTorsApp,WholeEPTorsHessian,(void*)&user); CHKERRQ(info);
154:
155: }
157: info = DAAppSetBeforeMonitor(EPTorsApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
158: info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
159: if (flg){
160: info = DAAppPrintInterpolationError(EPTorsApp); CHKERRQ(info);
161: info = DAAppPrintStageTimes(EPTorsApp); CHKERRQ(info);
162: }
163: info = TaoAppSetRelativeTolerance(EPTorsApp,1.0e-6); CHKERRQ(info);
164: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
165: info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);
167: info = TaoAppGetKSP(EPTorsApp,&ksp);CHKERRQ(info);
168: info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
169: info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100);CHKERRQ(info);
170: info = KSPGetPC(ksp,&pc);CHKERRQ(info);
171: info = PCSetType(pc,PCBJACOBI);CHKERRQ(info);
173: /* Check for any tao command line options */
174: info = TaoSetOptions(EPTorsApp, tao); CHKERRQ(info);
176: info = DAAppGetSolution(EPTorsApp,0,&X); CHKERRQ(info);
177: info = FormInitialGuess(DAarray[0],X); CHKERRQ(info);
178: info = DAAppSetInitialSolution(EPTorsApp,X); CHKERRQ(info);
179:
180: /* SOLVE THE APPLICATION */
181: info = TaoDAAppSolve(EPTorsApp, tao); CHKERRQ(info);
183: /* Get information on termination */
184: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
185: if (reason <= 0 ){
186: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
187: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
188: }
190: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
191: if (flg){
192: info = DAAppGetSolution(EPTorsApp,nlevels-1,&X); CHKERRQ(info);
193: info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
194: }
196: /* To View TAO solver information */
197: // info = TaoView(tao); CHKERRQ(info);
199: /* Free TAO data structures */
200: info = TaoDestroy(tao); CHKERRQ(info);
201: info = TaoAppDestroy(EPTorsApp); CHKERRQ(info);
203: /* Free PETSc data structures */
204: for (iter=0;iter<nlevels;iter++){
205: info = DADestroy(DAarray[iter]); CHKERRQ(info);
206: }
208: PreLoadEnd();
210: /* Finalize TAO */
211: TaoFinalize();
212: PetscFinalize();
214: return 0;
215: } /* main */
219: /*----- The following two routines
220: MyGridMonitorBefore MyGridMonitorAfter
221: help diplay info of iterations at every grid level
222: */
225: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx) {
227: AppCtx *user = (AppCtx*)ctx;
228: PetscInt mx,my;
229: int info;
231: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
232: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
233: user->hx = user->u1 / (mx - 1);
234: user->hy = user->u2 / (my - 1);
235: user->area = 0.5 * user->hx * user->hy;
236: user->fgctx.hx = user->hx;
237: user->fgctx.hy = user->hy;
238: user->fgctx.area = user->area;
239: user->fgctx.param = user->param;
241: PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my);
243: return 0;
244: }
247: /*------- USER-DEFINED: initialize the application context information -------*/
250: /*
251: AppCtxInitialize - Sets initial values for the application context parameters
253: Input:
254: ptr - void user-defined application context
256: Output:
257: ptr - user-defined application context with the default or user-provided
258: parameters
259: */
260: static int AppCtxInitialize(void *ptr) {
262: AppCtx *user = (AppCtx*)ptr;
263: PetscTruth flg; /* flag for PETSc calls */
264: int info;
266: /* Specify default parameters */
267: user->param = 25.0;
268: user->u1 = user->u2 = 1.0;
270: /* Check for command line arguments that override defaults */
271: info = PetscOptionsGetReal(TAO_NULL, "-par", &user->param, &flg); CHKERRQ(info);
272: info = PetscOptionsGetReal(TAO_NULL, "-u1", &user->u1, &flg); CHKERRQ(info);
273: info = PetscOptionsGetReal(TAO_NULL, "-u2", &user->u2, &flg); CHKERRQ(info);
275: return 0;
276: } /* AppCtxInitialize */
280: static int FormInitialGuess(DA da, Vec X)
281: {
282: int info;
283: PetscInt i, j, mx, my;
284: PetscInt xs, ys, xm, ym, xe, ye;
285: PetscReal hx, hy, temp, val;
286: double **x;
288: /* Get local mesh boundaries */
289: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
290: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
291: hx = 1.0/(mx-1); hy = 1.0/(my-1);
293: info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
294: xe = xs+xm; ye = ys+ym;
296: info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
297: /* Compute initial guess over locally owned part of mesh */
298: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
299: temp = PetscMin(j+1,my-j)*hy;
300: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
301: val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
302: x[j][i] = val;
303: }
304: }
305: info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);
307: return 0;
308: }
311: /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/
315: /*
316: FormBounds - Forms bounds on the variables
318: Input:
319: user - user-defined application context
321: Output:
322: XL - vector of lower bounds
323: XU - vector of upper bounds
324: */
325: static int DASetBounds(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr)
326: {
327: AppCtx *user = (AppCtx*)ptr;
328: int info;
329: PetscInt i, j, xs, xm, ys, ym;
330: double hx, hy, u1, u2, dist, d1, d2, hd, vd;
331: double **xl, **xu;
333: hx = user->hx;
334: hy = user->hy;
335: u1 = user->u1;
336: u2 = user->u2;
338: info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
339: info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
340: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
342: for (j = ys; j < ys+ym; j++){
343: for (i = xs; i < xs+xm; i++){
344: d1 = i * hx; d2 = u1 - d1; hd = PetscMin(d1,d2);
345: d1 = j * hy; d2 = u2 - d1; vd = PetscMin(d1,d2);
346: dist = PetscMin(hd,vd);
347: xl[j][i] = -dist;
348: xu[j][i] = dist;
349: }
350: }
352: info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
353: info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);
355: info = PetscLogFlops(xm * ym * 4); CHKERRQ(info);
356: return 0;
358: } /* DASetBounds */
363: /*
364: EPTorsLocalFunctionGradient - Evaluates function and gradient over the
365: local rectangular element
367: Input:
368: coor - vector with the indices of the position of current element
369: in the first, second and third directions
370: x - current point (values over the current rectangular element)
371: df - degrees of freedom at each point
372: ptr - user-defined application context
374: Output:
375: f - value of the objective funtion at the local rectangular element
376: g - gradient of the local function
377: */
378: static int EPTorsLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {
380: AppCtx *user = (AppCtx*)ptr;
382: double fquad, flin;
383: double hx, hy, dvdx, dvdy, area;
384: double cdiv3, cnt;
386: cdiv3 = user->param / 3.0;
387: hx = user->hx;
388: hy = user->hy;
389: area = user->area;
390: cnt = area * cdiv3;
392: /* lower triangle contribution */
393: dvdx = (x[0] - x[1]) / hx;
394: dvdy = (x[0] - x[2]) / hy;
395: fquad = dvdx * dvdx + dvdy * dvdy;
396: flin = x[0] + x[1] + x[2];
398: dvdx = 0.5 * dvdx * hy;
399: dvdy = 0.5 * dvdy * hx;
400: g[0] = dvdx + dvdy - cnt;
401: g[1] = -dvdx - 2.0 * cnt;
402: g[2] = -dvdy - 2.0 * cnt;
404: /* upper triangle contribution */
405: dvdx = (x[3] - x[2]) / hx;
406: dvdy = (x[3] - x[1]) / hy;
407: fquad += dvdx * dvdx + dvdy * dvdy;
408: flin += x[1] + x[2] + x[3];
410: dvdx = 0.5 * dvdx * hy;
411: dvdy = 0.5 * dvdy * hx;
412: g[1] += -dvdy;
413: g[2] += -dvdx;
414: g[3] = dvdx + dvdy - cnt;
416: *f = area * (0.5 * fquad - flin * cdiv3);
418: return 0;
419: } /* EPTorsLocalFunctionGradient */
423: /*------- USER-DEFINED: routine to evaluate the Hessian
424: at a local (rectangular element) level -------*/
427: /*
428: EPTorsLocalHessian - Computes the Hessian of the local (partial) function
429: defined over the current rectangle
431: Input:
432: coor - vector with the indices of the position of current element
433: in the first, second and third directions
434: x - current local solution (over the rectangle only)
435: df - degrees of freedom at each point
436: ptr - user-defined application context
438: Output:
439: H - Hessian matrix of the local function (wrt the four
440: points of the rectangle only)
441: */
442: static int EPTorsLocalHessian(PetscInt coor[2], double x[4], double H[4][4], void *ptr) {
444: AppCtx *user = (AppCtx*)ptr;
445: double hx, hy, dxdy, dydx;
446: double diagxy, bandxy, bandyx;
448: hx = user->hx;
449: hy = user->hy;
450: dxdy = hx/hy;
451: dydx = hy/hx;
452: diagxy = 0.5 * (dxdy + dydx);
453: bandxy = -0.5 * dxdy;
454: bandyx = -0.5 * dydx;
456: /* Hessian contribution at 0,0 */
457: H[0][0] = diagxy;
458: H[0][1] = H[1][0] = bandyx;
459: H[0][2] = H[2][0] = bandxy;
460: H[0][3] = H[3][0] = 0.0;
462: /* Hessian contribution at 1,0 */
463: H[1][1] = diagxy;
464: H[1][2] = H[2][1] = 0.0;
465: H[1][3] = H[3][1] = bandxy;
467: /* Hessian contribution at 0,1 */
468: H[2][2] = diagxy;
469: H[2][3] = H[3][2] = bandyx;
471: /* Hessian contribution at 1,1 */
472: H[3][3] = diagxy;
474: return 0;
476: } /* EPTorsLocalHessian */
479: /*------- USER-DEFINED: routine to evaluate the function
480: and gradient at the whole grid -------*/
483: /*
484: WholeEPTorsFunctionGradient - Evaluates function and gradient over the
485: whole grid
487: Input:
488: daapplication - TAO application object
489: da - distributed array
490: X - the current point, at which the function and gradient are evaluated
491: ptr - user-defined application context
493: Output:
494: f - value of the objective funtion at X
495: G - gradient at X
496: */
497: static int WholeEPTorsFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {
499: AppCtx *user = (AppCtx*)ptr;
500: Vec localX, localG;
501: int info;
502: PetscInt i, j;
503: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
504: double **x, **g;
505: double floc = 0.0;
506: PetscScalar zero = 0.0;
508: double fquad, flin;
509: double hx, hy, dvdx, dvdy, area;
510: double cdiv3, cnt;
512: cdiv3 = user->param / 3.0;
513: hx = user->hx;
514: hy = user->hy;
515: area = user->area;
516: cnt = area * cdiv3;
518: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
519: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
520: info = VecSet(G, zero); CHKERRQ(info);
521: info = VecSet(localG, zero); CHKERRQ(info);
523: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
524: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
526: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
527: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
529: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
530: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
532: xe = gxs + gxm - 1;
533: ye = gys + gym - 1;
534: for (j = ys; j < ye; j++) {
535: for (i = xs; i < xe; i++) {
537: /* lower triangle contribution */
538: dvdx = (x[j][i] - x[j][i+1]) / hx;
539: dvdy = (x[j][i] - x[j+1][i]) / hy;
540: fquad = dvdx * dvdx + dvdy * dvdy;
541: flin = x[j][i] + x[j][i+1] + x[j+1][i];
543: dvdx = 0.5 * dvdx * hy;
544: dvdy = 0.5 * dvdy * hx;
545: g[j][i] += dvdx + dvdy - cnt;
546: g[j][i+1] += -dvdx - 2.0 * cnt;
547: g[j+1][i] += -dvdy - 2.0 * cnt;
549: /* upper triangle contribution */
550: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
551: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
552: fquad += dvdx * dvdx + dvdy * dvdy;
553: flin += x[j][i+1] + x[j+1][i] + x[j+1][i+1];
555: dvdx = 0.5 * dvdx * hy;
556: dvdy = 0.5 * dvdy * hx;
557: g[j][i+1] += -dvdy;
558: g[j+1][i] += -dvdx;
559: g[j+1][i+1] += dvdx + dvdy - cnt;
561: floc += area * (0.5 * fquad - flin * cdiv3);
563: }
564: }
566: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);
568: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
569: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
571: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
572: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
574: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
575: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
577: info = PetscLogFlops((xe-xs) * (ye-ys) * 47 + 2); CHKERRQ(info);
578: return 0;
579: } /* WholeEPTorsFunctionGradient */
582: /*------- USER-DEFINED: routine to evaluate the Hessian
583: at the whole grid -------*/
586: /*
587: WholeEPTorsHessian - Evaluates Hessian over the whole grid
589: Input:
590: daapplication - TAO application object
591: da - distributed array
592: X - the current point, at which the function and gradient are evaluated
593: ptr - user-defined application context
595: Output:
596: H - Hessian at X
597: */
598: static int WholeEPTorsHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {
600: AppCtx *user = (AppCtx*)ptr;
601: int info;
602: PetscInt i, j, ind[4];
603: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
604: double smallH[4][4];
606: double hx, hy, dxdy, dydx;
607: double diagxy, bandxy, bandyx;
608: PetscTruth assembled;
610: hx = user->hx;
611: hy = user->hy;
612: dxdy = hx/hy;
613: dydx = hy/hx;
614: diagxy = 0.5 * (dxdy + dydx);
615: bandxy = -0.5 * dxdy;
616: bandyx = -0.5 * dydx;
618: info = MatAssembled(H,&assembled); CHKERRQ(info);
619: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
622: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
623: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
625: xe = gxs + gxm - 1;
626: ye = gys + gym - 1;
627: for (j = ys; j < ye; j++) {
628: for (i = xs; i < xe; i++) {
630: /* Hessian contribution at 0,0 */
631: smallH[0][0] = diagxy;
632: smallH[0][1] = smallH[1][0] = bandyx;
633: smallH[0][2] = smallH[2][0] = bandxy;
634: smallH[0][3] = smallH[3][0] = 0.0;
636: /* Hessian contribution at 1,0 */
637: smallH[1][1] = diagxy;
638: smallH[1][2] = smallH[2][1] = 0.0;
639: smallH[1][3] = smallH[3][1] = bandxy;
641: /* Hessian contribution at 0,1 */
642: smallH[2][2] = diagxy;
643: smallH[2][3] = smallH[3][2] = bandyx;
645: /* Hessian contribution at 1,1 */
646: smallH[3][3] = diagxy;
648: ind[0] = (j-gys) * gxm + (i-gxs);
649: ind[1] = ind[0] + 1;
650: ind[2] = ind[0] + gxm;
651: ind[3] = ind[2] + 1;
652: info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);
654: }
655: }
657: info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
658: info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
659: info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);
662: info = PetscLogFlops(6); CHKERRQ(info);
663: return 0;
665: } /* WholeEPTorsHessian */