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)(int[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, i, j, coor[2];
 20:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
 21:   double floc = 0.0, smallf;
 22:   PetscScalar **x, **g;
 23:   PetscScalar zero = 0.0;
 24:   PetscScalar smallX[4], smallG[4];

 27:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
 28:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
 29:   info = VecSet(G, zero); CHKERRQ(info);
 30:   info = VecSet(localG, zero); CHKERRQ(info);

 32:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
 33:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

 35:   info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
 36:   info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);

 38:   info = DAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL); CHKERRQ(info);
 39:   info = DAGetGhostCorners(da, &gxs, &gys, PETSC_NULL, &gxm, &gym, PETSC_NULL); CHKERRQ(info);

 41:   xe = gxs + gxm - 1;
 42:   ye = gys + gym - 1;
 43:   for (j = ys; j < ye; j++) {
 44:     for (i = xs; i < xe; i++) {

 46:         smallX[0] = x[j][i];
 47:         smallX[1] = x[j][i+1];
 48:         smallX[2] = x[j+1][i];
 49:         smallX[3] = x[j+1][i+1];
 50:         coor[0] = i; coor[1] = j;

 52:         info = fgctx->computeelementfunctiongradient(coor,smallX,&smallf,smallG,fgctx->elementfgctx);

 54:         floc += smallf;

 56:         g[j][i] += smallG[0];
 57:         g[j][i+1] += smallG[1];
 58:         g[j+1][i] += smallG[2];
 59:         g[j+1][i+1] += smallG[3];
 60:     }
 61:   }

 63:   info = PetscLogFlops((ye-ys)*(xe-xs)*(fgctx->elementfgflops + 5)); CHKERRQ(info);

 65:   info = PetscObjectGetComm((PetscObject)X,&comm); CHKERRQ(info);
 66:   info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, comm); CHKERRQ(info);

 68:   info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
 69:   info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);

 71:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
 72:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

 74:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
 75:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

 77:   return(0);
 78: } /* TaoDA2dLoopFunctionGradient  */

 82: /*@C
 83:    DAAppSetElementObjectiveAndGradientRoutine - Set routine that evaluates the
 84:    local part of a function on a 2-dimensional DA with 1 degree of freedom. 

 86:    Collective on TAO_APPLICATION

 88:    Input Parameters:
 89: +  daapp - the TAO_APPLICATION solver context
 90: .  funcgrad - local function gradient routine
 91: .  flops - the number of flops done performed in the funcgrad routine
 92: -  fgctx - [optional] user-defined context for private data for the evaluation.

 94:    Calling sequence of funcgrad:
 95: $     int funcgrad(int coordinates[2], PetscScalar x[4], double *f, PetscScalar g[4], void* ctx)

 97: +    coord - the global coordinates [i j] in each direction of the DA
 98: .    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)
 99: .    f - the local function value
100: .    g - the gradient of this local function for with respect to each variable
101: -    ctx - user defined context

103:    Fortran Note:
104:    If your Fortran compiler does not recognize symbols over 31 characters in length, then
105:    use the identical routine with the shortened name DAAppSetElementObjectiveAndGrad()
106:    

108:    Level: intermediate

110: .keywords: DA, Object Function, Gradient

112: .seealso: DAAppSetObjectiveAndGradientRoutine();
113: @*/
114: int DAAppSetElementObjectiveAndGradientRoutine(TAO_APPLICATION daapplication, int (*funcgrad)(int[2],PetscScalar[4],double*,PetscScalar[4],void*),
115:                                                int flops, void *ctx){
116:   int i,n,info;
117:   int dim,dof,s;
118:   DAStencilType st;
119:   TaoDA2D1DOFFGCtx *fgctx;
120:   DA da;

123:   info=DAAppGetNumberOfDAGrids(daapplication,&n); CHKERRQ(info);
124:   for (i=0;i<n;i++){
125:     info = DAAppGetDA(daapplication, i, &da); CHKERRQ(info);
126:     info = DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,&s,0,&st); CHKERRQ(info);
127:     if (dim!=2){
128:       SETERRQ(1,"TAO DA ERROR: DA must have dimension of 2");}
129:     if (dof!=1){
130:       SETERRQ(1,"TAO DA ERROR: DA must have exactly 1 degree of freedom per node");}
131:     if (s!=1){
132:       SETERRQ(1,"TAO DA ERROR: DA stencil width must equal 1"); }
133:     if (st!=DA_STENCIL_BOX){
134:       SETERRQ(1,"TAO DA ERROR: DA stencil must be DA_STENCIL_BOX");}
135:   }
136:   PetscNew(TaoDA2D1DOFFGCtx,&fgctx);
137:   fgctx->computeelementfunctiongradient=funcgrad;
138:   fgctx->elementfgctx=ctx;
139:   fgctx->elementfgflops = flops;
140:   info = DAAppSetObjectiveAndGradientRoutine(daapplication, TaoDA2dLoopFunctionGradient, (void*)fgctx); CHKERRQ(info);
141:   info = TaoAppSetDestroyRoutine(daapplication,TaoApplicationFreeMemory, (void*)fgctx); CHKERRQ(info);
142:   info = PetscLogInfo((daapplication,"Set objective function pointer for TAO_APPLICATION object.\n")); CHKERRQ(info);
143:   return(0);
144: }


