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: */