Actual source code: combustion3.c
1: /*$Id: combustion.c,v 1.1 2002/08/08 9:39 lopezca@mauddib.mcs.anl.gov $*/
3: /* Program usage: mpirun -np <proc> combustion [-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[] = "Steady-State Combustion.\n\
16: We solve the Steady-State Combustion problem (MINPACK-2 test suite) in a 2D\n\
17: rectangular domain, using distributed arrays (DAs) to partition the parallel grid.\n\
18: The command line options include:\n\
19: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
20: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
21: -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
22: -byelement, if computation is made by functions on rectangular elements\n\
23: -adic, if AD is used (AD is not used by default)\n\
24: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
25: parameter lambda (0 <= par <= 6.81)\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: TaoAppSetOptions();
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: */
51: typedef struct {
52: InactiveDouble param;
53: InactiveDouble hx, hy; /* increment size in both directions */
54: InactiveDouble area; /* area of the triangles */
55: } ADFGCtx;
57: typedef struct {
58: PetscReal param; /* nonlinearity parameter */
59: double hx, hy, area; /* increments and area of the triangle */
60: int mx, my; /* discretization including boundaries */
61: ADFGCtx fgctx; /* Used only when an ADIC generated gradient is used */
62: } AppCtx;
63: int ad_CombLocalFunction(int[2], DERIV_TYPE[4], DERIV_TYPE*, void*);
65: /* User-defined routines foun in this file */
66: static int AppCtxInitialize(void *ptr);
67: static int FormInitialGuess(DA, Vec, AppCtx*);
69: static int CombLocalFunctionGradient(int[3], double x[4], double *f, double g[4], void *ptr);
70: static int WholeCombFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);
72: static int CombLocalHessian(int[3], double x[4], double H[4][4], void *ptr);
73: static int WholeCombHessian(TAO_APPLICATION,DA,Vec,Mat,void*);
75: static int DAFixBoundary(TAO_APPLICATION, DA, Vec, Vec, void*);
77: static int MyGridMonitorBefore(TAO_APPLICATION, DA, int, void *);
82: int main( int argc, char **argv ) {
84: int info,iter; /* used to check for functions returning nonzeros */
85: int Nx,Ny;
86: int nlevels; /* multigrid levels */
87: double ff,gnorm;
88: DA DAarray[20];
89: Vec X;
90: PetscTruth flg, PreLoad = PETSC_TRUE; /* flags */
91: AppCtx user; /* user-defined work context */
92: TaoMethod method = "tao_tron"; /* minimization method */
93: TAO_SOLVER tao; /* TAO_SOLVER solver context */
94: TAO_APPLICATION CombApp; /* The PETSc application */
95: TaoTerminateReason reason;
97: /* Initialize TAO */
98: PetscInitialize(&argc, &argv, (char *)0, help);
99: TaoInitialize(&argc, &argv, (char *)0, help);
101: PreLoadBegin(PreLoad,"Solve");
102:
103: info = AppCtxInitialize((void*)&user); CHKERRQ(info);
105: nlevels=5;
106: info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
107: if (PreLoadIt == 0) {
108: nlevels = 1; user.mx = 11; user.my = 11;}
110: PetscPrintf(MPI_COMM_WORLD,"\n---- Steady-State Combustion Problem -----\n\n");
112: /* Let PETSc determine the vector distribution */
113: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
115: /* Create distributed array (DA) to manage parallel grid and vectors */
116: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
117: user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
118: for (iter=1;iter<nlevels;iter++){
119: info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
120: }
122: /* Create TAO solver and set desired solution method */
123: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
125: info = TaoApplicationCreate(PETSC_COMM_WORLD, &CombApp); CHKERRQ(info);
126: info = TaoAppSetDAApp(CombApp,DAarray,nlevels); CHKERRQ(info);
128: /* Sets routines for function, gradient and bounds evaluation */
129: info = DAAppSetVariableBoundsRoutine(CombApp,DAFixBoundary,(void *)&user); CHKERRQ(info);
131: info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
132: if (flg) {
134: /* Sets routines for function and gradient evaluation, element by element */
135: info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
136: if (flg) {
137: info = DAAppSetADElementFunctionGradient(CombApp,ad_CombLocalFunction,228,(void *)&user.fgctx); CHKERRQ(info);
138: } else {
139: info = DAAppSetElementObjectiveAndGradientRoutine(CombApp,CombLocalFunctionGradient,51,(void *)&user); CHKERRQ(info);
140: }
141: /* Sets routines for Hessian evaluation, element by element */
142: info = DAAppSetElementHessianRoutine(CombApp,CombLocalHessian,21,(void*)&user); CHKERRQ(info);
144: } else {
146: /* Sets routines for function and gradient evaluation, all in one routine */
147: info = DAAppSetObjectiveAndGradientRoutine(CombApp,WholeCombFunctionGradient,(void *)&user); CHKERRQ(info);
149: /* Sets routines for Hessian evaluation, all in one routine */
150: info = DAAppSetHessianRoutine(CombApp,WholeCombHessian,(void*)&user); CHKERRQ(info);
151:
152: }
154: info = DAAppSetBeforeMonitor(CombApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
155: info = DAAppPrintStageTimes(CombApp); CHKERRQ(info);
156: info = DAAppPrintInterpolationError(CombApp); CHKERRQ(info);
158: info = TaoAppSetRelativeTolerance(CombApp,1.0e-8); CHKERRQ(info);
159: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
160: info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);
162: /* Check for any tao command line options */
163: info = TaoSetOptions(CombApp, tao); CHKERRQ(info);
165: info = DAAppGetSolution(CombApp,0,&X); CHKERRQ(info);
166: info = FormInitialGuess(DAarray[0],X,&user); CHKERRQ(info);
167: info = DAAppSetInitialSolution(CombApp,X); CHKERRQ(info);
169: /* SOLVE THE APPLICATION */
170: info = TaoDAAppSolve(CombApp, tao); CHKERRQ(info);
172: /* Get information on termination */
173: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
174: if (reason <= 0 ){
175: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
176: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
177: }
179: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
180: if (flg){
181: info = DAAppGetSolution(CombApp,nlevels-1,&X); CHKERRQ(info);
182: info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
183: }
185: /* To View TAO solver information */
186: // info = TaoView(tao); CHKERRQ(info);
188: /* Free TAO data structures */
189: info = TaoDestroy(tao); CHKERRQ(info);
190: info = TaoAppDestroy(CombApp); CHKERRQ(info);
192: /* Free PETSc data structures */
193: for (iter=0;iter<nlevels;iter++){
194: info = DADestroy(DAarray[iter]); CHKERRQ(info);
195: }
197: PreLoadEnd();
199: /* Finalize TAO */
200: TaoFinalize();
201: PetscFinalize();
203: return 0;
204: } /* main */
208: /*----- The following two routines
209:
210: MyGridMonitorBefore MyGridMonitorAfter
212: help diplay info of iterations at every grid level -------*/
216: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, int level, void *ctx) {
218: AppCtx *user = (AppCtx*)ctx;
219: int info,mx,my;
221: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
222: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
223: user->mx = mx;
224: user->my = my;
225: user->hx = 1.0 / (user->mx - 1);
226: user->hy = 1.0 / (user->my - 1);
227: user->area = 0.5 * user->hx * user->hy;
228: user->fgctx.hx = user->hx;
229: user->fgctx.hy = user->hy;
230: user->fgctx.area = user->area;
231: user->fgctx.param = user->param;
233: PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my);
235: return 0;
236: }
238: /*------- USER-DEFINED: initialize the application context information -------*/
242: /*
243: AppCtxInitialize - Sets initial values for the application context parameters
245: Input:
246: ptr - void user-defined application context
248: Output:
249: ptr - user-defined application context with the default or user-provided
250: parameters
251: */
252: static int AppCtxInitialize(void *ptr) {
254: AppCtx *user = (AppCtx*)ptr;
255: PetscReal LambdaMax = 6.81, LambdaMin = 0.0; /* bounds on parameter lambda */
256: PetscTruth flg; /* flag for PETSc calls */
257: int info;
259: /* Specify dimension of the problem */
260: user->param = 5.0;
261: user->mx = 11;
262: user->my = 11;
264: /* Check for any command line arguments that override defaults */
265: info = PetscOptionsGetReal(TAO_NULL, "-par", &user->param, &flg); CHKERRQ(info);
266: if (user->param >= LambdaMax || user->param <= LambdaMin) {
267: SETERRQ(1,"Lambda is out of range.");
268: }
269: info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user->mx,&flg); CHKERRQ(info);
270: info = PetscOptionsGetInt(PETSC_NULL,"-my",&user->my,&flg); CHKERRQ(info);
272: user->hx = 1.0 / (user->mx - 1);
273: user->hy = 1.0 / (user->my - 1);
274: user->area = 0.5 * user->hx * user->hy;
275: info = PetscLogFlops(6); CHKERRQ(info);
277: return 0;
278: } /* AppCtxInitialize */
284: static int FormInitialGuess(DA da, Vec X, AppCtx *ctx)
285: {
286: int info, i, j, mx, my;
287: int xs, ys, xm, ym, xe, ye;
288: PetscReal hx, hy, temp, val, lambda;
289: double **x;
291: lambda = ctx->param;
292: lambda = lambda/(lambda+1.0);
294: /* Get local mesh boundaries */
295: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
296: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
297: hx = 1.0/(mx-1); hy = 1.0/(my-1);
299: info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
300: xe = xs+xm; ye = ys+ym;
302: info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
303: /* Compute initial guess over locally owned part of mesh */
304: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
305: temp = PetscMin(j+1,my-j)*hy;
306: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
307: val = lambda*sqrt(PetscMin((PetscMin(i+1,mx-i))*hx,temp));
308: x[j][i] = val;
309: }
310: }
311: info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);
313: return 0;
314: }
317: /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/
321: /*
322: FormBounds - Forms bounds on the variables
324: Input:
325: user - user-defined application context
327: Output:
328: XL - vector of lower bounds
329: XU - vector of upper bounds
331: */
332: static int DAFixBoundary(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr)
333: {
334: AppCtx *user = (AppCtx*)ptr;
335: int i, j, mx, my, info, xs, xm, ys, ym;
336: double lb = -TAO_INFINITY;
337: double ub = TAO_INFINITY;
338: double **xl, **xu;
340: mx = user->mx; /* number of points including the boundary */
341: my = user->my;
343: info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
344: info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
345: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
347: /* Compute initial guess over locally owned part of the grid */
348: for (j = ys; j < ys+ym; j++){
349: for (i = xs; i < xs+xm; i++){
350: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
351: xl[j][i] = xu[j][i] = 0.0;
352: } else {
353: xl[j][i] = lb;
354: xu[j][i] = ub;
355: }
356: }
357: }
359: info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
360: info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);
362: return 0;
363: } /* DAFixBoundary */
366: /*------- USER-DEFINED: routine to evaluate the function and gradient
367: at a local (rectangular element) level -------*/
371: /*
372: CombLocalFunctionGradient - Evaluates function and gradient over the
373: local rectangular element
375: Input:
376: coor - vector with the indices of the position of current element
377: in the first, second and third directions
378: x - current point (values over the current rectangular element)
379: ptr - user-defined application context
381: Output:
382: f - value of the objective funtion at the local rectangular element
383: g - gradient of the local function
385: */
386: static int CombLocalFunctionGradient(int coor[2], double x[4], double *f, double g[4], void *ptr) {
388: AppCtx *user = (AppCtx*)ptr;
390: double lambdad3, hx, hy, area;
391: double fquad, fexp, dvdx, dvdy;
393: lambdad3 = user->param / 3.0;
394: hx = user->hx;
395: hy = user->hy;
396: area = user->area;
398: /* lower triangle contribution */
399: dvdx = (x[0] - x[1]) / hx;
400: dvdy = (x[0] - x[2]) / hy;
401: fquad = dvdx * dvdx + dvdy * dvdy;
402: fexp = exp(x[0]) + exp(x[1]) + exp(x[2]);
404: dvdx = 0.5 * dvdx * hy;
405: dvdy = 0.5 * dvdy * hx;
406: g[0] = dvdx + dvdy - exp(x[0]) * area * lambdad3;
407: g[1] = -dvdx - 2.0 * exp(x[1]) * area * lambdad3;
408: g[2] = -dvdy - 2.0 * exp(x[2]) * area * lambdad3;
410: /* upper triangle contribution */
411: dvdx = (x[3] - x[2]) / hx;
412: dvdy = (x[3] - x[1]) / hy;
413: fquad += dvdx * dvdx + dvdy * dvdy;
414: fexp += exp(x[1]) + exp(x[2]) + exp(x[3]);
416: dvdx = 0.5 * dvdx * hy;
417: dvdy = 0.5 * dvdy * hx;
418: g[1] += -dvdy;
419: g[2] += -dvdx;
420: g[3] = dvdx + dvdy - exp(x[3]) * area * lambdad3;
423: *f = area * (0.5 * fquad - lambdad3 * fexp);
425: return 0;
426: } /* CombLocalFunctionGradient */
429: /*------- USER-DEFINED: routine to evaluate the Hessian
430: at a local (rectangular element) level -------*/
434: /*
435: CombLocalHessian - Computes the Hessian of the local (partial) function
436: defined over the current rectangle
438: Input:
439: coor - vector with the indices of the position of current element
440: in the first, second and third directions
441: x - current local solution (over the rectangle only)
442: ptr - user-defined application context
444: Output:
445: H - Hessian matrix of the local function (wrt the four
446: points of the rectangle only)
448: */
449: static int CombLocalHessian(int coor[2], double x[4], double H[4][4],void *ptr) {
451: AppCtx *user = (AppCtx*)ptr;
452: double hx, hy, lambdad3, area, dxdy, dydx;
453: double diagxy, dexp, bandxy, bandyx;
455: hx = user->hx;
456: hy = user->hy;
457: lambdad3 = user->param / 3.0;
458: area = user->area;
459: dxdy = hx/hy;
460: dydx = hy/hx;
461: diagxy = 0.5 * (dxdy + dydx);
462: bandxy = -0.5 * dxdy;
463: bandyx = -0.5 * dydx;
465: /* Hessian contribution at 0,0 */
466: dexp = exp(x[0]) * area * lambdad3;
467: H[0][0] = diagxy - dexp;
468: H[0][1] = H[1][0] = bandyx;
469: H[0][2] = H[2][0] = bandxy;
470: H[0][3] = H[3][0] = 0.0;
472: /* Hessian contribution at 1,0 */
473: dexp = exp(x[1]) * area * 2.0 * lambdad3;
474: H[1][1] = diagxy - dexp;
475: H[1][2] = H[2][1] = 0.0;
476: H[1][3] = H[3][1] = bandxy;
478: /* Hessian contribution at 0,1 */
479: dexp = exp(x[2]) * area * 2.0 * lambdad3;
480: H[2][2] = diagxy - dexp;
481: H[2][3] = H[3][2] = bandyx;
483: /* Hessian contribution at 1,1 */
484: dexp = exp(x[3]) * area * lambdad3;
485: H[3][3] = diagxy - dexp;
487: return 0;
488: } /* CombLocalHessian */
491: /*------- USER-DEFINED: routine to evaluate the function
492: and gradient at the whole grid -------*/
496: /*
497: WholeCombFunctionGradient - Evaluates function and gradient over the
498: whole grid
500: Input:
501: daapplication - TAO application object
502: da - distributed array
503: X - the current point, at which the function and gradient are evaluated
504: ptr - user-defined application context
506: Output:
507: f - value of the objective funtion at X
508: G - gradient at X
509: */
510: static int WholeCombFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {
512: AppCtx *user = (AppCtx*)ptr;
513: Vec localX, localG;
514: int info, i, j;
515: int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
516: double **x, **g;
517: double floc = 0.0;
518: PetscScalar zero = 0.0;
520: double lambdad3, hx, hy, area;
521: double fquad, fexp, dvdx, dvdy;
523: lambdad3 = user->param / 3.0;
524: hx = user->hx;
525: hy = user->hy;
526: area = user->area;
528: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
529: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
530: info = VecSet(G, zero); CHKERRQ(info);
531: info = VecSet(localG, zero); CHKERRQ(info);
533: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
534: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
536: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
537: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
539: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
540: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
542: xe = gxs + gxm - 1;
543: ye = gys + gym - 1;
544: for (j = ys; j < ye; j++) {
545: for (i = xs; i < xe; i++) {
547: /* lower triangle contribution */
548: dvdx = (x[j][i] - x[j][i+1]) / hx;
549: dvdy = (x[j][i] - x[j+1][i]) / hy;
550: fquad = dvdx * dvdx + dvdy * dvdy;
551: fexp = exp(x[j][i]) + exp(x[j][i+1]) + exp(x[j+1][i]);
553: dvdx = 0.5 * dvdx * hy;
554: dvdy = 0.5 * dvdy * hx;
555: g[j][i] += dvdx + dvdy - exp(x[j][i]) * area * lambdad3;
556: g[j][i+1] += -dvdx - 2.0 * exp(x[j][i+1]) * area * lambdad3;
557: g[j+1][i] += -dvdy - 2.0 * exp(x[j+1][i]) * area * lambdad3;
559: /* upper triangle contribution */
560: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
561: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
562: fquad += dvdx * dvdx + dvdy * dvdy;
563: fexp += exp(x[j][i+1]) + exp(x[j+1][i]) + exp(x[j+1][i+1]);
565: dvdx = 0.5 * dvdx * hy;
566: dvdy = 0.5 * dvdy * hx;
567: g[j][i+1] += -dvdy;
568: g[j+1][i] += -dvdx;
569: g[j+1][i+1] += dvdx + dvdy - exp(x[j+1][i+1]) * area * lambdad3;
571: floc += area * (0.5 * fquad - fexp * lambdad3);
573: }
574: }
576: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);
578: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
579: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
581: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
582: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
584: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
585: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
587: info = PetscLogFlops((xe-xs) * (ye-ys) * 55 + 1); CHKERRQ(info);
588: return 0;
590: } /* WholeCombFunctionGradient */
593: /*------- USER-DEFINED: routine to evaluate the Hessian
594: at the whole grid -------*/
597: /*
598: WholeCombHessian - Evaluates Hessian over the whole grid
600: Input:
601: daapplication - TAO application object
602: da - distributed array
603: X - the current point, at which the function and gradient are evaluated
604: ptr - user-defined application context
606: Output:
607: H - Hessian at X
608: */
609: static int WholeCombHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {
611: AppCtx *user = (AppCtx*)ptr;
612: Vec localX;
613: int info, i, j, ind[4];
614: int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
615: double smallH[4][4];
616: double **x;
618: double hx, hy, lambdad3, area, dxdy, dydx;
619: double diagxy, dexp, bandxy, bandyx;
620: PetscTruth assembled;
623: hx = user->hx;
624: hy = user->hy;
625: lambdad3 = user->param / 3.0;
626: area = user->area;
627: dxdy = hx/hy;
628: dydx = hy/hx;
629: diagxy = 0.5 * (dxdy + dydx);
630: bandxy = -0.5 * dxdy;
631: bandyx = -0.5 * dydx;
633: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
634: info = MatAssembled(H,&assembled); CHKERRQ(info);
635: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
637: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
638: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
640: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
642: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
643: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
645: xe = gxs + gxm - 1;
646: ye = gys + gym - 1;
647: for (j = ys; j < ye; j++) {
648: for (i = xs; i < xe; i++) {
650: /* Hessian contribution at 0,0 */
651: dexp = exp(x[j][i]) * area * lambdad3;
652: smallH[0][0] = diagxy - dexp;
653: smallH[0][1] = smallH[1][0] = bandyx;
654: smallH[0][2] = smallH[2][0] = bandxy;
655: smallH[0][3] = smallH[3][0] = 0.0;
657: /* Hessian contribution at 1,0 */
658: dexp = exp(x[j][i+1]) * area * 2.0 * lambdad3;
659: smallH[1][1] = diagxy - dexp;
660: smallH[1][2] = smallH[2][1] = 0.0;
661: smallH[1][3] = smallH[3][1] = bandxy;
663: /* Hessian contribution at 0,1 */
664: dexp = exp(x[j+1][i]) * area * 2.0 * lambdad3;
665: smallH[2][2] = diagxy - dexp;
666: smallH[2][3] = smallH[3][2] = bandyx;
668: /* Hessian contribution at 1,1 */
669: dexp = exp(x[j+1][i+1]) * area * lambdad3;
670: smallH[3][3] = diagxy - dexp;
672: ind[0] = (j-gys) * gxm + (i-gxs);
673: ind[1] = ind[0] + 1;
674: ind[2] = ind[0] + gxm;
675: ind[3] = ind[2] + 1;
676: info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);
678: }
679: }
681: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
683: info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
684: info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
685: info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);
687: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
689: info = PetscLogFlops((xe-xs) * (ye-ys) * 14 + 7); CHKERRQ(info);
690: return 0;
692: } /* WholeCombHessian */