147: typedef struct {
148:   /* Hesian over single element  */
149:   int  (*computeelementhessian)(int[2], PetscScalar[4],PetscScalar[4][4],void*);
150:   void *elementhctx;
151:   int elementhflops;
152: } TaoDA2D1DOFHCtx;


157: static int TaoDA2dLoopHessian(TAO_APPLICATION daapp, DA da, Vec X, Mat A, void* ctx) {

159:   TaoDA2D1DOFHCtx* hctx = (TaoDA2D1DOFHCtx*)ctx;
160:   Vec localX;
161:   int info, i, j, ind[4], coor[2];
162:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
163:   PetscScalar smallX[4];
164:   PetscScalar smallH[4][4];
165:   PetscScalar **x;
166:   PetscTruth assembled;


170:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
171:   info = MatAssembled(A,&assembled); CHKERRQ(info);
172:   if (assembled) {
173:     info = MatZeroEntries(A); CHKERRQ(info);
174:   }

176:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
177:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

179:   info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);

181:   info = DAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL); CHKERRQ(info);
182:   info = DAGetGhostCorners(da, &gxs, &gys, PETSC_NULL, &gxm, &gym, PETSC_NULL); CHKERRQ(info);

184:   xe = gxs + gxm - 1;
185:   ye = gys + gym - 1;
186:   for (j = ys; j < ye; j++) {
187:     for (i = xs; i < xe; i++) {

189:       smallX[0] = x[j][i];
190:       smallX[1] = x[j][i+1];
191:       smallX[2] = x[j+1][i];
192:       smallX[3] = x[j+1][i+1];
193:       coor[0] = i; coor[1] = j;

195:       info = hctx->computeelementhessian(coor,smallX,smallH,hctx->elementhctx); CHKERRQ(info);

197:       ind[0] = (j-gys) * gxm + (i-gxs);
198:       ind[1] = ind[0] + 1;
199:       ind[2] = ind[0] + gxm;
200:       ind[3] = ind[2] + 1;
201:       info = MatSetValuesLocal(A,4,ind,4,ind,(PetscScalar*)smallH[0],ADD_VALUES); CHKERRQ(info);
202:     }
203:   }

205:   info = PetscLogFlops((ye-ys)*(xe-xs)*hctx->elementhflops); CHKERRQ(info);

207:   info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);

209:   info = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
210:   info = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
211:   info = MatSetOption(A, MAT_SYMMETRIC); CHKERRQ(info);

213:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);

215:   return(0);
216: } /* Compute Hessian */


