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: PetscInt mx, my; /* discretization including boundaries */
61: ADFGCtx fgctx; /* Used only when an ADIC generated gradient is used */
62: } AppCtx;
63: int ad_CombLocalFunction(PetscInt[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(PetscInt[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(PetscInt[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, PetscInt, void *);
82: int main( int argc, char **argv ) {
84: int info; /* used to check for functions returning nonzeros */
85: PetscInt Nx,Ny,iter;
86: PetscInt nlevels; /* multigrid levels */
87: double ff,gnorm;
88: DA DAarray[20];
89: Vec X;
90: KSP ksp;
91: PetscTruth flg, PreLoad = PETSC_TRUE; /* flags */
92: AppCtx user; /* user-defined work context */
93: TaoMethod method = "tao_tron"; /* minimization method */
94: TAO_SOLVER tao; /* TAO_SOLVER solver context */
95: TAO_APPLICATION CombApp; /* The PETSc application */
96: TaoTerminateReason reason;
98: /* Initialize TAO */
99: PetscInitialize(&argc, &argv, (char *)0, help);
100: TaoInitialize(&argc, &argv, (char *)0, help);
102: PreLoadBegin(PreLoad,"Solve");
103:
104: info = AppCtxInitialize((void*)&user); CHKERRQ(info);
106: nlevels=5;
107: info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
108: if (PreLoadIt == 0) {
109: nlevels = 1; user.mx = 11; user.my = 11;}
111: PetscPrintf(MPI_COMM_WORLD,"\n---- Steady-State Combustion Problem -----\n\n");
113: /* Let PETSc determine the vector distribution */
114: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
116: /* Create distributed array (DA) to manage parallel grid and vectors */
117: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
118: user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
119: for (iter=1;iter<nlevels;iter++){
120: info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
121: }
123: /* Create TAO solver and set desired solution method */
124: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
126: info = TaoApplicationCreate(PETSC_COMM_WORLD, &CombApp); CHKERRQ(info);
127: info = TaoAppSetDAApp(CombApp,DAarray,nlevels); CHKERRQ(info);
129: /* Sets routines for function, gradient and bounds evaluation */
130: info = DAAppSetVariableBoundsRoutine(CombApp,DAFixBoundary,(void *)&user); CHKERRQ(info);
132: info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
133: if (flg) {
135: /* Sets routines for function and gradient evaluation, element by element */
136: info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
137: if (flg) {
138: info = DAAppSetADElementFunctionGradient(CombApp,ad_CombLocalFunction,228,(void *)&user.fgctx); CHKERRQ(info);
139: } else {
140: info = DAAppSetElementObjectiveAndGradientRoutine(CombApp,CombLocalFunctionGradient,51,(void *)&user); CHKERRQ(info);
141: }
142: /* Sets routines for Hessian evaluation, element by element */
143: info = DAAppSetElementHessianRoutine(CombApp,CombLocalHessian,21,(void*)&user); CHKERRQ(info);
145: } else {
147: /* Sets routines for function and gradient evaluation, all in one routine */
148: info = DAAppSetObjectiveAndGradientRoutine(CombApp,WholeCombFunctionGradient,(void *)&user); CHKERRQ(info);
150: /* Sets routines for Hessian evaluation, all in one routine */
151: info = DAAppSetHessianRoutine(CombApp,WholeCombHessian,(void*)&user); CHKERRQ(info);
152:
153: }
155: info = DAAppSetBeforeMonitor(CombApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
156: info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
157: if (flg){
158: info = DAAppPrintInterpolationError(CombApp); CHKERRQ(info);
159: info = DAAppPrintStageTimes(CombApp); CHKERRQ(info);
160: }
163: info = TaoAppSetRelativeTolerance(CombApp,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(CombApp,&ksp);CHKERRQ(info);
168: info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
169: info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100);CHKERRQ(info);
171: /* Check for any tao command line options */
172: info = TaoSetOptions(CombApp, tao); CHKERRQ(info);
174: info = DAAppGetSolution(CombApp,0,&X); CHKERRQ(info);
175: info = FormInitialGuess(DAarray[0],X,&user); CHKERRQ(info);
176: info = DAAppSetInitialSolution(CombApp,X); CHKERRQ(info);
178: /* SOLVE THE APPLICATION */
179: info = TaoDAAppSolve(CombApp, tao); CHKERRQ(info);
181: /* Get information on termination */
182: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
183: if (reason <= 0 ){
184: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
185: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
186: }
188: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
189: if (flg){
190: info = DAAppGetSolution(CombApp,nlevels-1,&X); CHKERRQ(info);
191: info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
192: }
194: /* To View TAO solver information */
195: // info = TaoView(tao); CHKERRQ(info);
197: /* Free TAO data structures */
198: info = TaoDestroy(tao); CHKERRQ(info);
199: info = TaoAppDestroy(CombApp); CHKERRQ(info);
201: /* Free PETSc data structures */
202: for (iter=0;iter<nlevels;iter++){
203: info = DADestroy(DAarray[iter]); CHKERRQ(info);
204: }
206: PreLoadEnd();
208: /* Finalize TAO */
209: TaoFinalize();
210: PetscFinalize();
212: return 0;
213: } /* main */
217: /*----- The following two routines
218:
219: MyGridMonitorBefore MyGridMonitorAfter
221: help diplay info of iterations at every grid level -------*/
225: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx) {
227: AppCtx *user = (AppCtx*)ctx;
228: int info;
229: PetscInt mx,my;
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->mx = mx;
234: user->my = my;
235: user->hx = 1.0 / (user->mx - 1);
236: user->hy = 1.0 / (user->my - 1);
237: user->area = 0.5 * user->hx * user->hy;
238: user->fgctx.hx = user->hx;
239: user->fgctx.hy = user->hy;
240: user->fgctx.area = user->area;
241: user->fgctx.param = user->param;
243: PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my);
244: return 0;
245: }
247: /*------- USER-DEFINED: initialize the application context information -------*/
251: /*
252: AppCtxInitialize - Sets initial values for the application context parameters
254: Input:
255: ptr - void user-defined application context
257: Output:
258: ptr - user-defined application context with the default or user-provided
259: parameters
260: */
261: static int AppCtxInitialize(void *ptr) {
263: AppCtx *user = (AppCtx*)ptr;
264: PetscReal LambdaMax = 6.81, LambdaMin = 0.0; /* bounds on parameter lambda */
265: PetscTruth flg; /* flag for PETSc calls */
266: int info;
268: /* Specify dimension of the problem */
269: user->param = 5.0;
270: user->mx = 11;
271: user->my = 11;
273: /* Check for any command line arguments that override defaults */
274: info = PetscOptionsGetReal(TAO_NULL, "-par", &user->param, &flg); CHKERRQ(info);
275: if (user->param >= LambdaMax || user->param <= LambdaMin) {
276: SETERRQ(1,"Lambda is out of range.");
277: }
278: info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user->mx,&flg); CHKERRQ(info);
279: info = PetscOptionsGetInt(PETSC_NULL,"-my",&user->my,&flg); CHKERRQ(info);
281: user->hx = 1.0 / (user->mx - 1);
282: user->hy = 1.0 / (user->my - 1);
283: user->area = 0.5 * user->hx * user->hy;
284: info = PetscLogFlops(6); CHKERRQ(info);
286: return 0;
287: } /* AppCtxInitialize */
293: static int FormInitialGuess(DA da, Vec X, AppCtx *ctx)
294: {
295: int info;
296: PetscInt i, j, mx, my;
297: PetscInt xs, ys, xm, ym, xe, ye;
298: PetscReal hx, hy, temp, val, lambda;
299: double **x;
301: lambda = ctx->param;
302: lambda = lambda/(lambda+1.0);
304: /* Get local mesh boundaries */
305: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
306: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
307: hx = 1.0/(mx-1); hy = 1.0/(my-1);
309: info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
310: xe = xs+xm; ye = ys+ym;
312: info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
313: /* Compute initial guess over locally owned part of mesh */
314: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
315: temp = PetscMin(j+1,my-j)*hy;
316: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
317: val = lambda*sqrt(PetscMin((PetscMin(i+1,mx-i))*hx,temp));
318: x[j][i] = val;
319: }
320: }
321: info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);
323: return 0;
324: }
327: /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/
331: /*
332: FormBounds - Forms bounds on the variables
334: Input:
335: user - user-defined application context
337: Output:
338: XL - vector of lower bounds
339: XU - vector of upper bounds
341: */
342: static int DAFixBoundary(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr)
343: {
344: AppCtx *user = (AppCtx*)ptr;
345: int info;
346: PetscInt i, j, mx, my, xs, xm, ys, ym;
347: double lb = -TAO_INFINITY;
348: double ub = TAO_INFINITY;
349: double **xl, **xu;
351: mx = user->mx; /* number of points including the boundary */
352: my = user->my;
354: info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
355: info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
356: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
358: /* Compute initial guess over locally owned part of the grid */
359: for (j = ys; j < ys+ym; j++){
360: for (i = xs; i < xs+xm; i++){
361: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
362: xl[j][i] = xu[j][i] = 0.0;
363: } else {
364: xl[j][i] = lb;
365: xu[j][i] = ub;
366: }
367: }
368: }
370: info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
371: info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);
373: return 0;
374: } /* DAFixBoundary */
377: /*------- USER-DEFINED: routine to evaluate the function and gradient
378: at a local (rectangular element) level -------*/
382: /*
383: CombLocalFunctionGradient - Evaluates function and gradient over the
384: local rectangular element
386: Input:
387: coor - vector with the indices of the position of current element
388: in the first, second and third directions
389: x - current point (values over the current rectangular element)
390: ptr - user-defined application context
392: Output:
393: f - value of the objective funtion at the local rectangular element
394: g - gradient of the local function
396: */
397: static int CombLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {
399: AppCtx *user = (AppCtx*)ptr;
401: double lambdad3, hx, hy, area;
402: double fquad, fexp, dvdx, dvdy;
404: lambdad3 = user->param / 3.0;
405: hx = user->hx;
406: hy = user->hy;
407: area = user->area;
409: /* lower triangle contribution */
410: dvdx = (x[0] - x[1]) / hx;
411: dvdy = (x[0] - x[2]) / hy;
412: fquad = dvdx * dvdx + dvdy * dvdy;
413: fexp = exp(x[0]) + exp(x[1]) + exp(x[2]);
415: dvdx = 0.5 * dvdx * hy;
416: dvdy = 0.5 * dvdy * hx;
417: g[0] = dvdx + dvdy - exp(x[0]) * area * lambdad3;
418: g[1] = -dvdx - 2.0 * exp(x[1]) * area * lambdad3;
419: g[2] = -dvdy - 2.0 * exp(x[2]) * area * lambdad3;
421: /* upper triangle contribution */
422: dvdx = (x[3] - x[2]) / hx;
423: dvdy = (x[3] - x[1]) / hy;
424: fquad += dvdx * dvdx + dvdy * dvdy;
425: fexp += exp(x[1]) + exp(x[2]) + exp(x[3]);
427: dvdx = 0.5 * dvdx * hy;
428: dvdy = 0.5 * dvdy * hx;
429: g[1] += -dvdy;
430: g[2] += -dvdx;
431: g[3] = dvdx + dvdy - exp(x[3]) * area * lambdad3;
434: *f = area * (0.5 * fquad - lambdad3 * fexp);
436: return 0;
437: } /* CombLocalFunctionGradient */
440: /*------- USER-DEFINED: routine to evaluate the Hessian
441: at a local (rectangular element) level -------*/
445: /*
446: CombLocalHessian - Computes the Hessian of the local (partial) function
447: defined over the current rectangle
449: Input:
450: coor - vector with the indices of the position of current element
451: in the first, second and third directions
452: x - current local solution (over the rectangle only)
453: ptr - user-defined application context
455: Output:
456: H - Hessian matrix of the local function (wrt the four
457: points of the rectangle only)
459: */
460: static int CombLocalHessian(PetscInt coor[2], double x[4], double H[4][4],void *ptr) {
462: AppCtx *user = (AppCtx*)ptr;
463: double hx, hy, lambdad3, area, dxdy, dydx;
464: double diagxy, dexp, bandxy, bandyx;
466: hx = user->hx;
467: hy = user->hy;
468: lambdad3 = user->param / 3.0;
469: area = user->area;
470: dxdy = hx/hy;
471: dydx = hy/hx;
472: diagxy = 0.5 * (dxdy + dydx);
473: bandxy = -0.5 * dxdy;
474: bandyx = -0.5 * dydx;
476: /* Hessian contribution at 0,0 */
477: dexp = exp(x[0]) * area * lambdad3;
478: H[0][0] = diagxy - dexp;
479: H[0][1] = H[1][0] = bandyx;
480: H[0][2] = H[2][0] = bandxy;
481: H[0][3] = H[3][0] = 0.0;
483: /* Hessian contribution at 1,0 */
484: dexp = exp(x[1]) * area * 2.0 * lambdad3;
485: H[1][1] = diagxy - dexp;
486: H[1][2] = H[2][1] = 0.0;
487: H[1][3] = H[3][1] = bandxy;
489: /* Hessian contribution at 0,1 */
490: dexp = exp(x[2]) * area * 2.0 * lambdad3;
491: H[2][2] = diagxy - dexp;
492: H[2][3] = H[3][2] = bandyx;
494: /* Hessian contribution at 1,1 */
495: dexp = exp(x[3]) * area * lambdad3;
496: H[3][3] = diagxy - dexp;
498: return 0;
499: } /* CombLocalHessian */
502: /*------- USER-DEFINED: routine to evaluate the function
503: and gradient at the whole grid -------*/
507: /*
508: WholeCombFunctionGradient - Evaluates function and gradient over the
509: whole grid
511: Input:
512: daapplication - TAO application object
513: da - distributed array
514: X - the current point, at which the function and gradient are evaluated
515: ptr - user-defined application context
517: Output:
518: f - value of the objective funtion at X
519: G - gradient at X
520: */
521: static int WholeCombFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {
523: AppCtx *user = (AppCtx*)ptr;
524: Vec localX, localG;
525: int info;
526: PetscInt i, j;
527: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
528: double **x, **g;
529: double floc = 0.0;
530: PetscScalar zero = 0.0;
532: double lambdad3, hx, hy, area;
533: double fquad, fexp, dvdx, dvdy;
535: lambdad3 = user->param / 3.0;
536: hx = user->hx;
537: hy = user->hy;
538: area = user->area;
540: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
541: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
542: info = VecSet(G, zero); CHKERRQ(info);
543: info = VecSet(localG, zero); CHKERRQ(info);
545: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
546: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
548: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
549: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
551: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
552: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
554: xe = gxs + gxm - 1;
555: ye = gys + gym - 1;
556: for (j = ys; j < ye; j++) {
557: for (i = xs; i < xe; i++) {
559: /* lower triangle contribution */
560: dvdx = (x[j][i] - x[j][i+1]) / hx;
561: dvdy = (x[j][i] - x[j+1][i]) / hy;
562: fquad = dvdx * dvdx + dvdy * dvdy;
563: fexp = exp(x[j][i]) + exp(x[j][i+1]) + exp(x[j+1][i]);
565: dvdx = 0.5 * dvdx * hy;
566: dvdy = 0.5 * dvdy * hx;
567: g[j][i] += dvdx + dvdy - exp(x[j][i]) * area * lambdad3;
568: g[j][i+1] += -dvdx - 2.0 * exp(x[j][i+1]) * area * lambdad3;
569: g[j+1][i] += -dvdy - 2.0 * exp(x[j+1][i]) * area * lambdad3;
571: /* upper triangle contribution */
572: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
573: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
574: fquad += dvdx * dvdx + dvdy * dvdy;
575: fexp += exp(x[j][i+1]) + exp(x[j+1][i]) + exp(x[j+1][i+1]);
577: dvdx = 0.5 * dvdx * hy;
578: dvdy = 0.5 * dvdy * hx;
579: g[j][i+1] += -dvdy;
580: g[j+1][i] += -dvdx;
581: g[j+1][i+1] += dvdx + dvdy - exp(x[j+1][i+1]) * area * lambdad3;
583: floc += area * (0.5 * fquad - fexp * lambdad3);
585: }
586: }
588: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);
590: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
591: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
593: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
594: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
596: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
597: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
599: info = PetscLogFlops((xe-xs) * (ye-ys) * 55 + 1); CHKERRQ(info);
600: return 0;
602: } /* WholeCombFunctionGradient */
605: /*------- USER-DEFINED: routine to evaluate the Hessian
606: at the whole grid -------*/
609: /*
610: WholeCombHessian - Evaluates Hessian over the whole grid
612: Input:
613: daapplication - TAO application object
614: da - distributed array
615: X - the current point, at which the function and gradient are evaluated
616: ptr - user-defined application context
618: Output:
619: H - Hessian at X
620: */
621: static int WholeCombHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {
623: AppCtx *user = (AppCtx*)ptr;
624: Vec localX;
625: int info;
626: PetscInt i, j, ind[4];
627: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
628: double smallH[4][4];
629: double **x;
631: double hx, hy, lambdad3, area, dxdy, dydx;
632: double diagxy, dexp, bandxy, bandyx;
633: PetscTruth assembled;
636: hx = user->hx;
637: hy = user->hy;
638: lambdad3 = user->param / 3.0;
639: area = user->area;
640: dxdy = hx/hy;
641: dydx = hy/hx;
642: diagxy = 0.5 * (dxdy + dydx);
643: bandxy = -0.5 * dxdy;
644: bandyx = -0.5 * dydx;
646: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
647: info = MatAssembled(H,&assembled); CHKERRQ(info);
648: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
650: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
651: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
653: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
655: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
656: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
658: xe = gxs + gxm - 1;
659: ye = gys + gym - 1;
660: for (j = ys; j < ye; j++) {
661: for (i = xs; i < xe; i++) {
663: /* Hessian contribution at 0,0 */
664: dexp = exp(x[j][i]) * area * lambdad3;
665: smallH[0][0] = diagxy - dexp;
666: smallH[0][1] = smallH[1][0] = bandyx;
667: smallH[0][2] = smallH[2][0] = bandxy;
668: smallH[0][3] = smallH[3][0] = 0.0;
670: /* Hessian contribution at 1,0 */
671: dexp = exp(x[j][i+1]) * area * 2.0 * lambdad3;
672: smallH[1][1] = diagxy - dexp;
673: smallH[1][2] = smallH[2][1] = 0.0;
674: smallH[1][3] = smallH[3][1] = bandxy;
676: /* Hessian contribution at 0,1 */
677: dexp = exp(x[j+1][i]) * area * 2.0 * lambdad3;
678: smallH[2][2] = diagxy - dexp;
679: smallH[2][3] = smallH[3][2] = bandyx;
681: /* Hessian contribution at 1,1 */
682: dexp = exp(x[j+1][i+1]) * area * lambdad3;
683: smallH[3][3] = diagxy - dexp;
685: ind[0] = (j-gys) * gxm + (i-gxs);
686: ind[1] = ind[0] + 1;
687: ind[2] = ind[0] + gxm;
688: ind[3] = ind[2] + 1;
689: info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);
691: }
692: }
694: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
696: info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
697: info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
698: info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);
700: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
702: info = PetscLogFlops((xe-xs) * (ye-ys) * 14 + 7); CHKERRQ(info);
703: return 0;
705: } /* WholeCombHessian */