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