221: /*@C
222:    DAAppSetElementHessianRoutine - Set routine that evaluates the
223:    local part of a Hessian matrix on a 2-dimensional DA with 1 degree of freedom. 

225:    Collective on TAO_APPLICATION

227:    Input Parameters:
228: +  daapp - the TAO_APPLICATION solver context
229: .  hess - local function gradient routine
230: .  flops - the number of flops done performed in the hess routine
231: -  fgctx - [optional] user-defined context for private data for the evaluation.

233:    Calling sequence of funcgrad:
234: $     int hess(int coordinates[2], PetscScalar x[4], PetscScalar h[4][4], void* ctx)

236: +    coord - the global coordinates [i j] in each direction of the DA
237: .    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)
238: .    h - the hessian of this local function for with respect to each variable
239: -    ctx - user defined context

241:    Level: intermediate

243: .keywords: DA, gradient

245: .seealso: DAAppSetHessianRoutine()
246: @*/
247: int DAAppSetElementHessianRoutine(TAO_APPLICATION daapplication, int (*hess)(int[2],PetscScalar[4],PetscScalar[4][4],void*), int flops, void *ctx){
248:   int i,n,info;
249:   int dim,dof,s;
250:   DAStencilType st;
251:   TaoDA2D1DOFHCtx *hctx;
252:   DA da;

255:   info=DAAppGetNumberOfDAGrids(daapplication,&n); CHKERRQ(info);
256:   for (i=0;i<n;i++){
257:     info = DAAppGetDA(daapplication, i, &da); CHKERRQ(info);
258:     info = DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,&s,0,&st); CHKERRQ(info);
259:     if (dim!=2){
260:       SETERRQ(1,"TAO DA ERROR: DA must have dimension of 2");}
261:     if (dof!=1){
262:       SETERRQ(1,"TAO DA ERROR: DA must have exactly 1 degree of freedom per node");}
263:     if (s!=1){
264:       SETERRQ(1,"TAO DA ERROR: DA stencil width must equal 1"); }
265:     if (st!=DA_STENCIL_BOX){
266:       SETERRQ(1,"TAO DA ERROR: DA stencil must be DA_STENCIL_BOX");}
267:   }
268:   PetscNew(TaoDA2D1DOFHCtx,&hctx);
269:   hctx->computeelementhessian=hess;
270:   hctx->elementhctx=ctx;
271:   hctx->elementhflops = flops;
272:   info = DAAppSetHessianRoutine(daapplication, TaoDA2dLoopHessian,  (void*)hctx);
273:   info = TaoAppSetDestroyRoutine(daapplication,TaoApplicationFreeMemory, (void*)hctx); CHKERRQ(info);
274:   info = PetscLogInfo((daapplication,"Set Hessian for TAO_APPLICATION object and create matrix.\n")); CHKERRQ(info);
275:   return(0);
276: }

278: /*
279: typedef struct {
280:   / * Function Evaluation over entire single element  * /
281:   int  (*computeelementfunction)(int[2], double*,int, double*,void*); 
282:   int  (*computeelementgradient)(int[2], double*,int, double*,void*); 
283:   void *elementfctx;
284:   void *elementgctx;
285:   int elementfflops;
286:   int elementgflops;
287: } TaoPetscFDHessianAppCtx;
288: */
289: /*
292: int DAAppSetElementFunction(TAO_APPLICATION daapplication, int (*func)(int,int,double*,int,double*,void*), int flops, void *ctx){
293:   int info;
295:   daapplication->computeelementfunction=func;
296:   daapplication->usrfctx=ctx;
297:   info = PetscLogInfo((daapplication->grid[0].da,"Set objective function pointer for TAO_APPLICATION object.\n")); CHKERRQ(info);
298:   return(0);
299: } */



303: /*

307: int DAAppSetElementGradient(TAO_APPLICATION daapplication, int (*grad)(int,int,double*,int,double*,void*), int flops, void *ctx){
308:   int info;
310:   daapplication->computeelementgradient=grad;
311:   daapplication->usrgctx=ctx;
312:   info = PetscLogInfo((daapplication->grid[0].da,"Set gradient function pointer for TAO_APPLICATION object.\n")); CHKERRQ(info);
313:   return(0);
314: }
315: */