Actual source code: daelement.c
1: #include "taodaapplication.h" /* taodaapplication.h */
2: #include "taoapp.h"
4: typedef struct {
5: /* Function Gradient Evaluation over single element */
6: int (*computeelementfunctiongradient)(PetscInt[2], PetscScalar[4],double*,PetscScalar[4],void*);
7: void *elementfgctx;
8: int elementfgflops;
9: } TaoDA2D1DOFFGCtx;
14: static int TaoDA2dLoopFunctionGradient(TAO_APPLICATION daapp, DA da, Vec X, double *f, Vec G, void * ctx) {
16: TaoDA2D1DOFFGCtx* fgctx = (TaoDA2D1DOFFGCtx*) ctx;
17: Vec localX, localG;
18: MPI_Comm comm;
19: int info;
20: PetscInt i, j;
21: PetscInt coor[2];
22: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
23: double floc = 0.0, smallf;
24: PetscScalar **x, **g;
25: PetscScalar zero = 0.0;
26: PetscScalar smallX[4], smallG[4];
29: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
30: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
31: info = VecSet(G, zero); CHKERRQ(info);
32: info = VecSet(localG, zero); CHKERRQ(info);
34: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
35: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
37: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
38: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
40: info = DAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL); CHKERRQ(info);
41: info = DAGetGhostCorners(da, &gxs, &gys, PETSC_NULL, &gxm, &gym, PETSC_NULL); CHKERRQ(info);
43: xe = gxs + gxm - 1;
44: ye = gys + gym - 1;
45: for (j = ys; j < ye; j++) {
46: for (i = xs; i < xe; i++) {
48: smallX[0] = x[j][i];
49: smallX[1] = x[j][i+1];
50: smallX[2] = x[j+1][i];
51: smallX[3] = x[j+1][i+1];
52: coor[0] = i; coor[1] = j;
54: info = fgctx->computeelementfunctiongradient(coor,smallX,&smallf,smallG,fgctx->elementfgctx);
56: floc += smallf;
58: g[j][i] += smallG[0];
59: g[j][i+1] += smallG[1];
60: g[j+1][i] += smallG[2];
61: g[j+1][i+1] += smallG[3];
62: }
63: }
65: info = PetscLogFlops((ye-ys)*(xe-xs)*(fgctx->elementfgflops + 5)); CHKERRQ(info);
67: info = PetscObjectGetComm((PetscObject)X,&comm); CHKERRQ(info);
68: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, comm); CHKERRQ(info);
70: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
71: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
73: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
74: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
76: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
77: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
79: return(0);
80: } /* TaoDA2dLoopFunctionGradient */
84: /*@C
85: DAAppSetElementObjectiveAndGradientRoutine - Set routine that evaluates the
86: local part of a function on a 2-dimensional DA with 1 degree of freedom.
88: Collective on TAO_APPLICATION
90: Input Parameters:
91: + daapp - the TAO_APPLICATION solver context
92: . funcgrad - local function gradient routine
93: . flops - the number of flops done performed in the funcgrad routine
94: - fgctx - [optional] user-defined context for private data for the evaluation.
96: Calling sequence of funcgrad:
97: $ int funcgrad(int coordinates[2], PetscScalar x[4], double *f, PetscScalar g[4], void* ctx)
99: + coord - the global coordinates [i j] in each direction of the DA
100: . x - the variables on the DA ( da[j][i], da[j][j+1], da[j+1][i], da[j+1][i+1] ) (bottom left, bottom right, top left, top right)
101: . f - the local function value
102: . g - the gradient of this local function for with respect to each variable
103: - ctx - user defined context
105: Fortran Note:
106: If your Fortran compiler does not recognize symbols over 31 characters in length, then
107: use the identical routine with the shortened name DAAppSetElementObjectiveAndGrad()
108:
110: Level: intermediate
112: .keywords: DA, Object Function, Gradient
114: .seealso: DAAppSetObjectiveAndGradientRoutine();
115: @*/
116: int DAAppSetElementObjectiveAndGradientRoutine(TAO_APPLICATION daapplication, int (*funcgrad)(PetscInt[2],PetscScalar[4],double*,PetscScalar[4],void*),
117: PetscInt flops, void *ctx){
118: int info;
119: PetscInt n,i;
120: PetscInt dim,dof,s;
121: DAStencilType st;
122: TaoDA2D1DOFFGCtx *fgctx;
123: DA da;
126: info=DAAppGetNumberOfDAGrids(daapplication,&n); CHKERRQ(info);
127: for (i=0;i<n;i++){
128: info = DAAppGetDA(daapplication, i, &da); CHKERRQ(info);
129: info = DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,&s,0,&st); CHKERRQ(info);
130: if (dim!=2){
131: SETERRQ(1,"TAO DA ERROR: DA must have dimension of 2");}
132: if (dof!=1){
133: SETERRQ(1,"TAO DA ERROR: DA must have exactly 1 degree of freedom per node");}
134: if (s!=1){
135: SETERRQ(1,"TAO DA ERROR: DA stencil width must equal 1"); }
136: if (st!=DA_STENCIL_BOX){
137: SETERRQ(1,"TAO DA ERROR: DA stencil must be DA_STENCIL_BOX");}
138: }
139: PetscNew(TaoDA2D1DOFFGCtx,&fgctx);
140: fgctx->computeelementfunctiongradient=funcgrad;
141: fgctx->elementfgctx=ctx;
142: fgctx->elementfgflops = flops;
143: info = DAAppSetObjectiveAndGradientRoutine(daapplication, TaoDA2dLoopFunctionGradient, (void*)fgctx); CHKERRQ(info);
144: info = TaoAppSetDestroyRoutine(daapplication,TaoApplicationFreeMemory, (void*)fgctx); CHKERRQ(info);
145: info = PetscInfo(daapplication,"Set objective function pointer for TAO_APPLICATION object.\n"); CHKERRQ(info);
146: return(0);
147: }
150: typedef struct {
151: /* Hesian over single element */
152: int (*computeelementhessian)(PetscInt[2], PetscScalar[4],PetscScalar[4][4],void*);
153: void *elementhctx;
154: PetscInt elementhflops;
155: } TaoDA2D1DOFHCtx;
160: static int TaoDA2dLoopHessian(TAO_APPLICATION daapp, DA da, Vec X, Mat A, void* ctx) {
162: TaoDA2D1DOFHCtx* hctx = (TaoDA2D1DOFHCtx*)ctx;
163: Vec localX;
164: int info;
165: PetscInt i,j,coor[2];
166: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye,ind[4];
167: PetscScalar smallX[4];
168: PetscScalar smallH[4][4];
169: PetscScalar **x;
170: PetscTruth assembled;
174: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
175: info = MatAssembled(A,&assembled); CHKERRQ(info);
176: if (assembled) {
177: info = MatZeroEntries(A); CHKERRQ(info);
178: }
180: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
181: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
183: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
185: info = DAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL); CHKERRQ(info);
186: info = DAGetGhostCorners(da, &gxs, &gys, PETSC_NULL, &gxm, &gym, PETSC_NULL); CHKERRQ(info);
188: xe = gxs + gxm - 1;
189: ye = gys + gym - 1;
190: for (j = ys; j < ye; j++) {
191: for (i = xs; i < xe; i++) {
193: smallX[0] = x[j][i];
194: smallX[1] = x[j][i+1];
195: smallX[2] = x[j+1][i];
196: smallX[3] = x[j+1][i+1];
197: coor[0] = i; coor[1] = j;
199: info = hctx->computeelementhessian(coor,smallX,smallH,hctx->elementhctx); CHKERRQ(info);
201: ind[0] = (j-gys) * gxm + (i-gxs);
202: ind[1] = ind[0] + 1;
203: ind[2] = ind[0] + gxm;
204: ind[3] = ind[2] + 1;
205: info = MatSetValuesLocal(A,4,ind,4,ind,(PetscScalar*)smallH[0],ADD_VALUES); CHKERRQ(info);
206: }
207: }
209: info = PetscLogFlops((ye-ys)*(xe-xs)*hctx->elementhflops); CHKERRQ(info);
211: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
213: info = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
214: info = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
215: info = MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);
218: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
220: return(0);
221: } /* Compute Hessian */
226: /*@C
227: DAAppSetElementHessianRoutine - Set routine that evaluates the
228: local part of a Hessian matrix on a 2-dimensional DA with 1 degree of freedom.
230: Collective on TAO_APPLICATION
232: Input Parameters:
233: + daapp - the TAO_APPLICATION solver context
234: . hess - local function gradient routine
235: . flops - the number of flops done performed in the hess routine
236: - fgctx - [optional] user-defined context for private data for the evaluation.
238: Calling sequence of funcgrad:
239: $ int hess(int coordinates[2], PetscScalar x[4], PetscScalar h[4][4], void* ctx)
241: + coord - the global coordinates [i j] in each direction of the DA
242: . x - the variables on the DA ( da[j][i], da[j][j+1], da[j+1][i], da[j+1][i+1] ) (bottom left, bottom right, top left, top right)
243: . h - the hessian of this local function for with respect to each variable
244: - ctx - user defined context
246: Level: intermediate
248: .keywords: DA, gradient
250: .seealso: DAAppSetHessianRoutine()
251: @*/
252: int DAAppSetElementHessianRoutine(TAO_APPLICATION daapplication, int (*hess)(PetscInt[2],PetscScalar[4],PetscScalar[4][4],void*), PetscInt flops, void *ctx){
253: int info;
254: PetscInt n,i;
255: PetscInt dim,dof,s;
256: DAStencilType st;
257: TaoDA2D1DOFHCtx *hctx;
258: DA da;
261: info=DAAppGetNumberOfDAGrids(daapplication,&n); CHKERRQ(info);
262: for (i=0;i<n;i++){
263: info = DAAppGetDA(daapplication, i, &da); CHKERRQ(info);
264: info = DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,&s,0,&st); CHKERRQ(info);
265: if (dim!=2){
266: SETERRQ(1,"TAO DA ERROR: DA must have dimension of 2");}
267: if (dof!=1){
268: SETERRQ(1,"TAO DA ERROR: DA must have exactly 1 degree of freedom per node");}
269: if (s!=1){
270: SETERRQ(1,"TAO DA ERROR: DA stencil width must equal 1"); }
271: if (st!=DA_STENCIL_BOX){
272: SETERRQ(1,"TAO DA ERROR: DA stencil must be DA_STENCIL_BOX");}
273: }
274: PetscNew(TaoDA2D1DOFHCtx,&hctx);
275: hctx->computeelementhessian=hess;
276: hctx->elementhctx=ctx;
277: hctx->elementhflops = flops;
278: info = DAAppSetHessianRoutine(daapplication, TaoDA2dLoopHessian, (void*)hctx);
279: info = TaoAppSetDestroyRoutine(daapplication,TaoApplicationFreeMemory, (void*)hctx); CHKERRQ(info);
280: info = PetscInfo(daapplication,"Set Hessian for TAO_APPLICATION object and create matrix.\n"); CHKERRQ(info);
281: return(0);
282: }
284: /*
285: typedef struct {
286: / * Function Evaluation over entire single element * /
287: int (*computeelementfunction)(int[2], double*,int, double*,void*);
288: int (*computeelementgradient)(int[2], double*,int, double*,void*);
289: void *elementfctx;
290: void *elementgctx;
291: int elementfflops;
292: int elementgflops;
293: } TaoPetscFDHessianAppCtx;
294: */
295: /*
298: int DAAppSetElementFunction(TAO_APPLICATION daapplication, int (*func)(int,int,double*,int,double*,void*), int flops, void *ctx){
299: int info;
301: daapplication->computeelementfunction=func;
302: daapplication->usrfctx=ctx;
303: info = PetscLogInfo((daapplication->grid[0].da,"Set objective function pointer for TAO_APPLICATION object.\n")); CHKERRQ(info);
304: return(0);
305: } */
309: /*
313: int DAAppSetElementGradient(TAO_APPLICATION daapplication, int (*grad)(int,int,double*,int,double*,void*), int flops, void *ctx){
314: int info;
316: daapplication->computeelementgradient=grad;
317: daapplication->usrgctx=ctx;
318: info = PetscLogInfo((daapplication->grid[0].da,"Set gradient function pointer for TAO_APPLICATION object.\n")); CHKERRQ(info);
319: return(0);
320: }
321: */